/*
Authors
Martin Schlather, martin.schlather@cu.lu
Collection of auxiliary functions
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 PURSE. 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.
*/
#include <math.h>
#include <stdio.h>
#include "RFsimu.h"
#include "MPPFcts.h"
#include <assert.h>
Real EXTREMES_STANDARDMAX = 3.0;
void SetExtremes(int *action, Real *standardGausMax)
{
switch(*action) {
case 0 :
if (*standardGausMax<=0.0) {
if (GENERAL_PRINTLEVEL>0) PRINTF("\nERROR! `standardGausMax' not a positive number. Ignored.");
} else EXTREMES_STANDARDMAX = *standardGausMax;
break;
case 1 :
*standardGausMax= EXTREMES_STANDARDMAX;
if (GetNotPrint) break;
case 2 :
PRINTF("\nMax-stable Random fields\n========================\nstandardGausMax=%e\n",
EXTREMES_STANDARDMAX);
break;
default : if (GENERAL_PRINTLEVEL>0) PRINTF("unknown action\n");
}
}
typedef struct extremes_storage{
Real *rf;
Real inv_mean_pos, assumedmax;
} extremes_storage;
void extremes_destruct(void **SX){
if (*SX!=NULL) {
extremes_storage *x;
x = *((extremes_storage**) SX);
if (x->rf != NULL) free(x->rf);
free(*SX);
*SX = NULL;
}
}
void InitMaxStableRF(Real *x, Real *T, int *dim, int *lx, int *grid,
int *Time,
int *covnr, Real *ParamList, int *nParam,
Real *mean,
int *ncov, int *anisotropy, int *op,
int *method,
int *distr,
int *keyNr,
int *error)
/* for all parameters except precision see InitSimulateRF in the
* package rf!
*
* assumedmax: assumed maximum value of a Gaussian random field;
* should be about 3 or 4
*
*/
{
Real sigma,meanDsigma;
extremes_storage *es;
int init_method[MAXCOV] ;
int gauss_distr = DISTR_GAUSS;
int i;
key_type *key;
bool storing;
storing = GENERAL_STORING;
GENERAL_STORING = true;
key = NULL;
es = NULL;
// I cannot see why a time component might be useful
// hence cathegorically do not allow for the use!
// if you allow for it check if it use is consistent with the
// rest
if (*Time) {*error=ERRORNOTPROGRAMMED; goto ErrorHandling;}
if (method[0] == (int) MaxMpp) {
/* ***********************************************************
MPP FUNCTIONS
* *********************************************************** */
if (*ncov!=1) {*error=ERRORFAILED; goto InternalErrorHandling;}
for (i=0; i<*covnr; i++) init_method[i] = (int) AdditiveMpp;
InitSimulateRF(x, T, dim, lx, grid, Time, covnr, ParamList, nParam,
mean, ncov, anisotropy, op,
init_method, distr, keyNr, error);
key = &(KEY[*keyNr]);
assert(*error>=0);
if (*error!=NOERROR) {goto ErrorHandling;}
key->method[0] = MaxMpp;
assert( (key->SX==NULL) && (key->destructX==NULL));
} else {
/* ***********************************************************
EXTREMAL GAUSSIAN
* *********************************************************** */
InitSimulateRF(x, T, dim, lx, grid, Time, covnr, ParamList, nParam,
mean, ncov, anisotropy, op,
method, &gauss_distr, keyNr, error);
key = &(KEY[*keyNr]);
assert(*error>=0);
if (*error!=NOERROR) {goto ErrorHandling;}
key->distribution = DISTR_MAXSTABLE;
if (key->SX==NULL) {
key->destructX = extremes_destruct;
if ((key->SX=
(extremes_storage*) malloc(sizeof(extremes_storage)))==NULL) {
*error=ERRORMEMORYALLOCATION; goto InternalErrorHandling;
}
es = (extremes_storage*) key->SX;
if ((es->rf =
(Real*) malloc(sizeof(Real) * key->totalpoints))==NULL) {
*error=ERRORMEMORYALLOCATION; goto InternalErrorHandling;
}
} else {es = (extremes_storage*) key->SX;}
sigma = sqrt(CovFct(ZERO,*dim,covnr,op,key->param,*ncov,*anisotropy));
if (sigma==0) {*error = ERRORSILLNULL; goto InternalErrorHandling;}
meanDsigma = key->mean / sigma;
es->inv_mean_pos =
1.0 / (sigma * INVSQRTTWOPI * exp(-0.5 * meanDsigma * meanDsigma)
+ key->mean * /* correct next lines the 2nd factor is given */
// pnorm(meanDsigma, key->mean, sigma,1,0));
pnorm(0.0, key->mean, sigma,1,0)
);
for (i=0;i<key->totalpoints;i++) es->rf[i]=0.0;
es->assumedmax = EXTREMES_STANDARDMAX * sigma +
((key->mean>0) ? key->mean : 0);
*error=NOERROR;
}
GENERAL_STORING = key->storing = storing;
return;
InternalErrorHandling:
ErrorMessage(Nothing,*error);
ErrorHandling:
GENERAL_STORING = storing;
key->active=false;
return;
}
void DoMaxStableRF(int *keyNr, Real *res, int *error)
{
key_type *key;
extremes_storage *es;
long totalpoints,i,control;
Real poisson, invpoisson, threshold, *zw, *RES;
long counter=0; // just for printing messages
bool next=false, storing;
*error=0;
zw = NULL;
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
*error=ERRORKEYNROUTOFRANGE; goto ErrorHandling;
}
key = &(KEY[*keyNr]);
if (!key->active) {*error=ERRORNOTINITIALIZED; goto ErrorHandling;}
storing = key->storing;
key->storing = true;
if (key->method[0] == (int) MaxMpp) {
/* ***********************************************************
MPP FUNCTIONS
* *********************************************************** */
// what to do with the mean?? -> ignore it
// what to do with the variance?? -> ignore it
int dim, d, m, v;
Real min[MAXDIM], max[MAXDIM], factor;
long segment[MAXDIM+1];
mpp_storage *s;
mppmodel model;
MPPRandom MppFct;
GetRNGstate();
assert(key->ncov>0);
if (key->ncov>1) {
// not used
if ((zw=(Real *)malloc(sizeof(Real)*key->totalpoints))==0){
*error=ERRORMEMORYALLOCATION;goto ErrorHandling;}
RES = zw;
for (i=0; i<key->totalpoints; i++) res[i]=zw[i]=0.0;
} else {
for (i=0; i<key->totalpoints; i++) res[i]=0.0;
RES = res;
}
s = (mpp_storage*) key->S[0];
dim = s->timespacedim;
assert(dim>0);
segment[0] = 1;
for (d=0; d<dim; d++)
segment[d+1] = segment[d] * key->length[d];
totalpoints = key->totalpoints;
for (v=0; v<s->actcov; v++) {
MppFct = s->MppFct[v];
factor = s->effectivearea[v] / s->integralpos[v];
control = 0;
poisson = -log(1.0 - UNIFORM_RANDOM); /* it is important to write 1-U
as the random generator returns
0 inclusively; the latter would
lead to an error !
*/
invpoisson = 1.0 / poisson;
while (control<totalpoints) {
// the following is essentially the same as the code in MPP.cc!!
// how to reorganise it to get rid of the duplication of the code??
//
MppFct(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; d<dim; d++) {
if (next = ((min[d] > key->x[d][XEND]) ||
(max[d] < key->x[d][XSTART]))) { break;}
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] = (int) ((max[d] - key->x[d][XSTART]) /
key->x[d][XSTEP]);
if (end[d] > key->length[d]) end[d] = key->length[d];
}
if (!next) {
// prepare coord starting value and the segment vectors for res
resindex = 0;
for (d=0; d<dim; 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 = dim - 1;
d = 0; // only needed for assert(d<dim)
while (index[dimM1]<end[dimM1]) {
register Real dummy;
assert(d<dim);
dummy = model(coord) * invpoisson;
if (RES[resindex] < dummy) RES[resindex]=dummy;
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];
}
}
} /* next */
} else { // not agrid
Real y[MAXDIM];
long j;
// the following algorithm can greatly be improved !
// but for ease, just the stupid algorithm
for (j=i=0; i<totalpoints; i++) {
for (d=0; d<dim; d++,j++) { y[d] = s->x[j]; }
register Real dummy;
dummy = model(y) * invpoisson;
if (RES[i] < dummy) RES[i]=dummy;
}
}
poisson -= log(1.0 - UNIFORM_RANDOM);
invpoisson = 1.0 / poisson;
threshold = s->maxheight[v] * invpoisson;
while ((control<totalpoints) && (RES[control]>=threshold)) {control++;}
if (GENERAL_PRINTLEVEL>=3) {
counter++;
if ((counter % 100==0) && (control<totalpoints))
PRINTF("%d: %d %f %f \n",counter,control,RES[control],
threshold);
}
} // while control < totalpoints
if (key->ncov>1) { // void -- maybe for future extensions
assert(false);
register Real dummy;
for (i=0; i<totalpoints;i++) {
if (RES[i] < (dummy=RES[i]*factor)) RES[i]=dummy;
RES[i] = 0.0;
}
} else {
for (i=0; i<totalpoints;i++) RES[i] *= factor;
}
} //v
//if (key->param[NUGGET]>0) {
// Real nugget;
// nugget = key->param[NUGGET];
// for (i=0; i<totalpoints; i++) {
// register Real dummy;
// dummy = - nugget / log(1.0 - UNIFORM_RANDOM);
// if (res[i] < dummy) res[i]=dummy;
// }
//}
PutRNGstate();
} else {
/* ***********************************************************
EXTREMAL GAUSSIAN
* *********************************************************** */
if (key->distribution!=DISTR_MAXSTABLE) {
*error=ERRORWRONGINIT;
goto InternalErrorHandling;
}
es = (extremes_storage*) key->SX;
assert(es!=NULL);
totalpoints = key->totalpoints;
control = 0;
DoSimulateRF(keyNr, es->rf, error);
// to get some consistency with GaussRF concerning Randomseed,
// DoSimulate must be called first; afterwards, it is called after
// creation of Poisson variable
if (*error) {goto ErrorHandling;}
if (!key->active){ *error=ERRORNOTINITIALIZED; goto InternalErrorHandling; }
GetRNGstate();
poisson = -log(1.0 - UNIFORM_RANDOM);
// it is important to write 1-U as the random generator returns
// 0 inclusively; the latter would lead to an error !
PutRNGstate();
invpoisson = 1.0 / poisson;
while (true) {
for (i=0; i<totalpoints; i++) {// 0, not control, since we approximate
// and 0 as starting control increases precision
es->rf[i] *= invpoisson;
if (res[i] < es->rf[i]) res[i]=es->rf[i];
}
GetRNGstate();
poisson -= log(1.0 - UNIFORM_RANDOM); // !! -= !!
PutRNGstate();
invpoisson = 1.0 / poisson;
threshold = es->assumedmax * invpoisson;
while ((control<totalpoints) && (res[control]>=threshold)) {control++;}
if (control>=totalpoints) break;
DoSimulateRF(keyNr, es->rf, error);
// first time no error; no reason to assume that an error occurs now;
// so error is not checked
if (GENERAL_PRINTLEVEL>=3) {
counter++;
if (counter % 10==0)
PRINTF("%d: %d %e %e \n",counter,control,invpoisson,threshold);
}
} // while
for (i=0; i<key->totalpoints; i++)
res[i] *= es->inv_mean_pos;
}
key->storing = storing;
key->active=GENERAL_STORING && key->storing;
if (zw!=NULL) free(zw);
return;
InternalErrorHandling:
if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,*error);
ErrorHandling:
if (zw!=NULL) free(zw);
key->active=false;
}
