/* Code for generating perfectly random configurations from
 * the Strauss point process.  Uses read-once CFTP and Murdoch's
 * method, and the coupling method described in Wilson (1999).
 * (C) 1999, David B. Wilson
 */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
void choke ()
{
  fprintf(stderr,"Choke!  Bug or memory shortage.\n");
  exit(-1);
}

inline double uniform() {
  return (((((double)random()) / 2147483648e0) + (double)random()) / 2147483648e0);
}

/* Strauss point process parameters
 */
typedef struct {
  double width, height, lambda, radius, repulsion;
} Strauss;
/* Point set utilities.  For large data sets it would be good
 * to do geometric binning, but I don't have the patience to
 * write code for that just now.
 */
typedef struct {
  double x,y;
} Point;
typedef struct {
  int num, size;
  Point *point;
} Points;
Points *EmptyPointSet ()
{
  Points *points;
  points = (Points*)malloc(sizeof(Points));
  if (!points) choke();
  points->num = 0;
  points->size = 1;
  points->point = (Point*)malloc(sizeof(Point));
  if (!points->point) choke();
  return points;
}
void DestroyPointSet (Points *points)
{
  free(points->point);
  free(points);
}
void AddPoint (Points *points, Point point)
{
  if (points->num==points->size) {
    points->size *= 2;
    points->point = realloc(points->point,points->size*sizeof(Point));
    if (!points->point) choke();
  }
  if (points->num>=points->size) choke();
  points->point[points->num++] = point;
}
void MergePointSet (Points *pts, Points *newpts)
{
  int i;
  for (i=0; i<newpts->num; ++i)
    AddPoint(pts,newpts->point[i]);
}
Points *CopyPointSet (Points *pts)
{
  Points *cpy;
  cpy = EmptyPointSet();
  MergePointSet(cpy,pts);
  return cpy;
}
void RemovePoint (Points *points, int pt)
{
  if (pt<0 || pt>= points->num) choke();
  points->point[pt] = points->point[--points->num];
}
void RemoveRandomPoint (Points *points)
{
  RemovePoint(points,(int) (points->num * uniform()));
}
inline double square (double x) {return x*x;}
double NewPointLikelihood (Points *points, Point point, Strauss ptp)
{
  double r = square(ptp.radius);
  double x = point.x;
  double y = point.y;
  Point *pts = points->point;
  int close = 0;
  int i;
  for (i=(points->num)-1; i>=0; --i)
    close += (square(pts[i].x - x) + square(pts[i].y - y) < r);
  return pow(ptp.repulsion,(double)close);
}
int NumClosePoints (Points *points, Strauss ptp)
{
  double r = square(ptp.radius);
  int i,j, close;
  close = 0;
  for (i=0; i<points->num; ++i)
    for (j=0; j<i; ++j)
      if (square(points->point[i].x - points->point[j].x) +
	  square(points->point[i].y - points->point[j].y) < r)
	++close;
  return close;
}

/* Early routine for testing.
 * Evolves a point process configuration for continuous time t.
 */
void EvolveConfig (Points *points, Strauss ptp, double t)
{
  double birthrate = ptp.lambda * ptp.width * ptp.height;
  int deathrate;
  double clock=0;
  while (1) {
    deathrate = points->num;
    clock += -log(1-uniform())/(birthrate+deathrate);
    if (clock>t) return;
    if (uniform() < deathrate/(deathrate+birthrate)) {
      RemoveRandomPoint(points);
    } else {
      Point pt;
      pt.x = ptp.width * uniform();
      pt.y = ptp.height * uniform();
      if (uniform() < NewPointLikelihood(points,pt,ptp))
	AddPoint(points,pt);
    }
  }
}

/* Produces a Poisson random variable with parameter lambda.
 * If lambda is too large, the procedure may underflow.
 * Other numerical problems are conceivable, e.g. the numerical
 * sum of probabilities could be marginally less than 1, in which
 * case there's a chance of the procedure not terminating.
 */
int PoissonOld (double lambda)
{
  double prob, u;
  int i;
  u = uniform();
  prob = exp(-lambda);
  for (i=0; ; ) {
    u -= prob;
    if (u<0) return i;
    prob *= lambda/(++i);
    if (prob==0) choke(); /* check for underflow */
  }
}
/* Slow for large lambda, but acceptable for present purposes.
 */
int Poisson (double lambda)
{
  int i;
  for (i=0; ; ++i) {
    lambda += log(1-uniform());
    if (lambda<=0)  return i;
  }
}

/* Wilson (1999) suggests taking the proposal distribution to be a Poisson
 * process with density 2*K*lambda.  For the Strauss process K=1, so the
 * intensity ratio is 2.
 */
#define IntensityRatio ((double)2.0)
Points *MHproposal (Strauss ptp)
{
  int i;
  Points *points;
  points = EmptyPointSet();
  for (i=Poisson(IntensityRatio * ptp.lambda * ptp.width*ptp.height); i; --i) {
    Point pt;
    pt.x = ptp.width * uniform();
    pt.y = ptp.height * uniform();
    AddPoint(points,pt);
  }
  return points;
}

