https://github.com/cran/RandomFields
Tip revision: 25ed5157dd19e7ca108a1d5ed2581a6a98c376d5 authored by Martin Schlather on 23 February 2011, 00:00:00 UTC
version 2.0.43
version 2.0.43
Tip revision: 25ed515
extremes.cc
/*
Authors
Martin Schlather, martin.schlather@math.uni-goettingen.de
simulation of max-stable random fields
Copyright (C) 2001 -- 2011 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 "RF.h"
#include "Covariance.h"
#include <assert.h>
#include "MPPstandard.h"
void SetExtremes(int *action, double *standardGausMax)
{
extremes_param *gp = &(GLOBAL.extremes);
if (*action) {
if (*standardGausMax<=0.0) {
if (PL>0)
PRINTF("\nERROR! `standardGausMax' not a positive number. Ignored.");
} else gp->standardmax = *standardGausMax;
} else {
*standardGausMax= gp->standardmax;
}
}
void extremes_destruct(void **SExtremes){
if (*SExtremes!=NULL) {
extremes_storage *x;
x = *((extremes_storage**) SExtremes);
if (x->rf != NULL) free(x->rf);
KEY_DELETE(&(x->key));
free(*SExtremes);
*SExtremes = NULL;
// do not free mpp since only pointed
}
}
void doMaxMpp(method_type *Emeth, res_type * res) {
// vassert(false);
/* ***********************************************************
MPP FUNCTIONS
* *********************************************************** */
// what to do with the mean?? -> ignore it
// what to do with the variance?? -> ignore it
GetRNGstate();
extremes_storage* Es = (extremes_storage*) (Emeth->S);
key_type *key = &(Es->key);
method_type *meth = key->meth;
// assert(false);
// jump leading dollars
// note that the variance parameter is and be ignored here, in contrast
// the additive models !
while (meth->nr != -RandomCoin-1 && meth->nsub==1 && meth->sub[0]!=NULL) {
meth = meth->sub[0];
}
assert(meth->nr == -RandomCoin-1);
STANDARDDEFINITIONEN(meth->cov);
// mpp_model randomcoin = C->randomcoin;
// printf("%ld", s);
// PrintModelInfo(cov); // assert(false);
double poisson, invpoisson;
res_type threshold;
long control = 0;
// ntot= s->ntot,
int dim = cov->tsdim,
spatial = loc->totalpoints;
mpp_get mppget = CovList[meth->cov->nr].mppget;
// print("name =%s %d %d %ld %ld %ld\n",
// CovList[meth->cov->nr].name, meth->cov->nr, S2ISO,
// CovList[S2ISO].mppget, CovList[meth->cov->nr].mppget,
// mppget);
// assert(S2ISO == meth->cov->nr);
// assert(mppget != NULL);
STANDARDINIT;
// assert(false);
poiss_initu(s);
s->location = unif_u;
poisson = rexp(1.0);
invpoisson = 1.0 / poisson;
// assert(false);
for(n=0; control < total; n++) {
// n nur zur Unterhaltung; loop allein ueber control !
// printf("%d\n", n);
// printf("%f\n", invpoisson);
randomcoin(s, cov);
Radius = s->effectiveRadius;
if (loc->grid && meth->type <= TypeDiag) {
STANDARDGRID;
if (inside) {
while (true) {
res_type dummy;
// printf("mppget %ld %ld\n", mppget, mppget_Gauss);
dummy = (res_type) mppget(x, s->u, cov, s) * invpoisson;
// assert(false);
// printf("%f u=%f %f %f\n", x[0], s->u[0], dummy, invpoisson);
if (res[zaehler] < dummy) res[zaehler]= dummy;
STANDARDINKREMENT;
}
}
} else { // not a grid
double *x=meth->space;
// the following algorithm can greatly be improved !
// but for ease, just the stupid algorithm
for (zaehler=0; zaehler<spatial; zaehler++, x += dim) {
res_type dummy = (res_type) mppget(x, s->u, cov, s) * invpoisson;
if (res[zaehler] < dummy) res[zaehler] = dummy;
}
}
poisson += rexp(1.0);
invpoisson = 1.0 / poisson;
threshold = (res_type) (s->maxheight * invpoisson);
while ((control<total) && (res[control]>=threshold)) control++;
if (n >= nthreshold) {
nthreshold += deltathresh;
PRINTF("%d %d-th position: value=%f threshold=%f \n",
n, control, (double) res[control], (double) threshold);
}
R_CheckUserInterrupt();
} // while control < total
if ((factor = (res_type) s->factor) != 1.0)
for (zaehler=0; zaehler<total; zaehler++)
res[zaehler] *= factor;
PutRNGstate();
}
void doExtrGauss(method_type *Emeth, res_type * res) {
/* ***********************************************************
EXTREMAL GAUSSIAN
* *********************************************************** */
extremes_storage* Es = (extremes_storage*) (Emeth->S);
key_type *key = &(Es->key);
method_type *meth = key->meth;
location_type *loc = meth->loc;
globalparam *gp = meth->gp;
int i, err,
total = loc->totalpoints,
every = gp->general.every,
nthreshold = (every>0) ? every : MAXINT,
deltathresh = nthreshold;
long n = 0,
control=0;
double poisson, invpoisson;
res_type factor, threshold;
if ((err = internal_DoSimulateRF(&(Es->key), 1, Es->rf)) != NOERROR)
XERR(err);
// to get some consistency with GaussRF concerning Randomseed,
// DoSimulate must be called first; afterwards, it is called after
// creation of Poisson variable
GetRNGstate();
poisson = rexp(1.0);
PutRNGstate();
invpoisson = 1.0 / poisson;
while (true) {
for (i=0; i<total; i++) {// 0, not control, since we approximate
// and 0 as starting control increases precision
res_type dummy = (res_type) (Es->rf[i] * invpoisson);
if (res[i] < dummy) res[i]=dummy;
}
GetRNGstate();
poisson += rexp(1.0);
PutRNGstate();
invpoisson = 1.0 / poisson;
threshold = (res_type) (Es->assumedmax * invpoisson);
while ((control<total) && (res[control]>=threshold)) {control++;}
if (control>=total) break;
if ((err = internal_DoSimulateRF(&(Es->key), 1, Es->rf)) != NOERROR)
XERR(err);
if (n >= nthreshold) {
PRINTF("%d %d-th position: value=%f threshold=%f \n",
n, control, res[control], (double) threshold);
nthreshold += deltathresh;
}
n++;
} // while
if ((factor = (res_type) Es->inv_mean_pos) != 1.0)
for (i=0; i<total; i++) res[i] *= factor;
}
void InitMaxStableRF(double *x, double *T, int *dim, int *lx, int *grid,
int *Time,
int *distr,
int *keyNr,
int expected_number_simu,
int *err)
// for all parameters see InitSimulateRF in RFsimu.cc
{
key_type *key = KEY + *keyNr;
extremes_param *lp;
cov_model *cov = STORED_MODEL[MODEL_SIMU];
extremes_storage* Es;
globalparam newparam;
method_type *Emeth;
strcpy(ERROR_LOC, "MaxStable");
strcpy(PREF_FAILURE, "");
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
*err=ERRORREGISTER;
goto ErrorHandling;
}
// PrintModelInfo(cov);
// assert(false);
KEY_DELETE_NOTTREND(key);
KEY_NULLNOTTREND(key);
if (*Time) {*err=ERRORNOTPROGRAMMED; goto ErrorHandling;}
if (key->trend.TrendModus != TREND_MEAN) {
*err=ERRORTREND; goto ErrorHandling;}
key->simu.distribution=*distr;
*err = loc_set(x, T, *dim, *lx, (bool) *Time, (bool) *grid,
&(key->loc), &GLOBAL);
if (*err) goto ErrorHandling;
memcpy(&(key->gp), &(GLOBAL), sizeof(globalparam));
lp = &(key->gp.extremes);
Emeth = key->meth = (method_type*) malloc(sizeof(method_type));
METHOD_NULL(Emeth);
Emeth->destruct = extremes_destruct; // SET_DESTRUCT(extremes_destruct);
if ((Emeth->S =(extremes_storage*) malloc(sizeof(extremes_storage)))==NULL){
*err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
Es = (extremes_storage*) Emeth->S; // key->meth->S
KEY_NULL(&(Es->key));
Es->rf = NULL;
Emeth->nr = -1-(Es->meth = cov->user[MaxMpp] ? MaxMpp : ExtremalGauss);
// printf("%d %d %d %d\n", MaxMpp, ExtremalGauss, Emeth->nr, cov->user[MaxMpp]); assert(false);
// PrintModelInfo(cov);
// printf("%d %d\n", cov->pref[MaxMpp], cov->pref[ExtremalGauss]);
// assert(false);
// 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
//PrintModelInfo(STORED_MODEL[MODEL_SIMU]);
strcpy(ERROR_LOC, "");
strcpy(PREF_FAILURE, "");
memcpy(&newparam, &(key->gp), sizeof(globalparam));
newparam.general.storing = true;
*err = internal_InitSimulateRF(x, T, *dim, *lx, (bool) *grid,
(bool) *Time,
Es->meth == MaxMpp ? *distr : DISTR_GAUSS,
&(Es->key), &newparam, 1000,
STORED_MODEL + MODEL_SIMU // hier cov
// wuerde zur Katastrophe fuehren
// da STORED_MODEL + MODEL_SIMU von
// internal_InitSimulateRF nicht auf 0 gesetzt
);
cov = Es->key.cov;
covcpy(&(key->cov), cov, false, true);
if (*err!=NOERROR) goto ErrorHandling;
// genuinemeth = Es->key.meth;
// printf("%ld %ld %d\n", genuinemeth->S, genuinemeth, genuinemeth->nr);
// printf("XYOH!!\n");
// PrintMethodInfo(key->meth); assert(false);
// printf("XOH!!\n");
// PrintMethodInfo(genuinemeth); assert(false);
if (Emeth->nr == -1-MaxMpp) {
Emeth->domethod = doMaxMpp;
mpp_storage *mpp;
method_type *genuinemeth = Es->key.meth;
// PrintMethodInfo(genuinemeth);
while (genuinemeth != NULL && genuinemeth->nr != -RandomCoin-1) {
genuinemeth = genuinemeth->sub[0];
}
if (genuinemeth == NULL) {
*err=ERRORNOGENUINEMETHOD;
goto ErrorHandling;
}
// while (genuinemeth->nr != -RandomCoin-1 && genuinemeth->nsub==1 && genuinemeth->sub[0]!=NULL) {
// if (genuinemeth->nr != DOLLAR) {
// PrintMethodInfo(meth);
// ERR("unexpected method for shape function -- please contact author");
// }
// // printf("%d\n", genuinemeth->nr);
// genuinemeth = genuinemeth->sub[0];
// }
// printf("%d %d %d %d\n", genuinemeth->nr, -RandomCoin, GATTER, LASTGATTER);
// if (genuinemeth->nr != -RandomCoin-1)
// ERR("shape function could not be detected -- please contact author");
mpp = (mpp_storage*) genuinemeth->S;
assert(mpp != NULL);
// printf("XOH!! %d\n", mpp);
// PrintMethodInfo(Es->key.meth);
// MaxMpp may not be mixed with other methods !
/* ***********************************************************
MPP FUNCTIONS
* *********************************************************** */
// for (m=0; m<=n; m++) key->meth[m].unimeth = MaxMpp;
// the basic functions in MPPFcts.cc are kept quite simple (e.g. height 1.0)
// hence height must be corrected to "unit volume"; effectiveara is
// involved since below only standard exponential variables are considered;
// not that effectivearea is not always the simulation window !
// xxx corrects for the max-superposition of the xxx maxstable rf
// printf("%d\n", genuinemeth->S);
// printf("%f\n", mpp->effectivearea);
// printf("%f\n", mpp->integralpos);
mpp->factor = mpp->effectivearea / (mpp->integralpos);
// printf("factor = %f %f %f\n",
// mpp->factor, mpp->effectivearea, (mpp->integralpos));
// assert(false);
// assert(false); // kkkk
// s = (mpp_storage*) key->meth[n].S;
// s->factor = s->effectivearea / (s->integralpos * (double) (key->m));
// for (v=0; v<n; v++) {
// mpp_storage* sp1;
// s = (mpp_storage*) key->meth[v].S;
// sp1 = (mpp_storage*) key->meth[v+1].S;
// s->factor = s->effectivearea / s->integralpos /
// (sp1->effectivearea / sp1->integralpos);
// }
} else {
/* ***********************************************************
EXTREMAL GAUSSIAN
* *********************************************************** */
// PrintModelInfo(cov);
// printf("%d\n", cov->user[ExtremalGauss]);
// assert(cov->user[ExtremalGauss]);
assert(cov->user[MaxMpp] == PREF_NONE);
assert(*distr == DISTR_MAXSTABLE);
Emeth->domethod = doExtrGauss;
double sigma, meanDsigma, mean;
key->simu.distribution = *distr;
CovList[cov->nr].cov(ZERO, cov, &sigma);
if (sigma==0) {*err = ERRORSILLNULL; goto ErrorHandling;}
sigma = sqrt(sigma);
mean = Es->key.trend.mean;
meanDsigma = mean / sigma;
Es->inv_mean_pos =
1.0 / (sigma * INVSQRTTWOPI * exp(-0.5 * meanDsigma * meanDsigma)
+ mean * /* correct next lines the 2nd factor is given */
// pnorm(meanDsigma, key->mean, sigma,1,0));
pnorm(0.0, mean, sigma, 1, 0)
);
Es->assumedmax = lp->standardmax * sigma + ((mean>0) ? mean : 0);
*err=NOERROR;
key->simu.active = true;
if ((Es->rf =
(res_type*) malloc(sizeof(res_type) * Es->key.loc.totalpoints))==NULL) {
*err=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
}
strcpy(ERROR_LOC, "");
key->simu.active = true;
return;
ErrorHandling:
char ErrM[500], ErrM2[800];
key->simu.active=false;
// printf("!!!!!!! %s", PREF_FAILURE);
errorMSG(*err, ErrM);
sprintf(ErrM2, "%s%s", PREF_FAILURE, ErrM);
ERR(ErrM2);
// return;
}
void DoMaxStableRF(int *keyNr, int *n, int *pairs, res_type *res, int *err)
{
key_type *key;
long total;
res_type *RES;
int ni, i;
method_type *meth;
*err=0;
key = NULL;
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
*err=ERRORKEYNROUTOFRANGE; goto ErrorHandling;
}
if (*pairs) {*err=ERRORPAIRS; goto ErrorHandling;}
key = &(KEY[*keyNr]);
if (!key->simu.active) {*err=ERRORNOTINITIALIZED; goto ErrorHandling;}
total = ((extremes_storage*)key->meth->S)->key.loc.totalpoints;
meth = key->meth;
for (RES = res, ni=0; ni<*n; ni++, RES += total) {
for (i=0; i<total; i++) RES[i]=0.0;
meth->domethod(meth, RES);
}
return;
ErrorHandling:
if (PL>0) ErrorMessage(Nothing,*err);
if (key!=NULL) key->simu.active=false;
}