https://github.com/cran/spatstat
Revision 3aca716ce2576a0dab83f08052acd47afed8ee6a authored by Adrian Baddeley on 29 February 2012, 00:00:00 UTC, committed by Gabor Csardi on 29 February 2012, 00:00:00 UTC
1 parent e567587
Tip revision: 3aca716ce2576a0dab83f08052acd47afed8ee6a authored by Adrian Baddeley on 29 February 2012, 00:00:00 UTC
version 1.25-4
version 1.25-4
Tip revision: 3aca716
geyer.c
#include <R.h>
#include <math.h>
#include <stdlib.h>
#include "methas.h"
#include "dist2.h"
void fexitc(const char *msg);
/*
Conditional intensity function for a Geyer saturation process.
*/
typedef struct Geyer {
/* model parameters */
double gamma;
double r;
double s;
/* transformations of the parameters */
double r2;
double loggamma;
int hard;
/* periodic distance */
double *period;
int per;
/* auxiliary counts */
int *aux;
} Geyer;
Cdata *geyerinit(state, model, algo)
State state;
Model model;
Algor algo;
{
int i, j;
double d2;
Geyer *geyer;
geyer = (Geyer *) R_alloc(1, sizeof(Geyer));
/* Interpret model parameters*/
geyer->gamma = model.ipar[0];
geyer->r = model.ipar[1]; /* not squared any more */
geyer->s = model.ipar[2];
geyer->r2 = geyer->r * geyer->r;
#ifdef MHDEBUG
Rprintf("Initialising Geyer gamma=%lf, r=%lf, sat=%lf\n",
geyer->gamma, geyer->r, geyer->s);
#endif
/* is the model numerically equivalent to hard core ? */
geyer->hard = (geyer->gamma < DOUBLE_EPS);
geyer->loggamma = (geyer->hard) ? 0 : log(geyer->gamma);
/* periodic boundary conditions? */
geyer->period = model.period;
geyer->per = (model.period[0] > 0.0);
/* allocate storage for auxiliary counts */
geyer->aux = (int *) R_alloc((size_t) state.npmax, sizeof(int));
/* Initialise auxiliary counts */
for(i = 0; i < state.npmax; i++)
geyer->aux[i] = 0;
for(i = 0; i < state.npts; i++) {
for(j = 0; j < state.npts; j++) {
d2 = dist2either(state.x[i], state.y[i], state.x[j], state.y[j],
geyer->period);
if(d2 < geyer->r2)
geyer->aux[i] += 1;
}
}
return((Cdata *) geyer);
}
double geyercif(prop, state, cdata)
Propo prop;
State state;
Cdata *cdata;
{
int ix, j, npts, tee;
double u, v, d2, r2, s;
double w, a, dd2, b, f, cifval;
double *x, *y;
int *aux;
Geyer *geyer;
geyer = (Geyer *) cdata;
npts = state.npts;
if(npts==0) return ((double) 1.0);
x = state.x;
y = state.y;
u = prop.u;
v = prop.v;
ix = prop.ix;
r2 = geyer->r2;
s = geyer->s;
aux = geyer->aux;
/*
tee = neighbour count at the point in question;
w = sum of changes in (saturated) neighbour counts at other points
*/
tee = w = 0.0;
if(prop.itype == BIRTH) {
if(geyer->per) {
/* periodic distance */
for(j=0; j<npts; j++) {
d2 = dist2(u,v,x[j],y[j],geyer->period);
if(d2 < r2) {
tee++;
f = s - aux[j];
if(f > 1) /* j is not saturated after addition of (u,v) */
w = w + 1; /* addition of (u,v) increases count by 1 */
else if(f > 0) /* j becomes saturated by addition of (u,v) */
w = w + f;
}
}
} else {
/* Euclidean distance */
for(j=0; j<npts; j++) {
d2 = pow(u - x[j], 2) + pow(v - y[j], 2);
if(d2 < r2) {
tee++;
f = s - aux[j];
if(f > 1) /* j is not saturated after addition of (u,v) */
w = w + 1; /* addition of (u,v) increases count by 1 */
else if(f > 0) /* j becomes saturated by addition of (u,v) */
w = w + f;
}
}
}
} else if(prop.itype == DEATH) {
tee = aux[ix];
if(geyer->per) {
/* Periodic distance */
for(j=0; j<npts; j++) {
if(j == ix) continue;
d2 = dist2(u,v,x[j],y[j],geyer->period);
if(d2 < r2) {
f = s - aux[j];
if(f > 0) /* j is not saturated */
w = w + 1; /* deletion of 'ix' decreases count by 1 */
else {
f = f+1;
if(f > 0) {
/* j is not saturated after deletion of 'ix'
(s must be fractional) */
w = w + f;
}
}
}
}
} else {
/* Euclidean distance */
for(j=0; j<npts; j++) {
if(j == ix) continue;
d2 = pow(u - x[j], 2) + pow(v - y[j], 2);
if(d2 < r2) {
f = s - aux[j];
if(f > 0) /* j was not saturated */
w = w + 1; /* deletion of 'ix' decreases count by 1 */
else {
f = f+1;
if(f > 0) {
/* j is not saturated after deletion of 'ix'
(s must be fractional) */
w = w + f;
}
}
}
}
}
} else if(prop.itype == SHIFT) {
/* Compute the cif at the new point, not the ratio of new/old */
for(j=0; j<npts; j++) {
if(j == ix) continue;
d2 = dist2either(u,v,x[j],y[j],geyer->period);
if(d2 < r2) {
tee++;
a = aux[j];
/* Adjust */
dd2 = dist2either(x[ix],y[ix], x[j],y[j],geyer->period);
if(dd2 < r2) a = a - 1;
b = a + 1;
if(a < s && s < b) {
w = w + s - a;
}
else if(s >= b) w = w + 1;
}
}
}
w = w + ((tee < s) ? tee : s);
if(geyer->hard) {
if(tee > 0) cifval = 0.0;
else cifval = 1.0;
}
else cifval = exp(geyer->loggamma*w);
return cifval;
}
void geyerupd(state, prop, cdata)
State state;
Propo prop;
Cdata *cdata;
{
/* Declare other variables */
int ix, npts, j;
double u, v, xix, yix, r2, d2, d2old, d2new;
double *x, *y;
int *aux;
Geyer *geyer;
geyer = (Geyer *) cdata;
aux = geyer->aux;
r2 = geyer->r2;
x = state.x;
y = state.y;
npts = state.npts;
if(prop.itype == BIRTH) {
/* Birth */
u = prop.u;
v = prop.v;
/* initialise auxiliary counter for new point */
aux[npts] = 0;
/* update all auxiliary counters */
if(geyer->per) {
/* periodic distance */
for(j=0; j < npts; j++) {
d2 = dist2(u,v,x[j],y[j],geyer->period);
if(d2 < r2) {
aux[j] += 1;
aux[npts] += 1;
}
}
} else {
/* Euclidean distance */
for(j=0; j < npts; j++) {
d2 = pow(u-x[j], 2) + pow(v-y[j], 2);
if(d2 < r2) {
aux[j] += 1;
aux[npts] += 1;
}
}
}
return;
}
if(prop.itype == DEATH) {
/* Death */
ix = prop.ix;
u = x[ix];
v = y[ix];
/* decrement auxiliary counter for each point */
if(geyer->per) {
/* periodic distance */
for(j=0; j<npts; j++) {
if(j==ix) continue;
d2 = dist2(u,v,x[j],y[j],geyer->period);
if(d2 < r2) {
if(j < ix) aux[j] -= 1;
else aux[j-1] = aux[j] - 1;
} else if(j >= ix) aux[j-1] = aux[j];
}
} else {
/* Euclidean distance */
for(j=0; j<npts; j++) {
if(j==ix) continue;
d2 = pow(u-x[j], 2) + pow(v-y[j], 2);
if(d2 < r2) {
if(j < ix) aux[j] -= 1;
else aux[j-1] = aux[j] - 1;
} else if(j >= ix) aux[j-1] = aux[j];
}
}
return;
}
if(prop.itype == SHIFT) {
/* Shift */
u = prop.u;
v = prop.v;
ix = prop.ix;
xix = x[ix];
yix = y[ix];
/* recompute auxiliary counter for point 'ix' */
aux[ix] = 0;
/* update auxiliary counters for other points */
if(geyer->per) {
for(j=0; j<npts; j++) {
if(j == ix) continue;
d2new = dist2(u,v,x[j],y[j],geyer->period);
d2old = dist2(xix,yix,x[j],y[j],geyer->period);
if(d2old >= r2 && d2new >= r2) continue;
if(d2new < r2) {
/* increment neighbour count for new point */
aux[ix] += 1;
if(d2old >= r2)
aux[j] += 1; /* point j gains a new neighbour */
} else if(d2old < r2)
aux[j] -= 1; /* point j loses a neighbour */
}
} else {
/* Euclidean distance */
for(j=0; j<npts; j++) {
if(j == ix) continue;
d2new = pow(u-x[j], 2) + pow(v-y[j], 2);
d2old = pow(x[ix]-x[j], 2) + pow(y[ix]-y[j],2);
if(d2old >= r2 && d2new >= r2) continue;
if(d2new < r2) {
/* increment neighbour count for new point */
aux[ix] += 1;
if(d2old >= r2)
aux[j] += 1; /* point j gains a new neighbour */
} else if(d2old < r2)
aux[j] -= 1; /* point j loses a neighbour */
}
}
return;
}
fexitc("Unrecognised transition type; bailing out.\n");
}
Cifns GeyerCifns = { &geyerinit, &geyercif, &geyerupd, FALSE};
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...