samer@0: /*
samer@0:  * Pseudorandom generators for SWI Prolog
samer@0:  * Samer Abdallah (2009)
samer@0:  * Some samplers adapted from code in Tom Minka's lightspeed Matlab toolbox.
samer@0:  * Other bits of code culled from SWI Prolog distribution random.c
samer@0:  *
samer@0:  * To do:
samer@0:  *    work out lower bound for skew stable distributions
samer@0:  * 	stream splitting with jumps
samer@0:  * 	reduce blob copying
samer@0:  * 	fast discrete distributions
samer@0:  * 	ziggurat for normals
samer@0:  */
samer@0: 
samer@0: #define _USE_MATH_DEFINES 1
samer@0: 
samer@0: #include <math.h>
samer@0: #include <float.h>
samer@0: #include "rndutils.h"
samer@0: 
samer@0: #if HAVE_GSL
samer@0: #include <gsl/gsl_sf_zeta.h>
samer@0: #else
samer@0: double gsl_sf_zeta(const double s) { return NAN; }
samer@0: double gsl_sf_hzeta(const double s, const double k) { return NAN; }
samer@0: #endif
samer@0: 
samer@0: #ifndef M_PI
samer@0: #define M_PI       3.14159265358979323846
samer@0: #endif
samer@0: 
samer@0: 
samer@0: // -----------------------------------------------------
samer@0: 
samer@0: // Returns a sample from Normal(0,1)
samer@0: double Normal(RndState *S)
samer@0: {
samer@0: 	double x=S->prev_normal;
samer@0: 
samer@0: 	if(!isnan(x)) {
samer@0: 		S->prev_normal=NAN;
samer@0: 	} else {
samer@0: 		double y,radius;
samer@0: 
samer@0: 		/* Generate a random point inside the unit circle */
samer@0: 		do {
samer@0: 			x = 2*Uniform(S)-1;
samer@0: 			y = 2*Uniform(S)-1;
samer@0: 			radius = (x*x)+(y*y);
samer@0: 		} while((radius >= 1.0) || (radius == 0.0));
samer@0: 		/* Box-Muller formula */
samer@0: 		radius = sqrt(-2*log(radius)/radius);
samer@0: 		x *= radius;
samer@0: 		y *= radius;
samer@0: 		S->prev_normal=y;
samer@0: 	}
samer@0: 	return x;
samer@0: }
samer@0: 
samer@0: /* Returns a sample from Gamma(a, 1).
samer@0:  * For Gamma(a,b), scale the result by b.
samer@0:  */
samer@0: double Gamma(RndState *S, double a)
samer@0: {
samer@0:   /* Algorithm:
samer@0:    * G. Marsaglia and W.W. Tsang, A simple method for generating gamma
samer@0:    * variables, ACM Transactions on Mathematical Software, Vol. 26, No. 3,
samer@0:    * Pages 363-372, September, 2000.
samer@0:    * http://portal.acm.org/citation.cfm?id=358414
samer@0:    */
samer@0:   double boost, d, c, v;
samer@0:   if(a < 1) {
samer@0:     /* boost using Marsaglia's (1961) method: gam(a) = gam(a+1)*U^(1/a) */
samer@0:     boost = exp(log(Uniform(S))/a);
samer@0:     a++;
samer@0:   } 
samer@0:   else boost = 1;
samer@0:   d = a-1.0/3; c = 1.0/sqrt(9*d);
samer@0:   while(1) {
samer@0:     double x,u;
samer@0:     do {
samer@0:       x = Normal(S);
samer@0:       v = 1+c*x;
samer@0:     } while(v <= 0);
samer@0:     v = v*v*v;
samer@0:     x = x*x;
samer@0:     u = Uniform(S);
samer@0:     if((u < 1-.0331*x*x) || 
samer@0:        (log(u) < 0.5*x + d*(1-v+log(v)))) break;
samer@0:   }
samer@0:   return( boost*d*v );
samer@0: }
samer@0: 
samer@0: /* Returns a sample from Dirichlet(n,a) where n is dimensionality */
samer@0: int Dirichlet(RndState *S, long n, double *a, double *x)
samer@0: {
samer@0: 	long i;
samer@0: 	double tot=0, z;
samer@0: 	for (i=0; i<n; i++) { z=Gamma(S,a[i]); tot+=z; x[i]=z; }
samer@0: 	for (i=0; i<n; i++) { x[i]/=tot; }
samer@0: 	return 1;
samer@0: }
samer@0: 
samer@0: /* Returns a sample from Beta(a,b) */
samer@0: double Beta(RndState *S, double a, double b)
samer@0: {
samer@0:   double g = Gamma(S,a);
samer@0:   return g/(g + Gamma(S,b));
samer@0: }
samer@0: 
samer@0: double Hyperbolic(RndState *S, double a)
samer@0: {
samer@0: 	return floor(pow(1-Uniform(S),-1/a));
samer@0: }
samer@0: 
samer@0: /* Sample from zeta distribution
samer@0:  *
samer@0:  * This is an inverse CDF method. It works by using  an 
samer@0:  * approximation of the Hurwitz Zeta function and then
samer@0:  * walking left or right until we get to the correct point.
samer@0:  * The approximation is good for large values so usually we
samer@0:  * only have to walk left or right a few steps for small 
samer@0:  * values.
samer@0:  */
samer@0: double Zeta(RndState *S, double s)
samer@0: {
samer@0: 	double v=(1-Uniform(S))*gsl_sf_zeta(s);
samer@0: 	double k0, k=round(pow(v*(s-1),1/(1-s))); // approx inverse of hzeta
samer@0: 	double zk0, zk=gsl_sf_hzeta(s,k);
samer@0: 
samer@0: 	if (zk>=v) {
samer@0: 		do {
samer@0: 			k0=k; zk0=zk;
samer@0: 			k=k0+1; zk=zk0-pow(k0,-s);
samer@0: 		} while (zk>=v && zk!=zk0);
samer@0: 		return k0;
samer@0: 	} else {
samer@0: 		do {
samer@0: 			k0=k; zk0=zk;
samer@0: 			k=k0-1; zk=zk0+pow(k,-s);
samer@0: 		} while (zk<v && zk!=zk0);
samer@0: 		return k;
samer@0: 	}
samer@0: }
samer@0: 
samer@0: /* Very fast binomial sampler. 
samer@0:  * Returns the number of successes out of n trials, with success probability p.
samer@0:  */
samer@0: int Binomial(RndState *S, double p, int n)
samer@0: {
samer@0:   int r = 0;
samer@0:   if(isnan(p)) return 0;
samer@0:   if(p < DBL_EPSILON) return 0;
samer@0:   if(p >= 1-DBL_EPSILON) return n;
samer@0:   if((p > 0.5) && (n < 15)) {
samer@0:     /* Coin flip method. This takes O(n) time. */
samer@0:     int i;
samer@0:     for(i=0;i<n;i++) {
samer@0:       if(Uniform(S) < p) r++;
samer@0:     }
samer@0:     return r;
samer@0:   }
samer@0:   if(n*p < 10) {
samer@0:     /* Waiting time method.  This takes O(np) time. */
samer@0:     double q = -log(1-p), e = -log(Uniform(S)), s;
samer@0:     r = n;
samer@0:     for(s = e/r; s <= q; s += e/r) {
samer@0:       r--;
samer@0:       if(r == 0) break;
samer@0:       e = -log(Uniform(S));
samer@0:     }
samer@0:     r = n-r;
samer@0:     return r;
samer@0:   }
samer@0:   if (1) {
samer@0:     /* Recursive method.  This makes O(log(log(n))) recursive calls. */
samer@0:     int i = (int)floor(p*(n+1));
samer@0:     double b = Beta(S,i, n+1-i);
samer@0:     if(b <= p) r = i + Binomial(S,(p-b)/(1-b), n-i);
samer@0:     else r = i - 1 - Binomial(S,(b-p)/b, i-1);
samer@0:     return r;
samer@0:   }
samer@0: }
samer@0: 
samer@0: double Exponential(RndState *S) { return -log(1-Uniform(S)); }
samer@0: 
samer@0: long Poisson(RndState *S, double lambda)
samer@0: {
samer@0: 	long r;
samer@0: 	if (lambda>=15) {
samer@0: 		double m=floor(lambda*7/8);
samer@0: 		double x=Gamma(S,m);
samer@0: 
samer@0: 		if (x>lambda) r=Binomial(S,lambda/x,m-1);
samer@0: 		else          r=m+Poisson(S,lambda-x);
samer@0: 	} else {
samer@0: 		double p;
samer@0: 		for (p=0, r=-1; p<lambda; p+=Exponential(S)) r++; 
samer@0: 	}
samer@0: 	return r;
samer@0: }
samer@0: 
samer@0: /* Returns a sample from stable distribution */
samer@0: double Stable(RndState *S, int param, double alpha, double beta)
samer@0: {
samer@0: 	double X, theta=M_PI*(Uniform(S)-0.5);
samer@0: 
samer@0: 	// These methods come from John Nolan's book about
samer@0: 	// stable distributions. They return standardised stable random 
samer@0: 	// numbers in the S(alpha,beta;0) parameterisation
samer@0: 
samer@0: 	if (beta==0) {
samer@0: 		if (alpha==1) {
samer@0: 			X=tan(theta);
samer@0: 		} else {
samer@0: 			double W=Exponential(S);
samer@0: 			X=sin(alpha*theta)/cos(theta)*pow(cos((alpha-1)*theta)/(W*cos(theta)),1/alpha-1);
samer@0: 		}
samer@0: 	} else {
samer@0: 		double W=Exponential(S);
samer@0: 		double zeta=beta*tan(alpha*M_PI/2);
samer@0: 		if (alpha==1) {
samer@0: 			double b=2*beta/M_PI;
samer@0: 			double bt=1+b*theta;
samer@0: 			X=bt*tan(theta) - b*log(W*cos(theta)/bt);
samer@0: 		} else {
samer@0: 			double ath = alpha*theta; 
samer@0: 			double a1th = theta-ath; 
samer@0: 			X = ((sin(ath) + zeta*cos(ath))/cos(theta))
samer@0: 				 * pow((cos(a1th) + zeta*sin(a1th))/(W*cos(theta)),1/alpha-1);
samer@0: 			if (param==0) X=X-zeta;
samer@0: 		}
samer@0: 	}
samer@0: 	return X;
samer@0: }
samer@0: 
samer@0: /* Returns a sample from discrete distribution */
samer@0: long Discrete(RndState *S, long n, double *p, double tot)
samer@0: {
samer@0: 	double u;
samer@0: 	long i;
samer@0: 
samer@0: 	for (i=0, u=tot*Uniform(S)-p[0]; u>=0; u-=p[++i]); 
samer@0: 	return i;
samer@0: }
samer@0: 
samer@0: double Uniform(RndState *S0)
samer@0: {
samer@0: 	return RngStream_Double(&S0->g,&S0->g);
samer@0: }
samer@0: 
samer@0: double Raw(RndState *S0)
samer@0: {
samer@0: 	return RngStream_Raw(&S0->g,&S0->g);
samer@0: }
samer@0: 
samer@0: void InitRndState(RndState *S)
samer@0: {
samer@0: 	RngStream_InitState(&S->g);
samer@0: 	S->prev_normal=NAN;
samer@0: }
samer@0: 
samer@0: 
samer@0: int spawn_gen(RndState *S0, RndState *S1)
samer@0: {
samer@0: 	unsigned seeds[6];
samer@0: 	int i;
samer@0: 
samer@0: 	for (i=0; i<6; i++) seeds[i]=(unsigned)RngStream_Raw(&S0->g,&S0->g);
samer@0: 	RngStream_SetState(&S1->g,seeds);
samer@0: 	for (i=0; i<3; i++) RngStream_Raw(&S1->g,&S1->g);
samer@0: 	S1->prev_normal=NAN;
samer@0: 	return 1;
samer@0: }
samer@0: 
samer@0: int randomise_from(FILE *p, RndState *S)
samer@0: {
samer@0: 	unsigned seeds[6];
samer@0: 	size_t r=fread((void *)seeds,sizeof(unsigned),6,p);
samer@0: 
samer@0: 	if (r==6) {
samer@0: 		int i;
samer@0: 		if (RngStream_SetState(&S->g,seeds)<0) return 0;
samer@0: 		for (i=0; i<3; i++) RngStream_Raw(&S->g,&S->g);
samer@0: 		S->prev_normal=NAN;
samer@0: 		return 1;
samer@0: 	} else return 0;
samer@0: }
samer@0: