triplets.c
#include <R.h>
#include <math.h>
#include "methas.h"
#include "dist2.h"
/* Conditional intensity computation for Triplets process */
/* Format for storage of parameters and precomputed/auxiliary data */
typedef struct Triplets {
double gamma;
double r;
double loggamma;
double r2;
double *period;
int hard;
int per;
int *neighbour; /* scratch list of neighbours of current point */
int Nmax; /* length of scratch space allocated */
} Triplets;
/* initialiser function */
Cdata *tripletsinit(state, model, algo)
State state;
Model model;
Algor algo;
{
/* create storage for model parameters */
Triplets *triplets;
triplets = (Triplets *) R_alloc(1, sizeof(Triplets));
/* create scratch space */
triplets->Nmax = 1024;
triplets->neighbour = (int *) R_alloc(1024, sizeof(int));
/* Interpret model parameters*/
triplets->gamma = model.ipar[0];
triplets->r = model.ipar[1]; /* No longer passed as r^2 */
triplets->r2 = triplets->r * triplets->r;
triplets->period = model.period;
#ifdef MHDEBUG
Rprintf("Initialising Triplets gamma=%lf, r=%lf\n",
triplets->gamma, triplets->r);
#endif
/* is the model numerically equivalent to hard core ? */
triplets->hard = (triplets->gamma < DOUBLE_EPS);
triplets->loggamma = (triplets->hard) ? 0 : log(triplets->gamma);
/* periodic boundary conditions? */
triplets->per = (model.period[0] > 0.0);
return((Cdata *) triplets);
}
/* conditional intensity evaluator */
double tripletscif(prop, state, cdata)
Propo prop;
State state;
Cdata *cdata;
{
int npts, kount, ix, j, k, nj, nk, N, Nmax, Nmore, N1;
int *neighbour;
double *x, *y;
double u, v;
double r2, d2, cifval;
Triplets *triplets;
triplets = (Triplets *) cdata;
r2 = triplets->r2;
u = prop.u;
v = prop.v;
ix = prop.ix;
x = state.x;
y = state.y;
npts = state.npts;
if(npts == 0)
return((double) 1.0);
neighbour = triplets->neighbour;
Nmax = triplets->Nmax;
N = 0;
/* compile list of neighbours */
for(j=0; j < npts; j++) {
if(j != ix) {
d2 = dist2either(u,v,x[j],y[j],triplets->period);
if(d2 < r2) {
/* add j to list of neighbours of current point */
if(N >= Nmax) {
/* storage space overflow: reallocate */
Nmore = 2 * Nmax;
triplets->neighbour = neighbour =
(int *) S_realloc((char *) triplets->neighbour,
Nmore, Nmax, sizeof(int));
triplets->Nmax = Nmax = Nmore;
}
neighbour[N] = j;
N++;
}
}
}
/* count r-close (ordered) pairs of neighbours */
kount = 0;
if(N > 1) {
N1 = N - 1;
for(j = 0; j < N1; j++) {
nj = neighbour[j];
for(k = j+1; k < N; k++) {
nk = neighbour[k];
if(nj != nk) {
d2 = dist2either(x[nj],y[nj],x[nk],y[nk],triplets->period);
if(d2 < r2) kount++;
}
}
}
}
if(triplets->hard) {
if(kount > 0) cifval = 0.0;
else cifval = 1.0;
} else cifval = exp((triplets->loggamma) * kount);
#ifdef MHDEBUG
Rprintf("triplet count=%d cif=%lf\n", kount, cifval);
#endif
return cifval;
}
Cifns TripletsCifns = { &tripletsinit, &tripletscif, (updafunptr) NULL, NO};