/* Authors Martin Schlather, martin.schlather@cu.lu Copyright (C) 2001 -- 2004 Martin Schlather, This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /* TO DO * noch fast alles auf key->x programmiert statt s->x * time component: aenderungen siehe TBM ! */ #include #include #include "RFsimu.h" #include "MPPFcts.h" Real MPP_RADIUS =0.0; Real MPP_APPROXZERO =0.001; Real ADDMPP_REALISATIONS = 100.0; // logically an integer, implementation // allows for positive real ! void SetMPP(int *action, Real *approxzero, Real *realisations, Real *radius) { switch(*action) { case 0 : if (*approxzero<=0.0) { if (GENERAL_PRINTLEVEL>0) PRINTF("\nERROR! `approx.zero' not a positive number. Ignored.\n"); } else MPP_APPROXZERO = *approxzero; if (*realisations<=0.0) { if (GENERAL_PRINTLEVEL>0) PRINTF("\nERROR! `realisation' not a positive number. Ignored.\n"); } else ADDMPP_REALISATIONS = *realisations; MPP_RADIUS = *radius; if ((MPP_RADIUS>0.0) && (GENERAL_PRINTLEVEL>1)) PRINTF("\nWARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect.\n"); break; case 1 : *approxzero = MPP_APPROXZERO; *realisations = ADDMPP_REALISATIONS; *radius = MPP_RADIUS; if (GetNotPrint) break; case 2 : PRINTF("\nmarked point process methods\n============================\napproxzero=%e,\nrealisations=%e\n", MPP_APPROXZERO,ADDMPP_REALISATIONS); break; default : if (GENERAL_PRINTLEVEL>0) PRINTF("unknown action\n"); } } void mpp_destruct(void **S) { // mpp_storage in RFsimu.h !! if (*S!=NULL) { free(*S); *S = NULL; } } void mpp_destructX(void **SX) { // mpp_storage in RFsimu.h !! if (*SX!=NULL) { free(*SX); *SX = NULL; } } int init_mpp(key_type * key, int m) { int error, d, i, v, start_param[MAXDIM], index_dim[MAXDIM]; bool time_exception[MAXCOV], no_last_comp /* , any_time_exception*/; Real max[MAXDIM]; mpp_storage *s; cov_fct *cov; char actcov; int covnr[MAXCOV]; if (key->anisotropy) { error=ERRORNOTPROGRAMMED; goto ErrorHandling; // RFtest.new.features.1.01.R: MaxStableRF is excluded by if (FALSE) // excludes also any time component // // even in the case of isotropy the method is not checked yet! } SET_DESTRUCT(mpp_destruct); assert(key->destructX==NULL); assert(key->S[m]==NULL); if ((key->S[m]=malloc(sizeof(mpp_storage)))==0){ error=ERRORMEMORYALLOCATION; goto ErrorHandling; } s = (mpp_storage*) key->S[m]; actcov=0; { int v; for (v=0; vncov; v++) { if ((key->method[v]==AdditiveMpp) && (key->left[v]) && (key->param[v][VARIANCE]>0)) { key->left[v]=false; assert(key->covnr[v]>=0 && key->covnr[v]covnr[v]; if (CovList[covnr[actcov]].add_mpp_scl==NULL) {error=ERRORNOTDEFINED; goto ErrorHandling;} memcpy(s->param[actcov], key->param[v], sizeof(Real) * key->totalparam); if ((vncov-1 && key->op[v]) || (v>0 && key->op[v-1])) { // no multiplicative models are allowed error=ERRORNOMULTIPLICATION; goto ErrorHandling; } actcov++; break; } } } if (actcov==0) { if (key->traditional) error=ERRORNOTINITIALIZED; else error=NOERROR_ENDOFLIST; goto ErrorHandling; } s->actcov = actcov; // transform points using anisotropy and get s->timespacedim //if ((param[0][key->totalparam-1]==0.0) && any_time_exception) { // param[0][key->totalparam-1]=1.0; //} // see tbm s->grid = key->grid && !key->anisotropy; GetTrueDim(key->anisotropy, key->timespacedim, s->param[0], &s->timespacedim, &no_last_comp, start_param, index_dim); if (error=Transform2NoGrid(key, s->param[0], s->timespacedim, start_param, &(s->x))) goto ErrorHandling; // are methods and parameters fine ? for (v=0; vTime) && no_last_comp && (cov->isotropic==SPACEISOTROPIC)) { error= ERRORWITHOUTTIME;goto ErrorHandling;} else if ((cov->check!=NULL) && (error=cov->check(s->param[v], s->timespacedim, AdditiveMpp))) goto ErrorHandling; } // determine the minimum area where random field values are to be generated if (s->grid) { for (d=0; dtimespacedim; d++) { s->min[d] = key->x[d][XSTART]; s->length[d]= key->x[d][XSTEP] * (Real) (key->length[d] - 1); // x[XSTART] and x[XSTEP] are assumed to be precise, but // not x[XEND] } } else { int ix; for (d=0;dmin[d]=RF_INF; max[d]=RF_NEGINF;} for (ix=i=0;itotalpoints;i++,ix+=s->timespacedim) { for (d=0; dtimespacedim; d++) { if (s->x[ix+d] < s->min[d]) s->min[d] = s->x[ix+d]; if (s->x[ix+d] > max[d]) max[d] = s->x[ix+d]; } } for (d=0;dlength[d] = max[d] - s->min[d]; } for (v=0; vparam[v][SCALE]=s->param[v][INVSCALE]=1.0; else assert(false); // no anisotropic function programmed yet; // neither considered what the extension should // look like } } for (i=0; iaddradius[i] = MPP_RADIUS; // must be set BEFORE the next command!! CovList[covnr[i]].add_mpp_scl(s, i); s->MppFct[i] = CovList[covnr[i]].add_mpp_rnd; } if ((MPP_RADIUS>0.0) && (GENERAL_PRINTLEVEL>=2)) PRINTF("Note: area has been enlarged by fixed value (%f)", s->addradius[0]); if (key->anisotropy) return NOERROR_REPEAT; else return 0; ErrorHandling: return error; } void do_addmpp(key_type *key, int m, Real *res ) { int i, d, dim, v; Real lambda[MAXDIM], min[MAXDIM], max[MAXDIM], factor, average, poisson; long segment[MAXDIM+1]; mpp_storage *s; mppmodel model; assert(key->active); { register long i,endfor; endfor = key->totalpoints; for (i=0;iS[m]; switch (key->distribution) { case DISTR_GAUSS : assert(s->actcov==1); lambda[0] = ADDMPP_REALISATIONS * s->effectivearea[0]; factor = sqrt(s->param[0][VARIANCE] / (ADDMPP_REALISATIONS *s->integralsq[0])); average = ADDMPP_REALISATIONS * s->integral[0]; break; case DISTR_POISSON : // not done yet -- it is only a fragment // on the R level, it is excluded that the program runs // into this part // // actually, only actcov==1 is possible at the moment! // // maybe nice to switch MPPFcts back to actcov=1 assert(s->actcov==1); for (v=0; vactcov; v++) { lambda[v] = s->param[v][VARIANCE] * s->effectivearea[v]; factor = average = 0.0 /* replace by true value !! */; assert(false); if(key->mean!=0) { // do not compare mean against VARIANCE since // VARIANCEs may differ with v // if (GENERAL_PRINTLEVEL>0) ErrorMessage(AdditiveMpp,ERRORVARMEAN); // print it always assert(false); // for debugging -- delete this later on // finally, the value of MEAN is ignored, and that of VARIANCE // is used for simulating the random field } } break; default : assert(false); } segment[0] = 1; // only used for s->grid! for (d=0; dtimespacedim; d++) { assert(dlength[d]; } for (i=0; itotalpoints; i++) res[i]=0.0; for (v=0; vactcov; v++) { for(poisson = -log(1.0 - UNIFORM_RANDOM); poisson < lambda[v]; poisson -= log(1.0 - UNIFORM_RANDOM)) { (s->MppFct[v])(s, v, min, max, &model ); if (s->grid) { int start[MAXDIM], end[MAXDIM], resindex, index[MAXDIM], segmentdelta[MAXDIM], dimM1; Real coord[MAXDIM], startcoord[MAXDIM]; // determine rectangle of indices, where the mpp function // is different from zero for (d=0; dtimespacedim; d++) { if (min[d]< key->x[d][XSTART]) {start[d]=0;} else start[d] = (int) ( (min[d] - key->x[d][XSTART]) / key->x[d][XSTEP]); // "end[d] = 1 + *" since inequalities are "* < end[d]" // not "* <= end[d]" ! end[d] = 1 + (int) ((max[d] - key->x[d][XSTART]) / key->x[d][XSTEP]); if (end[d] > key->length[d]) end[d] = key->length[d]; } // prepare coord starting value and the segment vectors for res resindex = 0; for (d=0; dtimespacedim; d++) { index[d]=start[d]; segmentdelta[d] = segment[d] * (end[d]-start[d]); resindex += segment[d] * start[d]; coord[d] = startcoord[d] = key->x[d][XSTART] + (Real)start[d] * key->x[d][XSTEP]; } // add mpp function to res dimM1 = s->timespacedim -1; d = 0; // only needed for assert(dtimespacedim); res[resindex] += model(coord); d=0; index[d]++; coord[d]+=key->x[d][XSTEP]; resindex++; while (index[d] >= end[d]) { if (d>=dimM1) break; // below not possible // as dim==1 will fail! index[d]=start[d]; coord[d]=startcoord[d]; resindex -= segmentdelta[d]; d++; // if (d>=dim) break; index[d]++; coord[d]+=key->x[d][XSTEP]; resindex += segment[d]; } } } else { // !s->grid Real y[MAXDIM]; long j; // the following algorithm can greatly be improved ! // but for ease, just the stupid algorithm for (j=i=0; itotalpoints; i++) { for (d=0; dtimespacedim; d++,j++) { y[d] = s->x[j]; } res[i] += model(y); } } } } if (key->distribution==DISTR_GAUSS) { for (i=0; itotalpoints;i++) { res[i] = factor * (res[i] - average); } } }