Revision 1a71709b2bc539256569ffb2c977827979825d64 authored by Adrian Baddeley on 08 April 2012, 05:33:48 UTC, committed by cran-robot on 08 April 2012, 05:33:48 UTC
1 parent 56d752c
PerfectStrauss.h
// ........................... Strauss process ..........................
// $Revision: 1.3 $ $Date: 2012/03/10 11:52:48 $
class StraussProcess : public PointProcess {
public:
double beta, gamma, R, Rsquared;
StraussProcess(double xmin, double xmax, double ymin, double ymax,
double b, double g, double Ri);
~StraussProcess(){}
void NewEvent(double *x, double *y, char *InWindow);
void GeneratePoisson(Point *headPoint,
long int *GeneratedPoints,
long int *LivingPoints,
long int *NoP);
double Interaction(double dsquared);
// void CalcBeta(long int xsidepomm, long int ysidepomm,
// double *betapomm);
// void CheckBeta(long int xsidepomm, long int ysidepomm,
// double *betapomm);
// double lnCondInt(struct Point2 *TempCell, Point2Pattern *p2p);
// void Beta(struct Point2 *TempCell);
// void CalcBeta(Point2Pattern *p2p);
};
StraussProcess::StraussProcess(double xmin, double xmax,
double ymin, double ymax,
double b, double g, double Ri) :
PointProcess(xmin, xmax, ymin, ymax){
beta = b; gamma = g; R = Ri;
Rsquared = R * R;
InteractionRange = R;
TotalBirthRate = beta*(xmax-xmin)*(ymax-ymin);
}
double StraussProcess::Interaction(double dsquared)
{
double rtn;
rtn = 1;
if(dsquared < Rsquared) rtn = gamma;
return(rtn);
}
void StraussProcess::NewEvent(double *x, double *y, char *InWindow)
{
double Xdim, Ydim;
Xdim = Xmax-Xmin;
Ydim = Ymax-Ymin;
*x = slumptal()*Xdim+Xmin;
*y = slumptal()*Ydim+Ymin;
*InWindow = 1;
}
void StraussProcess::GeneratePoisson(Point *headPoint,
long int *GeneratedPoints,
long int *LivingPoints,
long int *NoP)
{
int i;
double xtemp, ytemp, L, Xdim, Ydim;
struct Point *TempPoint;
Xdim = Xmax-Xmin;
Ydim = Ymax-Ymin;
L = beta*Xdim*Ydim;
*GeneratedPoints = poisson(L);
*LivingPoints = *GeneratedPoints;
for (i=1; i<=*GeneratedPoints ; i++){
//Rprintf("Generating StraussProcess Poisson 3\n");
//scanf("%f",&f1);
xtemp = slumptal()*Xdim+Xmin;
ytemp = slumptal()*Ydim+Ymin;
//
//Rprintf("Generating StraussProcess Poisson 3.2\n");
TempPoint = ALLOCATE(struct Point);
//
TempPoint->X = xtemp;
TempPoint->Y = ytemp;
TempPoint->No = i;
TempPoint->R = slumptal();
//Rprintf("Generating StraussProcess Poisson 3.6\n");
TempPoint->next = headPoint->next;
headPoint->next = TempPoint;
*NoP = *NoP + 1;
}
}
//void StraussProcess::CalcBeta(long int xsidepomm, long int ysidepomm,
// double *betapomm){
// long int i,j,k;
// k=0;
// // Rprintf("\ndiagnostic message: Strauss CalcBeta... %ld %ld\n",xsidepomm,ysidepomm);
// for(i=0; i<xsidepomm; i++){
// for(j=0; j<ysidepomm; j++){
// *(betapomm + i*ysidepomm + j) = this->beta;
// k++;
// }
// }
//}
//void StraussProcess::CheckBeta(long int xsidepomm, long int ysidepomm,
// double *betapomm){
// long int i,j,k;
// // double d1;
// k=0;
// // Rprintf("\ndiagnostic message: Strauss CalcBeta... %ld %ld\n",xsidepomm,ysidepomm);
// for(i=0; i<xsidepomm; i++){
// for(j=0; j<ysidepomm; j++){
// if((fabs(*(betapomm + i*ysidepomm + j)- beta)>0.001) && (k==0)){
// Rprintf("%f %f %f %ld %ld\n",fabs(*(betapomm + i*ysidepomm + j)- beta),
// *(betapomm + i*ysidepomm + j),beta,i,j);
// k++;
// // scanf("%lf",&d1);
// }
// }
// }
//}
//double StraussProcess::lnCondInt(struct Point2 *TempCell,
// Point2Pattern *p2p){
// double f1;
// long int xco,yco,xc,yc,fx,tx,fy,ty,ry,rx,k;
// double dy,dx, lnCI,dst2;
// struct Point2 *TempCell2;
//
// f1 = (TempCell->X-p2p->Xmin)/p2p->XCellDim; xc = int(f1);
// CLAMP(xc, 0, p2p->MaxXCell, "xc");
// f1 = (TempCell->Y-p2p->Ymin)/p2p->YCellDim; yc = int(f1);
// CLAMP(yc, 0, p2p->MaxYCell, "yc");
//
// dx = (Xmax-Xmin)/(double(p2p->MaxXCell+1));
// dy = (Ymax-Ymin)/(double(p2p->MaxYCell+1));
// rx = int(this->InteractionRange/dx+1.0);
// ry = int(this->InteractionRange/dy+1.0);
//
// lnCI = log(TempCell->Beta);
//
// k = 0;
//
// if((xc+rx)<=p2p->MaxXCell) tx=xc+rx; else tx = p2p->MaxXCell;
// if((yc+ry)<=p2p->MaxYCell) ty=yc+ry; else ty = p2p->MaxYCell;
// if((xc-rx)>=0) fx=xc-rx; else fx = 0;
// if((yc-ry)>=0) fy=yc-ry; else fy = 0;
//
// //Rprintf("MCI! %d %d %d %d\n",fx,tx,fy,ty);
//
// for(xco = fx; xco <= tx; xco++){
// for(yco = fy; yco <= ty; yco++){
// CHECK(p2p->headCell[xco][yco],
// "internal error: p2p->headCell[xco][yco] is null in lnCondInt()");
// TempCell2 = p2p->headCell[xco][yco]->next;
// CHECK(TempCell2, "internal error: TempCell2 is null in lnCondInt()");
// while(TempCell2!=TempCell2->next){
// if(TempCell2 != TempCell){
// k++;
// dst2 = pow(TempCell->X-TempCell2->X,2)+
// pow(TempCell->Y-TempCell2->Y,2);
// lnCI += log(Interaction(dst2));
// }
// TempCell2 = TempCell2->next;
// CHECK(TempCell2,
// "internal error: TempCell2 is null in lnCondInt() loop");
// }
// }
// }
// return(lnCI);
//}
//void StraussProcess::Beta(struct Point2 *TempCell){
// TempCell->Beta = beta;
//}
//void StraussProcess::CalcBeta(Point2Pattern *p2p){
// long int xco,yco;
// // double dy,dx;
// struct Point2 *TempMother;
//
// for(xco = 0; xco <= p2p->MaxXCell; xco++){
// for(yco = 0; yco <= p2p->MaxYCell; yco++){
// CHECK(p2p->headCell[xco][yco],
// "internal error: p2p->headCell[xco][yco] is null in CalcBeta()");
// TempMother = p2p->headCell[xco][yco]->next;
// CHECK(TempMother, "internal error: TempMother is null in CalcBeta()");
// while(TempMother!=TempMother->next){
// TempMother->Beta = this->beta;
// TempMother = TempMother->next;
// CHECK(TempMother,
// "internal error: TempMother is null in CalcBeta() loop");
// }
// }
// }
//}
// ........................... Interface to R ..........................
extern "C" {
SEXP PerfectStrauss(SEXP beta,
SEXP gamma,
SEXP r,
SEXP xrange,
SEXP yrange) {
// input parameters
double Beta, Gamma, R, Xmin, Xmax, Ymin, Ymax;
double *Xrange, *Yrange;
// internal
int xcells, ycells;
PointProcess *TheProcess;
long int EndTime, StartTime;
// output
int noutmax;
SEXP xout, yout, nout, out;
double *xx, *yy;
int *nn;
SEXP stout, etout;
int *ss, *ee;
// protect arguments from garbage collector
PROTECT(beta = AS_NUMERIC(beta));
PROTECT(gamma = AS_NUMERIC(gamma));
PROTECT(r = AS_NUMERIC(r));
PROTECT(xrange = AS_NUMERIC(xrange));
PROTECT(yrange = AS_NUMERIC(yrange));
// that's 5 protected objects
// extract arguments
Beta = *(NUMERIC_POINTER(beta));
Gamma = *(NUMERIC_POINTER(gamma));
R = *(NUMERIC_POINTER(r));
Xrange = NUMERIC_POINTER(xrange);
Xmin = Xrange[0];
Xmax = Xrange[1];
Yrange = NUMERIC_POINTER(yrange);
Ymin = Yrange[0];
Ymax = Yrange[1];
// compute cell array size
xcells = (int) floor((Xmax-Xmin)/ R);
if(xcells > 9) xcells = 9; if(xcells < 1) xcells = 1;
ycells = (int) floor((Ymax-Ymin)/ R);
if(ycells > 9) ycells = 9; if(ycells < 1) ycells = 1;
#ifdef DBGS
Rprintf("xcells %d ycells %d\n",xcells,ycells);
Rprintf("Initialising\n");
#endif
// Initialise Strauss point process
StraussProcess ExampleProcess(Xmin,Xmax,Ymin,Ymax, Beta, Gamma, R);
// Initialise point pattern
Point2Pattern ExamplePattern(Xmin,Xmax,Ymin,Ymax, xcells, ycells);
// parameters: min x, max x, min y, max y, "cells" in x and y direction
// used for speeding up neighbour counting, 9 is max here
#ifdef DBGS
Rprintf("Initialisation complete\n");
#endif
// Synchronise random number generator
GetRNGstate();
// Initialise perfect sampler
Sampler PerfectSampler(&ExampleProcess);
// Perform perfect sampling
PerfectSampler.Sim(&ExamplePattern, &StartTime, &EndTime);
// Synchronise random number generator
PutRNGstate();
// Get upper estimate of number of points
noutmax = ExamplePattern.UpperCount() + 1;
// Allocate space for output
PROTECT(xout = NEW_NUMERIC(noutmax));
PROTECT(yout = NEW_NUMERIC(noutmax));
PROTECT(nout = NEW_INTEGER(1));
PROTECT(stout = NEW_INTEGER(1));
PROTECT(etout = NEW_INTEGER(1));
xx = NUMERIC_POINTER(xout);
yy = NUMERIC_POINTER(yout);
nn = INTEGER_POINTER(nout);
ss = INTEGER_POINTER(stout);
ee = INTEGER_POINTER(etout);
// copy data into output storage
ExamplePattern.Return(xx, yy, nn, noutmax);
*ss = StartTime;
*ee = EndTime;
// pack up into output list
PROTECT(out = NEW_LIST(5));
SET_VECTOR_ELT(out, 0, xout);
SET_VECTOR_ELT(out, 1, yout);
SET_VECTOR_ELT(out, 2, nout);
SET_VECTOR_ELT(out, 3, stout);
SET_VECTOR_ELT(out, 4, etout);
// return
UNPROTECT(11); // 5 arguments plus xout, yout, nout, stout, etout, out
return(out);
}
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...