#include <R.h>
#include <math.h>
#include "methas.h"
#include "dist2.h"
/* for debugging code, include #define DEBUG 1 */
/* Conditional intensity computation for Multitype Hardcore process */
/* NOTE: types (marks) are numbered from 0 to ntypes-1 */
/* Storage of parameters and precomputed/auxiliary data */
typedef struct MultiHard {
int ntypes;
double *hc; /* hc[i,j] = hc[j+ntypes*i] for i,j = 0... ntypes-1 */
double *hc2; /* squared radii */
double range2; /* square of interaction range */
double *period;
int per;
} MultiHard;
/* initialiser function */
Cdata *multihardinit(state, model, algo)
State state;
Model model;
Algor algo;
{
int i, j, ntypes, n2;
double h, h2, range2;
MultiHard *multihard;
multihard = (MultiHard *) R_alloc(1, sizeof(MultiHard));
multihard->ntypes = ntypes = model.ntypes;
n2 = ntypes * ntypes;
#ifdef DEBUG
Rprintf("initialising space for %d types\n", ntypes);
#endif
/* Allocate space for parameters */
multihard->hc = (double *) R_alloc((size_t) n2, sizeof(double));
/* Allocate space for transformed parameters */
multihard->hc2 = (double *) R_alloc((size_t) n2, sizeof(double));
/* Copy and process model parameters*/
range2 = 0.0;
for(i = 0; i < ntypes; i++) {
for(j = 0; j < ntypes; j++) {
h = model.ipar[i + j*ntypes];
h2 = h * h;
MAT(multihard->hc, i, j, ntypes) = h;
MAT(multihard->hc2, i, j, ntypes) = h2;
if(range2 > h2) range2 = h2;
}
}
multihard->range2 = range2;
/* periodic boundary conditions? */
multihard->period = model.period;
multihard->per = (model.period[0] > 0.0);
#ifdef DEBUG
Rprintf("end initialiser\n");
#endif
return((Cdata *) multihard);
}
/* conditional intensity evaluator */
double multihardcif(prop, state, cdata)
Propo prop;
State state;
Cdata *cdata;
{
int npts, ntypes, ix, ixp1, j, mrk, mrkj;
int *marks;
double *x, *y;
double u, v;
double d2, range2, cifval;
double *period;
MultiHard *multihard;
DECLARE_CLOSE_D2_VARS;
multihard = (MultiHard *) cdata;
range2 = multihard->range2;
period = multihard->period;
u = prop.u;
v = prop.v;
mrk = prop.mrk;
ix = prop.ix;
x = state.x;
y = state.y;
marks = state.marks;
npts = state.npts;
#ifdef DEBUG
Rprintf("computing cif: u=%lf, v=%lf, mrk=%d\n", u, v, mrk);
#endif
cifval = 1.0;
if(npts == 0)
return(cifval);
ntypes = multihard->ntypes;
#ifdef DEBUG
Rprintf("scanning data\n");
#endif
ixp1 = ix+1;
/* If ix = NONE = -1, then ixp1 = 0 is correct */
if(multihard->per) { /* periodic distance */
if(ix > 0) {
for(j=0; j < ix; j++) {
if(CLOSE_PERIODIC_D2(u,v,x[j],y[j],period,range2,d2)) {
mrkj = marks[j];
if(d2 < MAT(multihard->hc2, mrk, mrkj, ntypes)) {
cifval = 0.0;
return(cifval);
}
}
}
}
if(ixp1 < npts) {
for(j=ixp1; j<npts; j++) {
if(CLOSE_PERIODIC_D2(u,v,x[j],y[j],period,range2,d2)) {
mrkj = marks[j];
if(d2 < MAT(multihard->hc2, mrk, mrkj, ntypes)) {
cifval = 0.0;
return(cifval);
}
}
}
}
} else { /* Euclidean distance */
if(ix > 0) {
for(j=0; j < ix; j++) {
if(CLOSE_D2(u, v, x[j], y[j], range2, d2)) {
mrkj = marks[j];
if(d2 < MAT(multihard->hc2, mrk, mrkj, ntypes)) {
cifval = 0.0;
return(cifval);
}
}
}
}
if(ixp1 < npts) {
for(j=ixp1; j<npts; j++) {
if(CLOSE_D2(u, v, x[j], y[j], range2, d2)) {
mrkj = marks[j];
if(d2 < MAT(multihard->hc2, mrk, mrkj, ntypes)) {
cifval = 0.0;
return(cifval);
}
}
}
}
}
#ifdef DEBUG
Rprintf("returning positive cif\n");
#endif
return cifval;
}
Cifns MultiHardCifns = { &multihardinit, &multihardcif, (updafunptr) NULL, YES};