Revision 28f4ee9e406fce786bb126198f84ce6ec43cb9c1 authored by Adrian Baddeley on 10 February 2011, 00:00:00 UTC, committed by Gabor Csardi on 10 February 2011, 00:00:00 UTC
1 parent c4ff53c
strauss.c
#include <R.h>
#include <math.h>
#include "methas.h"
#include "dist2.h"
/* Conditional intensity computation for Strauss process */
/* Format for storage of parameters and precomputed/auxiliary data */
typedef struct Strauss {
double beta;
double gamma;
double r;
double loggamma;
double r2;
double *period;
int hard;
int per;
} Strauss;
/* initialiser function */
Cdata *straussinit(state, model, algo)
State state;
Model model;
Algor algo;
{
/* create storage for model parameters */
Strauss *strauss;
strauss = (Strauss *) R_alloc(1, sizeof(Strauss));
/* Interpret model parameters*/
strauss->beta = model.par[0];
strauss->gamma = model.par[1];
strauss->r = model.par[2]; /* No longer passed as r^2 */
strauss->r2 = strauss->r * strauss->r;
strauss->period = model.period;
#ifdef MHDEBUG
Rprintf("Initialising Strauss beta=%lf, gamma=%lf, r=%lf\n",
strauss->beta, strauss->gamma, strauss->r);
#endif
/* is the model numerically equivalent to hard core ? */
strauss->hard = (strauss->gamma < DOUBLE_EPS);
strauss->loggamma = (strauss->hard) ? 0 : log(strauss->gamma);
/* periodic boundary conditions? */
strauss->per = (model.period[0] > 0.0);
return((Cdata *) strauss);
}
/* conditional intensity evaluator */
double strausscif(prop, state, cdata)
Propo prop;
State state;
Cdata *cdata;
{
int npts, kount, ix, ixp1, j;
double *period, *x, *y;
double u, v;
double r2, d2, a, cifval;
Strauss *strauss;
strauss = (Strauss *) cdata;
period = strauss->period;
r2 = strauss->r2;
u = prop.u;
v = prop.v;
ix = prop.ix;
x = state.x;
y = state.y;
npts = state.npts;
if(npts == 0)
return(strauss->beta);
kount = 0;
ixp1 = ix+1;
/* If ix = NONE = -1, then ixp1 = 0 is correct */
if(strauss->per) { /* periodic distance */
if(ix > 0) {
for(j=0; j < ix; j++) {
d2 = dist2(u,v,x[j],y[j],strauss->period);
if(d2 < r2) kount = kount+1;
}
}
if(ixp1 < npts) {
for(j=ixp1; j<npts; j++) {
d2 = dist2(u,v,x[j],y[j],strauss->period);
if(d2 < r2) kount = kount+1;
}
}
}
else { /* Euclidean distance */
if(ix > 0) {
for(j=0; j < ix; j++) {
a = r2 - pow(u - x[j], 2);
if(a > 0) {
a -= pow(v - y[j], 2);
if(a > 0) kount = kount+1;
}
}
}
if(ixp1 < npts) {
for(j=ixp1; j<npts; j++) {
a = r2 - pow(u - x[j], 2);
if(a > 0) {
a -= pow(v - y[j], 2);
if(a > 0) kount = kount+1;
}
}
}
}
if(strauss->hard) {
if(kount > 0) cifval = 0.0;
else cifval = strauss->beta;
}
else cifval = (strauss->beta) * exp((strauss->loggamma) * kount);
return cifval;
}
Cifns StraussCifns = { &straussinit, &strausscif, (updafunptr) NULL, FALSE};
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...