/*
 *    slice.cpp
 *    Unbounded slice sampling with ease.
 *    $Id: slice.cpp,v 1.2 2015/10/13 06:21:43 daichi Exp $
 *
 */
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include "slice.h"
#include "random.h"

using namespace std;

typedef double (*likfun)(double x, void *args);

static double A = 100.0;

static double
expand (double p)	/* p -> x */
{
	return - A * log (1 / p - 1);
}

static double
shrink (double x)	/* x -> p */
{
	return 1 / (1 + exp(- x / A));
}

static double
expandp (double p)	/* p -> (x > 0) */
{
	return p / (1 - p);
}

static double
shrinkp (double x)	/* (x > 0) -> p */
{
	return  x / (1 + x);
}

double
slice_sample_positive (double x, likfun loglik, void *arg)
{
	double st = 0, ed = 1;
	double r, rnew, slice, newlik;
	int iter, maxiter = 1000;

	r = shrinkp (x);
	slice = (*loglik)(x, arg) - 2 * log (1 - r) + log (RANDOM);
	
	for (iter = 0; iter < maxiter; iter++)
	{
		rnew = unif (st, ed);
		newlik = (*loglik)(expandp(rnew), arg) - 2 * log (1 - rnew);

		if (newlik > slice)
			return expandp (rnew);
		else if (rnew > r)
			ed = rnew;
		else if (rnew < r)
			st = rnew;
		else
			return x;
	}
	fprintf(stderr, "slice_sample: max iteration %d reached.\n", maxiter);
	return x;
}

double
slice_sample (double x, likfun loglik, void *arg)
{
	double st = 0, ed = 1;
	double r, rnew, slice, newlik;
	int iter, maxiter = 1000;

	r = shrink (x);
	slice = (*loglik)(x, arg) - log (A * r * (1 - r)) + log (RANDOM);
	
	for (iter = 0; iter < maxiter; iter++)
	{
		rnew = unif (st, ed);
		newlik = (*loglik)(expand(rnew), arg)
					- log (A * rnew * (1 - rnew));

		if (newlik > slice)
			return expand (rnew);
		else if (rnew > r)
			ed = rnew;
		else if (rnew < r)
			st = rnew;
		else
			return x;
	}
	fprintf(stderr, "slice_sample: max iteration %d reached.\n", maxiter);
	return x;
}

#if 1

double
f (double x, void *arg)
{
	return - x * (x - 1) * (x - 2) * (x - 3.4);
}

double
f2 (double x, void *arg)
{
	return - (x - 1000) * (x - 1000);
}

int
main (int argc, char *argv[])
{
	int i, N = atoi(argv[1]);
	double x = 1;

	for (i = 0; i < N; i++)
	{
		// x = slice_sample_positive (x, f, NULL);
		x = slice_sample (x, f, NULL);
		printf("%.4lf\n", x);
	}

}

#endif