/* From Wilson (1999):
 *   Our representation for a set of configurations of the point process
 *   will consist of an integer $k$ and two finite sets of points $\Delta$
 *   and $L$.  The set represented by $(k,\Delta,L)$ consists of all
 *   configurations that contain each point in $L$, possibly some of the
 *   points in $\Delta$, and at most $k$ additional points that could be
 *   anywhere.
 * Here we have
 * Delta = Delt1 U Delt2
 * configuration = pnts U Delt1 U Lower
 *   (at any time either pnts or Lower is empty)
 * lower bound   = Lower
 * upper bound   = Lower U Delt1 U Delt2 U "extra" additional points
 *
 * When pnts is NULL, do a coalescence timing experiment, and store the
 * (continous) time to coalescence in "limit".
 * When pnts is not NULL, evolve the configuration given by pnts by doing
 * a Metropolis-Hastings update with an independence sampler (a la Murdoch)
 * and then run the birth-and-death Markov chain for an amount of continuous
 * time given by "limit", while simultaneously testing for coalescence.
 * See section 9 of Wilson (1999).
 */
CompositeMap (Strauss ptp, Points *pnts, double *limit)
{
  Points *proposal, *Delt1, *Delt2, *Lower;
  int extra;
  double birthrate = ptp.lambda * ptp.width * ptp.height;
  double clock, u,rate;

  fprintf(stderr,"."); fflush(stderr);
  Delt1 = EmptyPointSet();
  Delt2 = EmptyPointSet();
  Lower = EmptyPointSet();
  proposal = MHproposal(ptp);
  /* Compute the max # points in a configuration after MH update. */
  extra = ceil(proposal->num - NumClosePoints(proposal,ptp)*
	       log(ptp.repulsion)/log(IntensityRatio));
  /* Maybe the pnts configuration gets changed to the proposal. */
  if (pnts && uniform()<
      pow(ptp.repulsion,(double)(NumClosePoints(proposal,ptp)-
				 NumClosePoints(pnts,ptp)))*
      pow(IntensityRatio,(double)(pnts->num-proposal->num)))
    {free(pnts->point); *pnts = *proposal;}
  else
    free(proposal->point);
  free(proposal);
  if (pnts && pnts->num>extra) choke();  /* sanity check */
  clock = 0;
  while (extra) {
    rate = birthrate + Delt1->num + Delt2->num + extra;
    clock += -log(1-uniform())/rate;
    if (pnts && clock>*limit) {
      MergePointSet(pnts,Delt1);
      DestroyPointSet(Lower); DestroyPointSet(Delt2); DestroyPointSet(Delt1);
      fprintf(stderr,"\n"); fflush(stderr);
      return 0;
    }
    u = uniform();
    if (u<birthrate/rate) {
      Point pt; pt.x = ptp.width * uniform(); pt.y = ptp.height * uniform();
      if (pnts &&
	  uniform()<NewPointLikelihood(pnts ,pt,ptp) &&
	  uniform()<NewPointLikelihood(Delt1,pt,ptp))
	AddPoint(Delt1,pt);
      else
	AddPoint(Delt2,pt);
    } else if (u<(birthrate+Delt1->num)/rate)
      RemoveRandomPoint(Delt1);
    else if (u<(birthrate+Delt1->num+Delt2->num)/rate)
      RemoveRandomPoint(Delt2);
    else {
      if (pnts && uniform()<((double)pnts->num)/extra)
	RemoveRandomPoint(pnts);
      --extra;
    }
  }
  if (pnts && pnts->num) choke();  /* sanity check */
  fprintf(stderr,"."); fflush(stderr);
  while (1) {
    if (!pnts && !(Delt1->num+Delt2->num)) {
      *limit = clock;
      DestroyPointSet(Lower); DestroyPointSet(Delt2); DestroyPointSet(Delt1);
      fprintf(stderr,"\n"); fflush(stderr);
      return -1;
    }
    rate = birthrate + Delt1->num + Delt2->num + Lower->num;
    clock += -log(1-uniform())/rate;
    if (pnts && clock>*limit) {
      int coalesced = !(Delt1->num+Delt2->num);
      MergePointSet(pnts,Lower); MergePointSet(pnts,Delt1);
      DestroyPointSet(Lower); DestroyPointSet(Delt2); DestroyPointSet(Delt1);
      fprintf(stderr,"\n"); fflush(stderr);
      return coalesced;
    }
    u = uniform();
    if (u<birthrate/rate) {
      Point pt; pt.x = ptp.width * uniform(); pt.y = ptp.height * uniform();
      if (uniform()<NewPointLikelihood(Lower,pt,ptp))
	if (uniform()<NewPointLikelihood(Delt1,pt,ptp))
	  if (uniform()<NewPointLikelihood(Delt2,pt,ptp))
	    AddPoint(Lower,pt);
	  else
	    AddPoint(Delt1,pt);
	else
	  AddPoint(Delt2,pt);
    } else if (u<(birthrate+Delt1->num)/rate)
      RemoveRandomPoint(Delt1);
    else if (u<(birthrate+Delt1->num+Delt2->num)/rate)
      RemoveRandomPoint(Delt2);
    else
      RemoveRandomPoint(Lower);
  }
}
ApplyCompositeMap (Strauss ptp, Points *pnts)
{
  double clock;
  int coalescence;
  /* timing experiment */
  coalescence = CompositeMap(ptp,NULL,&clock);
  if (coalescence != -1) choke();
  /* actual map        */
  coalescence = CompositeMap(ptp,pnts,&clock);
  if (coalescence != 0 && coalescence != 1) choke();
  return coalescence;
}
Points *state;
void InitializeROCFTP (Strauss ptp)
{
  state = EmptyPointSet();
  while (!ApplyCompositeMap(ptp, state));
}
Points *NextSampleROCFTP (Strauss ptp)
{
  Points *oldstate;
  
  while (1) {
    oldstate = CopyPointSet(state);
    if (ApplyCompositeMap(ptp, state))
      return oldstate;
    DestroyPointSet(oldstate);
  }
}
     
ShowPts (Points *points, Strauss ptp)
{
  FILE *ps;
  double scalex,scaley,scale;
  double x1, x2, y1, y2;
  int i;
  ps = fopen("/tmp/pt.ps","w");
  scalex = 6.5*72/ptp.width;
  scaley = 9*72/ptp.height;
  scale = (scalex>scaley) ? scaley : scalex;
  x1 = 8.5*36-ptp.width*scale/2;
  x2 = 8.5*36+ptp.width*scale/2;
  y1 = 11*36-ptp.height*scale/2;
  y2 = 11*36+ptp.height*scale/2;
  fprintf(ps,
"%%!PS-Adobe-2.0 
%%%%Title: picture of Strauss point process
%%%%Creator: strauss (written by dbwilson)
%%%%BoundingBox: %lf %lf %lf %lf
%%%%Pages: 1
%%%%EndComments
%%%%Page: strauss 1
%% %lf X %lf @ intensity %lf, radius %lf, repulsion %lf
10 dict begin
gsave
%lf %lf translate
/sc {%lf mul} def
0 0 moveto %lf sc 0 lineto 0 %lf sc rlineto -%lf sc 0 rlineto closepath gsave stroke grestore clip newpath
/pt {exch sc exch sc 1 index 1 index 1 0 360 arc fill 0.5 %lf mul sc 0 360 arc stroke} def
",
	  x1,y1,x2,y2,
	  ptp.width,ptp.height,ptp.lambda,ptp.radius,ptp.repulsion,
	  x1,y1,
	  scale,
	  ptp.width,ptp.height,ptp.width,
	  ptp.radius);
  for (i=0; i<points->num; ++i)
    fprintf(ps,"%lf %lf pt\n",points->point[i].x,points->point[i].y);
  fprintf(ps,"
grestore
showpage
%%%%PageTrailer
end
%%%%Trailer
%%%%EOF
");
  fclose(ps);
  /*  system("ghostview -geometry -0+0 /tmp/pt.ps");*/
}

main (int argc, char *argv[])
{
  Strauss ptp;
  Points *points;
  int seed=1;
  if (argc != 6 && argc != 7) {
    fprintf(stderr,"%s width height lambda radius repulsion [seed]\n",argv[0]);
    exit(-1);
  }
  ptp.width = atof(argv[1]);
  ptp.height = atof(argv[2]);
  ptp.lambda = atof(argv[3]);
  ptp.radius = atof(argv[4]);
  ptp.repulsion = atof(argv[5]);
  if (argc==7) seed = atoi(argv[6]);

  if (ptp.width<0 || ptp.height<0 || ptp.lambda<0 || ptp.radius<0 ||
      ptp.repulsion<0 || ptp.repulsion>1) {
    fprintf(stderr,"width, height, lambda, and radius must be nonnegative.\n");
    fprintf(stderr,"repulsion must be between 0 and 1 (inclusive)");
    exit(-1);
  }
  fprintf(stderr,"%lfx%lf @ %lf (%lf) *%lf, seed %d\n",
	  ptp.width,
	  ptp.height,
	  ptp.lambda,
	  ptp.radius,
	  ptp.repulsion,
	  seed);
  srandom(seed);
  if (ptp.repulsion>0) {
    InitializeROCFTP(ptp);
    points = NextSampleROCFTP(ptp);
  } else {
    ptp.repulsion = 1e-6;
    InitializeROCFTP(ptp);
    do {
      fprintf(stderr,"attempting repulsion=1e-6 and rejection sampling\n");
      points = NextSampleROCFTP(ptp);
    } while (NumClosePoints(points,ptp)!=0);
    ptp.repulsion = 0;
  }
  ShowPts (points, ptp);
  DestroyPointSet(points);
  /*
  {
    int i;
    for (i=0; i<=10; ++i) {
      ptp.radius = i*0.1;
      srandom(1);
      EvolveConfig(points, ptp, 1.0);
      ShowPts(points, ptp);
    }
  }
  */
}