https://github.com/cran/RandomFields
Tip revision: 6332d8d86088cebf6f828f1d29c71b8060f7e88b authored by Martin Schlather on 23 June 2016, 09:04:50 UTC
version 3.1.16
version 3.1.16
Tip revision: 6332d8d
extremes.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
simulation of max-stable random fields
Copyright (C) 2001 -- 2015 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 3
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 "Operator.h"
#define POISSON_INTENSITY 0
#define RANDOMCOIN_INTENSITY (COMMON_GAUSS + 1)
#define RANDOMCOIN_METHOD (COMMON_GAUSS + 2)
/* not used
double GetTotal(cov_model* cov) {
cov_fct *C = CovList + cov->nr;
double v = 1.0;
int i,
nsub = cov->nsub,
kappas = C->kappas;
if (C->Type == RandomType) {
if (cov->total_n >= 0) {
int VARIABLE_IS_NOT_USED old_n = cov->total_n;
if (cov->total_n < GLOBAL.distr.repetitions) {
double *r = (double*) MALLOC(cov->xdimown * sizeof(double));
while(cov->total_n < GLOBAL.distr.repetitions) {
DORANDOM(cov, r);
assert(cov->total_n > cov->total_n);
}
FREE(r);
}
v *= cov->total_sum / (double) cov->total_n;
}
} else {
for (i=0; i<nsub; i++) {
v *= GetTotal(cov->sub[i]);
}
for (i=0; i<kappas; i++) {
if (cov->kappasub[i] != NULL)
v *= GetTotal(cov->kappasub[i]);
}
}
return v;
}
*/
int addStandard(cov_model **Cov) {
cov_model *cov = (*Cov)->calling,
*next = *Cov;
int err,
dim = next->xdimprev,
vdim = next->vdim[0],
role = next->role;
assert(cov->vdim[0] == cov->vdim[1]);
addModel(Cov, STANDARD_SHAPE);
cov_model *key = *Cov;
SetLoc2NewLoc(key, PLoc(cov));
assert(key->calling == cov);
assert(key->sub[PGS_LOC] == NULL && key->sub[PGS_FCT] != NULL);
#define checkstandardnext\
if ((err = CHECK(key, dim, dim, PointShapeType, XONLY, \
CoordinateSystemOf(cov->isoown), \
vdim, role)) != NOERROR) return err
checkstandardnext;
if (!CallingSet(*Cov)) BUG;
if (hasPoissonRole(cov)) {
addModel(key, PGS_LOC, UNIF);
assert(key->nsub == 2);
} else {
if ((err = STRUCT(key, key->sub + PGS_LOC)) != NOERROR) return err;
key->sub[PGS_LOC]->calling = key;
}
if (!CallingSet(*Cov)) BUG;
checkstandardnext;
return NOERROR;
}
void range_mpp(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
// assert(cov != NULL);
range->min[GEV_XI] = RF_NEGINF;
range->max[GEV_XI] = RF_INF;
range->pmin[GEV_XI] = -10;
range->pmax[GEV_XI] = 10;
range->openmin[GEV_XI] = true;
range->openmax[GEV_XI] = true;
range->min[GEV_MU] = RF_NEGINF;
range->max[GEV_MU] = RF_INF;
range->pmin[GEV_MU] = -1000;
range->pmax[GEV_MU] = 1000;
range->openmin[GEV_MU] = true;
range->openmax[GEV_MU] = true;
range->min[GEV_S] = 0;
range->max[GEV_S] = RF_INF;
range->pmin[GEV_S] = 0.001;
range->pmax[GEV_S] = 1000;
range->openmin[GEV_S] = true;
range->openmax[GEV_S] = true;
}
int init_mpp(cov_model *cov, gen_storage *S) {
cov_model *sub = cov->key != NULL ? cov->key :
cov->sub[0] != NULL ? cov->sub[0] : cov->sub[1];
if (sub == NULL) SERR("substructure could be detected");
location_type *loc = Loc(cov);
//window_info *w = &(S->window);
int err,
// dim = cov->tsdim,
role = cov->role,
maxstable = hasMaxStableRole(cov);
pgs_storage *pgs;
if (!maxstable && !hasPoissonRole(cov)) ILLEGAL_ROLE;
if (!isPointShape(sub)) SERR1("%s is not shape/point function", NICK(sub));
if (loc->distances) return ERRORFAILED;
if ((err = INIT(sub, maxstable ? 1 : role == ROLE_POISSON ? 0 : 2, S))
!= NOERROR) return err;
pgs = sub->Spgs;
assert(pgs != NULL);
GetDiameter(loc, pgs->supportmin, pgs->supportmax, pgs->supportcentre);
if (maxstable) {
if (!R_FINITE(sub->mpp.mMplus[1]) || sub->mpp.mMplus[1] <= 0.0) {
// APMI(sub);
// printf("mplus %s %f\n", NICK(sub), sub->mpp.mMplus[1]);
// crash();
SERR("integral of positive part of submodel unkown");
}
//PMI(sub, "sub");
if (!R_FINITE(pgs->zhou_c) && !R_FINITE(pgs->sum_zhou_c))
SERR("maximal height of submodel unkown or the set of locations exceeds possibilities of the programme");
} else if (role == ROLE_POISSON_GAUSS) {
if (ISNAN(sub->mpp.mM[2]) || !R_FINITE(sub->mpp.mM[2] || sub->mpp.mM[2] <=0.0)){
SERR("second moment unkown");
}
} // else role == ROLE_POISSON -- nothing to do
if ((err = FieldReturn(cov)) != NOERROR) { // must be later than INIT !
return err;
}
cov->simu.active = err == NOERROR;
return err;
}
#define SET_SUB { \
sub = (cov->key != NULL) ? cov->key : /* if changed, change do_max_pgs! */ \
cov->sub[0] != NULL ? cov->sub[0] : cov->sub[1]; \
if (sub == NULL) ERR("substructure could be detected"); \
pgs = sub->Spgs; \
assert(pgs != NULL); \
gridlen = pgs->gridlen; \
end = pgs->end; \
start = pgs->start; \
delta = pgs->delta; \
nx = pgs->nx; \
xstart = pgs->xstart; /* nach DO gesetzt */ \
x = pgs->x; /* nach DO gesetzt */ \
inc = pgs->inc;}
void dompp(cov_model *cov, gen_storage *s, double *simuxi) {
// printf("here\n");
assert(cov->simu.active);
location_type *loc = Loc(cov);
if (loc->caniso != NULL) BUG;
cov_model *sub = NULL;
pgs_storage *pgs = NULL;
long spatial = loc->totalpoints;
int d, err,
*gridlen = NULL,
*end = NULL,
*start = NULL,
*delta = NULL,
*nx = NULL,
dim = cov->tsdim,
every = GLOBAL.general.every, //
nthreshold = (every>0) ? every : MAXINT,
// covrole = cov->role,
subrole = ROLE_FAILED,
deltathresh = nthreshold;
double logdens, factor, summand, Sign[2], value[2], logthreshold,
*xstart = NULL, // nach DO gesetzt
*x = NULL, // nach DO gesetzt
*inc = NULL,
*rf = NULL,
//M1 = RF_NA,
logM2 = RF_NA,
gumbel = RF_NA,
frechet = RF_NA,
threshold = RF_NA,
poisson=0.0;
double *res = cov->rf;
assert(res != NULL);
ext_bool fieldreturn;
bool
compoundpoisson, poissongauss,
maxstable = hasMaxStableRole(cov),
simugrid = loc->grid,
partialfield = false,
// AVERAGE braucht Sonderbehandlung:
randomcoin = cov->method == RandomCoin //only for ROLE_ POISSON_ GAUSS
;
long k, i, zaehler, n, cumgridlen[MAXMPPDIM +1], Total_n,
total_pts = loc->totalpoints,
vdim = cov->vdim[0],
vdimtot = total_pts * vdim,
control = 0;
assert(cov->vdim[0] == cov->vdim[1]);
cov_model *key = cov->key;
ext_bool loggiven;
double Minimum = RF_NEGINF; // -10 // TO DO
// printf("A\n");if (PL == 107) BUG;
if (vdim != 1)
ERR("Poisson point process based methods only work in the univariate case");
SET_SUB;
//printf("B\n");if (PL == 108) BUG;
if (!isPointShape(sub)) BUG;
subrole = sub->role;
compoundpoisson = subrole == ROLE_POISSON;
poissongauss = subrole == ROLE_POISSON_GAUSS;
Total_n = maxstable ? -1 : rpois(pgs->intensity);
assert(maxstable || Total_n > 0);
loggiven = key == NULL ? cov->sub[0]->loggiven : key->loggiven;
if (!loggiven) Minimum = exp(Minimum);
if (simugrid) {
cumgridlen[0] = 1;
for (d=0; d<dim; d++) {
inc[d] = loc->xgr[d][XSTEP];
gridlen[d] = (int) loc->xgr[d][XLENGTH];
cumgridlen[d+1] = gridlen[d] * cumgridlen[d];
// printf("set d=%d inc=%f gridlen=%d cum=%d\n", d, inc[d], gridlen[d], (int) cumgridlen[d]);
}
}
// printf("maxstable=%d\n", maxstable);if (PL == 109) BUG;
if (maxstable) {
if (loggiven) {
for (i=0; i<vdimtot; i++) res[i] = RF_NEGINF;
}
else for (i=0; i<vdimtot; i++) res[i] = 0.0;
if (sub->mpp.moments < 1 || !R_FINITE(sub->mpp.mMplus[1]) ||
sub->mpp.mMplus[1] == 0.0)
ERR("integral of positive part of the submodel is unknown or not finite");
// M1 = sub->mpp.mMplus[1];
threshold = logthreshold = RF_INF;
pgs->globalmin = Minimum;
pgs->currentthreshold =
loggiven ? pgs->globalmin - threshold : pgs->globalmin / threshold;
if (simuxi != NULL) {
// to do : RPstudent, extremal t process
}
} else if (compoundpoisson){
assert(simuxi == NULL);
for (i=0; i<vdimtot; i++) res[i] = 0.0;
} else {
assert(simuxi == NULL);
for (i=0; i<vdimtot; res[i++] = 0.0);
logM2 = log(sub->mpp.mM[2]);
pgs->currentthreshold = GLOBAL.mpp.about_zero;
}
// APMI(sub);
for(n=0; ; n++) {
// if (n % 1000 == 0)
// printf("n=%d tot=%d contr=%d: %f < %f; res[0]=%f\n", (int) n, (int) Total_n, (int) control, res[control], threshold, res[0]);
// assert(n <= 10);
// printf("n=%ld\n", n);;if (PL == 110) BUG;
// assert(n < 1200);
// assert(res[0] < 1e-4);
// assert(res[0] < 80);
if (!maxstable && (n >= Total_n)) break;
// for (d=0; d<dim; d++) info->min[d] = info->max[d] = 0.0;
//assert(n < 2 );
// printf("XX.\n"); APMI(sub);
if (!sub->deterministic &&
(err = REINIT(sub, sub->mpp.moments, s)) != NOERROR) BUG;
DO(sub, s);
//printf("PL =%d\n", PL);;if (PL == 111) BUG;
// APMI(sub);
if (!sub->deterministic) SET_SUB;
// printf("here!!!\n");
// PMI(sub);
fieldreturn = sub->fieldreturn;
// MARCO: hattest Du dies reingenommen?
// if (subrole == ROLE_BROWNRESNICK) {
// assert(sub->SBR != NULL);
// n += sub->SBR->hatnumber;
// }
//printf("log_den %f\n", pgs->log_density);
logdens = pgs->log_density;
if (!R_FINITE(logdens)) {
// PMI(sub);
BUG;
// get_shape = CovList[sub->gatternr].cov;
// get_logshape = CovList[sub->gatternr].log;
}
if (compoundpoisson) {
summand = 0.0;
} else if (poissongauss) {
summand = -0.5 * (logdens + logM2);
} else { // max-stable field
assert(subrole == ROLE_SMITH || subrole == ROLE_SCHLATHER ||
subrole == ROLE_BROWNRESNICK);
if (n >= GLOBAL.extreme.maxpoints) {
PRINTF("'maxpoints' (= %d) exceeded. Simulation is likely to be a rough approximation only. Consider increasing 'maxpoints'.\n",
GLOBAL.extreme.maxpoints
);
break;
}
poisson += rexp(1.0);
frechet = pow(poisson, - 1.0 / pgs->alpha);
gumbel = log(frechet);
// MARCO, habe ich am 17.11.13 geaendert
// pgs->log_zhou_c = 0.0; // wrong! to do
summand = gumbel - logdens;
// threshold = sub->loggiven ? GumbelThreshold : FrechetThreshold;
threshold = logthreshold =
(double) (gumbel + log(sub->mpp.maxheights[0])); //Gumbel
//
if (!loggiven) {
threshold = exp(threshold); // Frechet
}
//
// APMI(sub);
// printf("Max = %f gumbel=%f %f %d logdens=%f\n", sub->mpp.maxheights[0], gumbel, threshold, loggiven, logdens);;if (PL == 112) BUG;
//printf("for thres=%e %d res=%e zhou=%f(%f) logden=%f gumbel=%f loggiven=%d summand=%f\n", threshold, (int) control, res[control], pgs->zhou_c, sub->mpp.maxheights[0], logdens, gumbel, loggiven, summand);// assert(false);
// assert(R_FINITE(threshold ));
// { int i; for (i=0; i<total_pts; i++) printf("%e ", res[i]); }
while ((control<total_pts) && (res[control]>=threshold)) {
// print("control: %d %f %f\n", control, res[control], threshold); assert(false);
control++;
}
if (control >= total_pts) {
break;
}
// printf("n=%d every =%d %d %d\n", (int) n, GLOBAL.extreme.check_every, PL, PL_STRUCTURE);
if (n % GLOBAL.extreme.check_every == 0) {
pgs->globalmin = res[control];
for (i=control + 1; i<total_pts; i++)
if (res[i] < pgs->globalmin) pgs->globalmin = res[i];
if (pgs->globalmin < Minimum) pgs->globalmin = Minimum;
if (PL >= PL_STRUCTURE)
PRINTF("control: %ld %f %f global=%f n=%ld every=%d %f log.max=%f\n",
control, res[control], threshold, pgs->globalmin,
n, GLOBAL.extreme.check_every, sub->mpp.maxheights[0],
log(sub->mpp.maxheights[0]));
}
// Marco:
pgs->currentthreshold = loggiven
? pgs->globalmin - gumbel
: pgs->globalmin / frechet;
//printf("loggiven =%d %e %f cur.thr=%e\n", loggiven, pgs->globalmin, gumbel, pgs->currentthreshold);;if (PL == 113) BUG;
} // maxstable
factor = exp(summand);
// printf("factor %4.4f sumd=%4.4f logdens=%4.4f logM2=%4.4f th=%4.4f %d<%d \n", factor, summand, logdens, logM2, threshold, (int) control, (int) total_pts); // assert(false);
// printf("dim=%d loggiven=%d maxstable=%d simugrid=%d randomcoin=%d\n", dim, loggiven, maxstable, simugrid, randomcoin);// assert(false);
//printf("A=%d\n", n);
zaehler = 0;
if (simugrid) {
partialfield = false;
for (d=0; d<dim; d++) {
double
step = inc[d] == 0.0 ? 1.0 : inc[d],
dummy = ceil((pgs->supportmin[d] - loc->xgr[d][XSTART]) / step
- 1e-8);
start[d] = dummy < 0 ? 0 : dummy > gridlen[d] ? gridlen[d] : (int)dummy;
partialfield |= start[d] > 0;
dummy = trunc(1.00000001 +
((pgs->supportmax[d] - loc->xgr[d][XSTART]) / step ));
end[d] = dummy < 0 ? 0 : dummy > gridlen[d] ? gridlen[d] : (int) dummy;
partialfield |= end[d] < gridlen[d];
// printf("window [%f %f] [%d %d] %d %d\n", pgs->supportmin[d], pgs->supportmax[d], start[d], end[d], gridlen[d], partialfield); //assert(n < 150); //APMI(cov);
if (start[d] >= end[d]) { //
// printf("broken %d %d %d supp=%f %f, loc.start=%f %f\n", d, start[d], end[d], pgs->supportmin[d], pgs->supportmax[d], loc->xgr[d][XSTART], inc[d]);
break;
}
delta[d] = (end[d] - start[d]) * cumgridlen[d];
nx[d] = start[d];
zaehler += start[d] * cumgridlen[d];
x[d] = xstart[d] = loc->xgr[d][XSTART] + (double) start[d] * inc[d];
if (false || zaehler < 0) {
PRINTF("d=%d start=%d, end=%d xstart=%f %f pgs=[%4.3f, %4.3f] xgr=%f %f %f inc=%3.2f\n",
d, start[d], end[d], xstart[d], pgs->xstart[d], pgs->supportmin[d], pgs->supportmax[d],
loc->xgr[d][XSTART], loc->xgr[d][XSTEP], loc->xgr[d][XLENGTH],
inc[d]);
// assert(false);
}
}
if (d < dim) {
//printf("cont: d=%d\n", d);
continue;
}
}
//
//printf("simugrid %d\n", simugrid);
// assert(false);
// APMI(cov);
#define STANDARDINKREMENT_ZAEHLER \
d = 0; \
nx[d]++; \
x[d] += inc[d]; \
zaehler += cumgridlen[d]; \
while (nx[d] >= end[d]) { \
nx[d] = start[d]; \
x[d] = xstart[d]; \
zaehler -= delta[d]; /* delta is positive */ \
if (++d >= dim) break; \
nx[d]++; \
x[d] += inc[d]; \
zaehler += cumgridlen[d]; \
} \
if (d >= dim) { break;}
// end define StandardInkrement
#define GET SHAPE(x, sub, value); /*printf("%f : %f fctr=%f %ld \n", x[0], *value, factor, zaehler); // */
#define LOGGET LOGSHAPE(x, sub, value, Sign);
#define LOGGETPOS LOGGET if (Sign[0] > 0)
#define GETRF *value = (double) (rf[zaehler]);
// printf("%d\n", total_pts); assert(false);
#define RFMAX(OP) { \
rf = sub->rf; \
for (i=0; i<total_pts; i++) { \
*value = (double) (OP); \
if (res[i] < *value) res[i]=*value; \
} \
}
#define RFSUM(OP) { \
rf = sub->rf; \
for (i=0; i<total_pts; i++) { \
*value = (double) (OP); \
if (res[i] < *value) res[i] += OP; \
} \
}
#define MAXAPPLY(GET, OP) \
assert(zaehler >=0 && zaehler<total_pts); \
/* if (zaehler<0 || zaehler>=total_pts) { PMI(cov); printf("//z=%d n=%d\n", zaehler, n); } */ \
GET { \
OP; \
/* assert( {if (*value > sub->mpp.maxheights[0] * exp(gumbel) * (1+1e-15)) { printf("//r=%f %f delta=%e %d %e\n", *value / exp(gumbel), sub->mpp.maxheights[0], *value / exp(gumbel) - sub->mpp.maxheights[0], loggiven, factor); };true;}); */ \
/*assert(loggiven || *value <= sub->mpp.maxheights[0] * exp(gumbel) * (1+1e-15)); */ \
if (res[zaehler] < *value) res[zaehler] = *value; \
}
#define SUMAPPLY(GET,OP) GET; res[zaehler] += OP;
#define APPLY(GET, OP, WHAT) { \
if (simugrid) { \
while(true) { \
WHAT(GET, OP); \
STANDARDINKREMENT_ZAEHLER; \
} \
} else { \
for(x=loc->x; zaehler<spatial; zaehler++, x += dim) { \
WHAT(GET, OP); \
} \
} \
}
#define OG_BASIC(OP) \
assert(zaehler >=0 && zaehler < total_pts); \
int jj, idx, index=0; \
for (jj=0; jj<dim; jj++) { \
idx = round((x[jj] - ogstart[jj]) / ogstep[jj]); \
if (idx < 0 || idx >= pgs->own_grid_len[jj]) break; \
index += pgs->own_grid_cumsum[jj] * idx; \
} \
if (jj >= dim) { OP; } else if (false) printf("%d %f %f %f len=%f %d %d\n", jj, x[jj], ogstart[jj], ogstep[jj], pgs->own_grid_len[jj], pgs->own_grid_cumsum[jj], idx); /* // */ \
#define OG_MAXAPPLY(NONE, OP) \
rf = sub->rf; \
OG_BASIC(*value = rf[index]; \
OP; if (res[zaehler] < *value) res[zaehler]=*value;)
#define OG_SUMAPPLY(NONE, OP) \
rf = sub->rf; \
OG_BASIC(*value = rf[index]; res[zaehler] += OP;)
// OWNGRIDAPPLY(OP, OG_SUMAPPLY_BASIC)
#define OWNGRIDAPPLY(NONE, OP, WHAT) { \
double *ogstart = pgs->own_grid_start, \
*ogstep = pgs->own_grid_step; \
APPLY(NONE, OP, WHAT); \
}
#define ALL_APPLY(APPLY, MGET, MAPPLY, SGET, SAPPLY, NONLOGGET) \
if (loggiven == true) { \
if (maxstable) APPLY(MGET, (*value) += summand, MAPPLY) \
else APPLY(SGET, Sign[0] * exp(*value + summand), SAPPLY); \
} else { /* not loggiven */ \
if (sub->loggiven){ \
if (maxstable) \
APPLY(MGET, *value=Sign[0] * exp(*value+summand), MAPPLY) \
else APPLY(SGET, Sign[0] * exp(*value + summand), SAPPLY); \
} else { \
assert(R_FINITE(factor)); \
if (maxstable) APPLY(NONLOGGET, (*value) *= factor, MAPPLY) \
else APPLY(NONLOGGET, *value * factor, SAPPLY); \
} \
}
#define AVEAPPLY(G,OP0,OP1) { \
warning(" * timescale ?!"); \
if (simugrid) { \
while (true) { \
long zt; \
double zw,inct, segt, ct, st, A, cit, sit; \
inct = sub->q[AVERAGE_YFREQ] * loc->xgr[dim][XSTEP]; \
cit = cos(inct); \
sit = sin(inct); \
G; \
segt = OP1; \
A = OP0; \
ct = A * cos(segt); \
st = A * sin(segt); \
for (zt = zaehler; zt < total_pts; zt += spatial) { \
res[zt] += (double) ct; \
zw = ct * cit - st * sit; \
st = ct * sit + st * cit; \
ct = zw; \
res[zaehler] += zw; \
} \
STANDARDINKREMENT_ZAEHLER; \
} \
} else { /* not a grid */ \
x=loc->x; \
/* the following algorithm can greatly be improved ! */ \
/* but for ease, just the stupid algorithm */ \
for (; zaehler<spatial; zaehler++, x += dim) { \
G; res[zaehler] += OP0 * cos(OP1); \
} \
} \
}
//printf("%d %d %d\n", fieldreturn, loggiven, sub->loggiven);
// todo : !simugrid! && Bereiche einengen!!
// *atomdependent rfbased/not; (( recurrent/dissipative: min/max gesteuert))
if (poissongauss && !randomcoin) { // AVERAGE
if (sub->loggiven) {
AVEAPPLY(LOGGET, Sign[0] * exp(value[0] + summand),
Sign[1] * exp(value[1]));
} else AVEAPPLY(GET, value[0] * factor, value[1]);
} else {
assert(subrole == ROLE_SMITH || subrole == ROLE_SCHLATHER ||
subrole ==ROLE_BROWNRESNICK || subrole==ROLE_POISSON ||
subrole==ROLE_POISSON_GAUSS);
if (fieldreturn == HUETCHEN_OWN_GRIDSIZE) { // atomdependent rfbased/not;
ALL_APPLY(OWNGRIDAPPLY, NULL, OG_MAXAPPLY, NULL, OG_SUMAPPLY, NULL);
} else if (fieldreturn == true) { // extended boolean
// todo : !simugrid! && Bereiche einengen!!
// note ! no Sign may be given currently when fieldreturn and loggiven!!
if (partialfield) {
//!! Windows
if (!simugrid) BUG;
ALL_APPLY(APPLY, GETRF, MAXAPPLY, GETRF, SUMAPPLY, GETRF);
} else { // !partialfield
if (sub->loggiven) {
if (maxstable) RFMAX(rf[i] + summand)
else RFSUM(exp(rf[i] + summand));
} else { // Linux
if (maxstable) RFMAX(rf[i] * factor)
else RFSUM(rf[i] * factor);
}
}
} else if (fieldreturn == false) { // not field return
ALL_APPLY(APPLY, LOGGETPOS, MAXAPPLY, LOGGET, SUMAPPLY, GET);
} else BUG; // neither HUETCHEN_OWN_GRIDSIZE, nor true nor false
}
if (n >= nthreshold) {
if (maxstable) {
LPRINT("%d-th pos.: value=%1.3f threshold=%1.3f gu=%1.3f %1.3f %d (%d %d)\n",
control, (double) res[control], (double) threshold,
gumbel, logdens, loggiven,
n, nthreshold); //, deltathresh);
nthreshold += deltathresh;
} else {
PRINTF("n=%d Total=%d\n", (int) (n / 1000), (int)(Total_n));
}
}
R_CheckUserInterrupt();
// printf("n=%d %f %f %f 0:%f\n", (int) n, res[0], res[1], res[2], res[3]);;if (PL == n) BUG;
// assert(res[0] < 80);
}// for(n,... -- while control < total_pts
//printf("n=%d %f %f %f 3::%f\n", n, res[0], res[1], res[2], res[3]);
// Trafo to margins
// printf("%d frechet\n", sub->loggiven);if (PL == 10000) BUG;
if (maxstable){
double meansq, var,
mean = RF_NA,
n_min = RF_INF,
eps = GLOBAL.extreme.eps_zhou;
if (pgs->alpha != 1.0) {
if (loggiven) for (k=0; k<total_pts; k++) res[k] *= pgs->alpha;
else for (k=0; k<total_pts; k++) res[k] = pow(res[k], pgs->alpha);
}
if (pgs->estimated_zhou_c) {
if (PL > 5) PRINTF("zhou_c: %ld %d",pgs->n_zhou_c, GLOBAL.extreme.min_n_zhou);
while (pgs->n_zhou_c < GLOBAL.extreme.min_n_zhou) {
if ((err = REINIT(sub, sub->mpp.moments, s)) != NOERROR) BUG;
DO(sub, s);
SET_SUB;
}
while (pgs->n_zhou_c < n_min) {
for (k=0; k<10; k++) {
if ((err = REINIT(sub, sub->mpp.moments, s)) != NOERROR) BUG;
DO(sub, s);
SET_SUB;
}
mean = (double) (pgs->sum_zhou_c / (long double) pgs->n_zhou_c);
meansq = mean * mean;
var = (double) (pgs->sq_zhou_c / (long double) pgs->n_zhou_c
- meansq);
n_min = var / (meansq * eps * eps);
if (n_min > GLOBAL.extreme.max_n_zhou) n_min =GLOBAL.extreme.max_n_zhou;
//
if (PL > 5) PRINTF(" n=%ld, min=%f %f mean=%f zhou=%f eps=%f\n", pgs->n_zhou_c, n_min, (double) pgs->sum_zhou_c, mean, pgs->zhou_c, eps);
}
//
} else {
mean = pgs->zhou_c;
}
// mean /= M1;
//printf(" n=%ld, min=%f %f mean=%f zhou=%f eps=%f %f max=%f\n", pgs->n_zhou_c, n_min, (double) pgs->sum_zhou_c, mean, pgs->zhou_c, eps, M1, sub->mpp.maxheights[0]);
// printf("loggiven = %d %f %f\n", loggiven, mean, pgs->zhou_c);
// mean = 250; // 2 statt 11
double xi = P0(GEV_XI);
// printf("xi=%f %d logmean=%d %f\n", xi, loggiven, pgs->logmean, mean);
if (loggiven) { // 3.6.13 sub->loggiven geaendert zu loggiven
if (!pgs->logmean) mean = log(mean);
for (k=0; k<total_pts; k++) res[k] += mean;
if (xi != 0.0)
for (k=0; k<total_pts; k++) res[k] = exp(xi * res[k]) - 1.0;
} else {
for (k=0; k<total_pts; k++) res[k] *= mean;
if (xi==0.0) for (k=0; k<total_pts; k++) res[k] = log(res[k]);
else if (xi == 1.0) for (k=0; k<total_pts; k++) res[k] -= 1.0;
else for (k=0; k<total_pts; k++) res[k] = pow(res[k], xi) - 1.0;
}
// muss das vorletzte sein !
if (cov->kappasub[GEV_S] != NULL) {
ERR1("'%s' cannot be chosen as a function yet.", KNAME(GEV_S));
// check whether s is positive !!
} else {
double ss = xi == 0 ? P0(GEV_S) : (P0(GEV_S) / xi);
//printf("s=%f\n", ss);
if (ss != 1.0) for (k=0; k<total_pts; k++) res[k] *= ss;
}
// muss das letzte sein !
if (cov->kappasub[GEV_MU] != NULL) {
ERR1("'%s' cannot be chosen as a function yet.", KNAME(GEV_MU));
} else {
double mu = P0(GEV_MU);
// printf("mu=%f xi=%f %f\n", mu, xi, s);
//printf("mu=%f\n", mu);
if (mu != 0.0) for (k=0; k<total_pts; k++) res[k] += mu;
}
}
//printf("n=%d %f %f %f 0:%f\n", n, res[0], res[1], res[2], res[3]);;if (PL == 116) BUG;
// for (k=0; k<total_pts; k++) printf("%f ", res[k]); printf("\n");
if (PL >= PL_DETAILSUSER) {
PRINTF("number of shape functions used: %ld\n", n);
}
return;
} // end dompp
void dompp(cov_model *cov, gen_storage *s) {
dompp(cov, s, NULL);
}
//void calculate_integral(int d, double R, cov_model *cov,
// integral_type *integral) {
//
//}
////////////////////////////////////////////////////////////////////
// Poisson
int check_poisson(cov_model *cov) {
cov_model *next=cov->sub[0],
*key = cov->key,
*sub = key == NULL ? next : key;
int err,
dim = cov->tsdim; // taken[MAX DIM],
mpp_param *gp = &(GLOBAL.mpp);
Types type = sub == key ? PointShapeType : ShapeType;
//print("eee\n");
cov->role = ROLE_POISSON;
kdefault(cov, POISSON_INTENSITY, gp->intensity[dim]);
if ((err = checkkappas(cov, true)) != NOERROR) return err;
if (cov->tsdim != cov->xdimprev || cov->tsdim != cov->xdimown)
return ERRORDIM;
if ((err = CHECK(sub, dim, dim, type, XONLY, CoordinateSystemOf(cov->isoown),
SUBMODEL_DEP, cov->role)) != NOERROR) return err;
setbackward(cov, sub);
return NOERROR;
}
void range_poisson(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
range->min[POISSON_INTENSITY] = RF_NEGINF;
range->max[POISSON_INTENSITY] = RF_INF;
range->pmin[POISSON_INTENSITY] = -10;
range->pmax[POISSON_INTENSITY] = 10;
range->openmin[POISSON_INTENSITY] = true;
range->openmax[POISSON_INTENSITY] = true;
}
int struct_poisson(cov_model *cov, cov_model **newmodel){
cov_model *next=cov->sub[0];
location_type *loc = Loc(cov);
ASSERT_NEWMODEL_NULL;
if (cov->role != ROLE_POISSON)
SERR1("'%s' not called as random coin", NICK(cov));
if (cov->key != NULL) COV_DELETE(&(cov->key));
if (loc->Time || (loc->grid && loc->caniso != NULL)) {
TransformLoc(cov, false, GRIDEXPAND_AVOID, false);
SetLoc2NewLoc(next, PLoc(cov)); // passt das?
}
if (!isPointShape(next)) {
int err;
if ((err = covCpy(&(cov->key), next)) != NOERROR) return err;
if ((err = addStandard(&(cov->key))) != NOERROR) return err;
}
return NOERROR;
}
int init_poisson(cov_model *cov, gen_storage *S) {
// location_type *loc = Loc(cov);
// mpp_storage *s = &(S->mpp);
cov_model *key=cov->key;
int err;
if ((err = init_mpp(cov, S)) != NOERROR) {
return err;
}
if (!isPointShape(key)) SERR("no definition of a shape function found");
key->Spgs->intensity = key->Spgs->totalmass * P0(POISSON_INTENSITY);
return err;
}
////////////////////////////////////////////////////////////////////
// Schlather
void extremalgaussian(double *x, cov_model *cov, double *v) {
// schlather process
cov_model *next = cov->sub[0];
if (cov->role == ROLE_SCHLATHER) COV(x, next, v)
else {
COV(x, next, v);
*v = 1.0 - sqrt(0.5 * (1.0 - *v));
}
}
//#define MPP_SHAPE 0
//#define MPP_TCF 1
int SetGEVetc(cov_model *cov, int role) {
int err;
extremes_param *ep = &(GLOBAL.extreme);
if (cov->role != ROLE_COV) cov->role = role;
if (cov->sub[MPP_TCF] != NULL && cov->sub[MPP_SHAPE]!=NULL)
SERR2("either '%s' or '%s' must be given",
SNAME(MPP_TCF), SNAME(MPP_SHAPE));
if ((err = checkkappas(cov, false)) != NOERROR) return err;
kdefault(cov, GEV_XI, ep->GEV_xi);
kdefault(cov, GEV_S, P0(GEV_XI) == 0.0 ? 1.0 : fabs(P0(GEV_XI)));
kdefault(cov, GEV_MU, P0(GEV_XI) == 0.0 ? 0.0 : 1.0);
// print("xi s mu %f %f %f\n", P0(GEV_XI), P0(GEV_S), P0(GEV_MU)); assert(false);
if ((err = checkkappas(cov, true)) != NOERROR) return err;
if (cov->tsdim != cov->xdimprev || cov->tsdim != cov->xdimown) {
return ERRORDIM;
}
return NOERROR;
}
int check_schlather(cov_model *cov) {
//print("check schlather\n");
if (!((cov->sub[MPP_TCF] == NULL) xor (cov->sub[MPP_SHAPE] == NULL)))
SERR("two submodels given instead of one.");
cov_model
*key = cov->key,
*next = cov->sub[cov->sub[MPP_TCF] != NULL ? MPP_TCF : MPP_SHAPE];
int err,
dim = cov->tsdim; // taken[MAX DIM],
// mpp_param *gp = &(GLOBAL.mpp);
double v;
bool schlather = CovList[cov->nr].Init == init_mpp;
// printf("A \n");
ASSERT_ONE_SUBMODEL(cov);
/// printf("B A \n");
if ((err = SetGEVetc(cov, ROLE_SCHLATHER)) != NOERROR) return err;
// printf("CA \n");
cov_model *sub = cov->key==NULL ? next : key;
// printf("GA \n");
if (key == NULL) {
// printf("key=NULL\n");
int role = isVariogram(sub) ? ROLE_COV //Marco, 29.5.13
: isShape(sub) && schlather ? ROLE_SCHLATHER
: isGaussProcess(sub) ? ROLE_GAUSS
: isBernoulliProcess(sub) && schlather ? ROLE_BERNOULLI
: ROLE_UNDEFINED;
ASSERT_ROLE_DEFINED(sub);
if ((err = isPosDef(next)
? CHECK(next, dim, dim, PosDefType, XONLY, //err
IsotropicOf(cov->isoown), SCALAR, role)
: CHECK(next, dim, dim, ProcessType, XONLY, //err,
CoordinateSystemOf(cov->isoown),
SCALAR, role)) != NOERROR) return err;
if (sub->vdim[0] != 1) SERR("only univariate processes are allowed");
assert(cov->vdim[0] == cov->vdim[1]);
setbackward(cov, sub);
if (role==ROLE_COV) {
if (next->pref[Nothing] == PREF_NONE) return ERRORPREFNONE;
COV(ZERO, next, &v);
if (v != 1.0)
SERR2("a correlation function is required as submodel, but '%s' has variance %f.", NICK(next), v);
}
} else { // key != NULL
//printf("key!=NULL\n");
if ( (err = CHECK(key, dim, dim, PointShapeType, XONLY,
CoordinateSystemOf(cov->isoown),
SUBMODEL_DEP, ROLE_SCHLATHER)) != NOERROR) {
//APMI(cov);
return err;
}
setbackward(cov, sub);
}
//print("end check schlather\n");
return NOERROR;
}
int struct_schlather(cov_model *cov, cov_model **newmodel){
cov_model
*sub = cov->sub[cov->sub[MPP_TCF] != NULL ? MPP_TCF : MPP_SHAPE];
int err, role, ErrNoInit;
bool schlather = CovList[cov->nr].Init == init_mpp; // else opitz
if (cov->role != ROLE_SCHLATHER) BUG;
ASSERT_NEWMODEL_NULL;
if (cov->key != NULL) COV_DELETE(&(cov->key));
if (cov->sub[MPP_TCF] != NULL) {
if ((err = STRUCT(sub, &(cov->key))) > NOERROR) return err;
cov->key->calling = cov;
} else {
if ((err = covCpy(&(cov->key), sub)) != NOERROR) return err;
}
assert(cov->key->calling == cov);
if (cov->key->nr != GAUSSPROC && !isBernoulliProcess(cov->key)) {
if (isVariogram(cov->key)) addModel(&(cov->key), GAUSSPROC);
else {
if (isGaussProcess(cov->key)) {
SERR("invalid model specification");
} else SERR2("'%s' currently only allowed for gaussian processes %s",
NICK(cov), schlather ? "and binary gaussian processes" : "");
}
}
assert(cov->key->calling == cov);
role = cov->key->nr == GAUSSPROC ? ROLE_GAUSS //Marco, 29.5.13
: isBernoulliProcess(cov->key) ? ROLE_BERNOULLI
: ROLE_UNDEFINED;
//PMI(cov->key, "role ");
ASSERT_ROLE_DEFINED(cov->key);
if ((err = CHECK(cov->key, cov->tsdim, cov->xdimown, ProcessType,
cov->domown,
cov->isoown, cov->vdim, role)) != NOERROR) return err;
// APMI(cov->key);
if ((ErrNoInit = STRUCT(cov->key, NULL)) > NOERROR) return ErrNoInit; //err
addModel(&(cov->key), STATIONARY_SHAPE);
if ((err = CHECK(cov->key, cov->tsdim, cov->xdimown, PointShapeType,
cov->domown,
cov->isoown, cov->vdim, ROLE_SCHLATHER)) != NOERROR)
return err;
return ErrNoInit;
}
#define OPITZ_ALPHA (LAST_MAXSTABLE + 1)
int init_opitzprocess(cov_model *cov, gen_storage *S) {
int err;
if ( (err = init_mpp(cov, S)) != NOERROR) return err;
cov_model *key = cov->key;
assert(key != NULL);
assert(key->mpp.moments == 1);
pgs_storage *pgs = key->Spgs;
assert(pgs != NULL);
key->mpp.mMplus[1] = INVSQRTTWOPI * pow(2, 0.5 * P0(OPITZ_ALPHA) - 0.5) *
gammafn(0.5 * P0(OPITZ_ALPHA) + 0.5);
pgs->zhou_c = 1.0 / key->mpp.mMplus[1];
pgs->alpha = P0(OPITZ_ALPHA);
return NOERROR;
}
void range_opitz(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
// assert(cov != NULL);
range_mpp(cov, range);
range->min[OPITZ_ALPHA] = 0;
range->max[OPITZ_ALPHA] = RF_INF;
range->pmin[OPITZ_ALPHA] = 0.1;
range->pmax[OPITZ_ALPHA] = 5;
range->openmin[OPITZ_ALPHA] = true;
range->openmax[OPITZ_ALPHA] = true;
}
////////////////////////////////////////////////////////////////////
// Smith
// siegburg 7:47, 7:54; 8:19, Bus 611 8:32; 8:41 Immenburg
void param_set_identical(cov_model *to, cov_model *from, int depth) {
assert(depth >= 0);
assert(to != NULL && from != NULL);
assert(strcmp(NICK(from), NICK(to)) == 0); // minimal check
int i;
assert(from->qlen == to->qlen);
if (from->q != NULL) MEMCOPY(to->q, from->q, sizeof(double) * from->qlen);
for (i=0; i<MAXPARAM; i++) { PCOPY(to, from, i); }
if (depth > 0) {
for (i=0; i<MAXSUB; i++) {
if (from->sub[i] != NULL)
param_set_identical(to->sub[i], from->sub[i], depth - 1);
}
}
}
int PointShapeLocations(cov_model *key, cov_model *shape) {
// fuegt im max-stabilen Fall diegleiche Funktion als
// intensitaetsfunktion fuer die Locationen ein; random scale
// wird mitbeachtet
//
// fuegt im coin fall die ge-powerte Shape-Fkt als
// intensitaetsfunktion fuer die Locationen ein; random scale
// wird mitbeachtet
//
// fuegt
assert(key != NULL);
// printf("fx \n");
int err, i,
dim = key->xdimown,
nr = key->nr;
assert(isPointShape(key));
if (key->sub[PGS_LOC] != NULL) return NOERROR;
if ((err = covCpy(key->sub + PGS_FCT, shape)) != NOERROR) return err;
if (nr == PTS_GIVEN_SHAPE) {
//printf("here %ld %d %d %d\n", key->sub[PGS_LOC], ScaleOnly(shape),
// !shape->deterministic, shape->sub[0]->deterministic);
key->nsub = 2;
if (key->sub[PGS_LOC] != NULL) BUG;
bool randomscale = ScaleOnly(shape) && !shape->deterministic &&
shape->sub[0]->deterministic;
if ((err = covCpyWithoutRandomParam(key->sub + PGS_LOC,
randomscale ? shape->sub[0] : shape))
!= NOERROR) return err;
if (shape->role == ROLE_POISSON_GAUSS) {
assert(dim <= MAXMPPDIM);
addModel(key, PGS_LOC, POW);
assert(key->nsub == 2);
kdefault(key->sub[PGS_LOC], POW_ALPHA, GLOBAL.mpp.shape_power);
addModel(key, PGS_LOC, SCATTER);
PARAMALLOC(key, SCATTER_MAX, dim, 1);
for (i=0; i<dim; i++)
PARAM(key, SCATTER_MAX)[i] = GLOBAL.mpp.scatter_max[i];
PARAMALLOC(key, SCATTER_STEP, dim, 1);
for (i=0; i<dim; i++)
PARAM(key, SCATTER_STEP)[i] = GLOBAL.mpp.scatter_step[i];
addModel(key, PGS_FCT, RANDOMSIGN);
} else if (shape->role != ROLE_MAXSTABLE) BUG;
if (!randomscale && !shape->deterministic) {
addSetDistr(key->sub + PGS_LOC, key->sub[PGS_FCT],
param_set_identical, true, MAXINT);
}
addModel(key, PGS_LOC, RECTANGULAR);
if (randomscale) {
addModel(key, PGS_LOC, LOC);
addSetDistr(key->sub + PGS_LOC, shape, ScaleDollarToLoc, true, 0);
}
} else if (nr == STANDARD_SHAPE) {
assert(key != NULL && shape != NULL);
if ((err = STRUCT(shape, key->sub + PGS_LOC)) != NOERROR) return err;
key->sub[PGS_LOC]->calling = key;
if (key->sub[PGS_LOC] == NULL)
SERR("simple enlarging of the simulation window does not work");
} else BUG;
// APMI(key);
return NOERROR;
}
int addPointShape(cov_model **Key, cov_model *shape, cov_model *pts,
cov_model *local_pts,
int dim, int vdim) {
// versucht die automatische Anpassung einer PointShapeType-Funktion;
// derzeit wird
// * PGS und
// * STANDARD_SHAPE (weiss nicht mehr wieso -> coins?)
// versucht; inklusive Init
#define PGS_N 2
int i,
err = NOERROR,
pgs[PGS_N] = {PTS_GIVEN_SHAPE, STANDARD_SHAPE};
#define msgN 200
char msg[PGS_N][LENERRMSG];
assert(shape != NULL);
assert(*Key == NULL);
// to do: strokorbball: raeumlich waere Standard/Uniform besser;
// praeferenzen programmieren?
for (i=0; i<PGS_N; i++) {
if (i > 0) {
FinalErrorMSG(err, msg[i-1]);
XERR(err);
}
// if (i > 0) XERR(err); assert(i ==0);
if (*Key != NULL) COV_DELETE(Key);
assert(shape->calling != NULL);
addModel(Key, pgs[i], shape->calling);
assert((*Key)->sub[PGS_LOC] == NULL && (*Key)->sub[PGS_FCT] == NULL);
// PMI(shape);
// PMI(pts);
if (pts != NULL) {
if ((err = covCpy((*Key)->sub + PGS_FCT, shape)) != NOERROR ||
(err = covCpy((*Key)->sub + PGS_LOC, pts)) != NOERROR
) return err;
Ssetcpy((*Key)->sub[PGS_FCT], (*Key)->sub[PGS_LOC], shape, pts);
Ssetcpy((*Key)->sub[PGS_LOC], (*Key)->sub[PGS_FCT], pts, shape);
// PMI((*Key)->sub[PGS_FCT]);
// APMI((*Key)->sub[PGS_LOC]);
} else {
//
if ((err = PointShapeLocations(*Key, shape)) != NOERROR) continue;
if (local_pts != NULL) { // only in plusmalS.cc
if ((*Key)->nr != PTS_GIVEN_SHAPE) continue;
cov_model *local, *dummy;
if ((err = covCpy(&local, false, local_pts, shape->prevloc, NULL,
true, true, false)) != NOERROR) return err;
local->calling = (*Key)->calling;
dummy = local;
while (dummy->sub[DOLLAR_SUB] != NULL) dummy = dummy->sub[DOLLAR_SUB];
if (dummy->nr != LOC) BUG;
dummy->sub[DOLLAR_SUB] = *Key;
(*Key)->calling = dummy;
}
}
(*Key)->calling = shape->calling;
(*Key)->sub[PGS_FCT]->calling = (*Key)->sub[PGS_LOC]->calling = *Key;
//PMI((*Key)->calling);
assert((*Key)->sub[PGS_LOC] != NULL && (*Key)->sub[PGS_FCT] != NULL);
(*Key)->nsub = 2;
// APMI(*Key);
// printf(">>>>>>> KEY start %d %s \n", i, NICK((*Key)));
if ((err = CHECK(*Key, dim, dim, PointShapeType, XONLY,
CoordinateSystemOf(shape->isoown),
vdim, ROLE_MAXSTABLE)) != NOERROR) {
// printf(">>>>>>> KEY break %d %s %d\n", i, Nick(*Key), dim);
//PMI(*Key);
XERR(err);
continue;
}
NEW_COV_STORAGE(*Key, gen);
if ((err = INIT(*Key, 1, (*Key)->Sgen)) == NOERROR) break;
} // for i_pgs
if (err != NOERROR) {
errorMSG(err, msg[i-1]);
sprintf(ERRORSTRING,
"no method found to simulate the given model:\n\tpgs: %s\n\tstandard: %s)",
msg[0], msg[1]);
return ERRORM;
}
return NOERROR;
}
int addPointShape(cov_model **Key, cov_model *shape, cov_model *pts,
int dim, int vdim) {
return addPointShape(Key, shape, pts, NULL, dim, vdim);
}
int struct_ppp_pts(cov_model **Key, cov_model *shape,
cov_model *calling,
int tsdim, int vdim) {
// Ruft S TRUCT(shape, Key) auf und verarbeitet weiter, je nachdem
// was nun in Key steht: pgs, shape, random, nichts.
cov_model
*dummy = NULL;
int err = NOERROR;
if ((err = STRUCT(shape, Key)) == NOERROR && *Key != NULL) {
// todo: IN WELCHEM FALL laeuft man hier rein?
(*Key)->calling = calling;
Types type = TypeConsistency(PointShapeType, *Key, 0) ? PointShapeType :
TypeConsistency(RandomType, *Key, 0) ? RandomType :
TypeConsistency(ShapeType, *Key, 0) ? ShapeType :
OtherType;
if (type == RandomType) { // random locations given;
// so, it must be of pgs type (or standard):
if ((err = CHECK_R(*Key, shape->tsdim)) != NOERROR) {
goto ErrorHandling;
}
dummy = *Key;
*Key = NULL;
if ((err = addPointShape(Key, shape, dummy, tsdim, vdim)) != NOERROR)
goto ErrorHandling;
if (*Key == NULL) BUG;
(*Key)->calling = calling;
// APMI(cov);
} else {
if ((err = CHECK(*Key, shape->tsdim, shape->xdimprev, type,
shape->domprev, shape->isoprev, shape->vdim,
ROLE_MAXSTABLE))
!= NOERROR) {
goto ErrorHandling;
}
if (type==PointShapeType) {
if ((err = PointShapeLocations(*Key, shape)) != NOERROR)
goto ErrorHandling;
} else if (type == ShapeType) {
dummy = *Key;
*Key = NULL;
// suche nach geeigneten locationen
if ((err = addPointShape(Key, dummy, NULL, tsdim, vdim)) != NOERROR)
goto ErrorHandling;
//printf("e ddd\n");
} else BUG;
} // ! randomtype
} else { // S-TRUCT does not return anything
int err2;
//assert(false);
// XERR(err);
// APMI(*Key);
if ((err2 = addPointShape(Key, shape, NULL, tsdim, vdim)) != NOERROR) {
if (err == NOERROR) err = err2;
goto ErrorHandling;
}
err = NOERROR;
}
ErrorHandling:
if (dummy != NULL) COV_DELETE(&dummy);
return err;
}
int check_smith(cov_model *cov) {
cov_model
*shape = cov->sub[MPP_SHAPE],
*TCF = cov->sub[MPP_TCF],
*next = shape != NULL ? shape : TCF,
*key = cov->key,
*sub = (key==NULL) ? next : key;
;
int err, role,
dim = cov->tsdim; // taken[MAX DIM],
//Types type;
ASSERT_ONE_SUBMODEL(cov);
if ((err = SetGEVetc(cov, ROLE_SMITH)) != NOERROR) return err;
if (cov->tsdim != cov->xdimprev || cov->tsdim !=cov->xdimown) return ERRORDIM;
if (key != NULL) {
if ((err = CHECK(key, dim, dim, PointShapeType, XONLY,
CoordinateSystemOf(cov->isoown),
SUBMODEL_DEP, ROLE_SMITH)) != NOERROR) return err;
} else { // key == NULL
if (next == TCF) {
if ((err = CHECK(next, dim, dim, TcfType, XONLY, ISOTROPIC,
SCALAR, ROLE_SMITH)) != NOERROR) {
return err;
}
if ((dim == 1 && next->rese_derivs < 1) ||
(dim >= 2 && dim <= 3 && next->rese_derivs < 2) ||
dim > 3) {
//printf("dim = %d\n", dim);
//APMI(next);
SERR("submodel does not have enough derivatives (programmed).");
}
} else { // shape
role = // (isVariogram(sub) && !isPosDef(sub)) ? ROLE_COV //Marco, 29.5.13
isShape(sub) ? ROLE_MAXSTABLE // 23.4.14
: isPointShape(sub) ? ROLE_SMITH
: isGaussProcess(sub) ? ROLE_GAUSS
: isBernoulliProcess(sub) ? ROLE_BERNOULLI
: ROLE_UNDEFINED;
ASSERT_ROLE_DEFINED(sub);
if ((err = CHECK(next, dim, dim, ShapeType, XONLY,
CoordinateSystemOf(cov->isoown),
SCALAR, role)) != NOERROR) {
return err;
}
if (next->full_derivs < 0)
SERR1("'%s' requires an explicit submodel.", NICK(cov));
}
}
setbackward(cov, next);
return NOERROR;
}
int struct_smith(cov_model *cov, cov_model **newmodel){
cov_model
*tmp_shape = NULL,
*shape = cov->sub[MPP_SHAPE],
*tcf = cov->sub[MPP_TCF],
*tcf_shape = NULL,
*sub = shape != NULL ? shape : tcf;
location_type *loc = Loc(cov);
int err = NOERROR;
//APMI(shape);
if (cov->role != ROLE_SMITH) BUG;
if (loc->Time || (loc->grid && loc->caniso != NULL)) {
TransformLoc(cov, false, GRIDEXPAND_AVOID, false);
SetLoc2NewLoc(sub, PLoc(cov));
}
if (cov->key != NULL) COV_DELETE(&(cov->key));
ASSERT_NEWMODEL_NULL;
if (tcf != NULL) {
// to do: ausbauen
if ((err = covCpy(&tcf_shape, sub)) != NOERROR) goto ErrorHandling;
addModel(&tcf_shape, STROKORB_MONO);
if ((err = CHECK(tcf_shape, tcf->tsdim, tcf->xdimprev, ShapeType,
tcf->domprev, tcf->isoprev, tcf->vdim,
ROLE_MAXSTABLE)) != NOERROR) goto ErrorHandling;
tmp_shape = tcf_shape;
} else {
tmp_shape = shape;
}
// APMI(tmp_shape);
if ((err = struct_ppp_pts(&(cov->key), tmp_shape, cov,
cov->tsdim, cov->vdim[0]))
!= NOERROR) goto ErrorHandling;
// APMI(cov);
err = NOERROR;
ErrorHandling:
if (tcf_shape != NULL && tmp_shape != NULL) COV_DELETE(&tmp_shape);
return err;
}
typedef double (*logDfct)(double *data, double gamma);
double HueslerReisslogD(double *data, double gamma) {
double
g = sqrt(2.0 * gamma),
logy2y1 = log(data[1] / data[0]);
return -pnorm(0.5 * g + logy2y1 / g, 0.0, 1.0, true, false) / data[0]
-pnorm(0.5 * g - logy2y1 / g, 0.0, 1.0, true, false) / data[1];
}
double schlatherlogD(double *data, double gamma) {
double
sum = data[0] + data[1],
prod = data[0] * data[1];
return 0.5 * sum / prod *
(1.0 + sqrt(1.0 - 2.0 * (2.0 - gamma) * prod / (sum * sum)));
}
#define LL_AUTO 0
#define LL_FULL 1
#define LL_COMPOSITE 2
#define LL_SELECTION 3
#define LL_UNKNOWN ERR("unknown value for 'likelihood' -- please contact author");
void loglikelihoodMaxstable(double *data, cov_model *cov, logDfct logD,
double *v) {
cov_model *sub = cov;
while (isProcess(sub))
sub = sub->key != NULL ? sub->key : sub->sub[0];
if (cov->q == NULL) {
location_type *loc = Loc(cov);
long len = loc->totalpoints;
QALLOC(len);
if (loc->grid || loc->Time) TransformLoc(sub, false, true, false);
}
location_type *loc = Loc(cov);
double *x, *y;
int
dim = cov->xdimown;
long i, j,
len = loc->totalpoints;
if (data != NULL) {
double
xi = P0(GEV_XI),
mu = P0(GEV_MU),
s = P0(GEV_S);
if (xi == 0) {
for (i=0; i<len; i++) cov->q[i] = exp((data[i] - mu) / s);
} else {
double
xiinv = 1.0 / xi,
xis = xi / s;
for (i=0; i<len; i++) cov->q[i] = pow(1.0 + xis * (data[i] - mu), xiinv);
}
}
switch(GLOBAL.fit.likelihood) {
case LL_AUTO: case LL_COMPOSITE:
double z[2], gamma, gamma0;
COV(ZERO, sub, &gamma0);
x = loc->x;
for (i=0; i<len; i++, x += dim) {
y = x + dim;
for (j=i+1; j<len; j++) {
NONSTATCOV(x, y, sub, &gamma);
z[0] = cov->q[i];
z[1] = cov->q[j];
*v += logD(z, gamma0 - gamma);
}
}
break;
case LL_FULL: ERR("full MLE not available for Brown Resnick");
break;
case LL_SELECTION: NotProgrammedYet("LL_SELECTION");
default : LL_UNKNOWN;
}
}
void loglikelihoodBR(double *data, cov_model *cov, double *v) {
loglikelihoodMaxstable(data, cov, HueslerReisslogD, v);
}
void loglikelihoodSchlather(double *data, cov_model *cov, double *v) {
loglikelihoodMaxstable(data, cov, schlatherlogD, v);
}
//////////////////////////////////////////////////////////////////////
// Random Coin
//////////////////////////////////////////////////////////////////////
int check_randomcoin(cov_model *cov) {
cov_model
*pdf = cov->sub[COIN_COV],
*shape = cov->sub[COIN_SHAPE],
*next = shape != NULL ? shape : pdf,
*key = cov->key;
int err,
dim = cov->tsdim; // taken[MAX DIM],
mpp_param *gp = &(GLOBAL.mpp);
//extremes_param *ep = &(GLOBAL.extreme);
SERR("'random coin' method does not work for the current version");
INTERNAL;
ASSERT_ONE_SUBMODEL(cov);
ROLE_ASSERT(ROLE_POISSON_GAUSS || (cov->role==ROLE_GAUSS && cov->key!=NULL));
kdefault(cov, RANDOMCOIN_INTENSITY, gp->intensity[dim]);
kdefault(cov, RANDOMCOIN_METHOD, 0);
if ((err = checkkappas(cov, false)) != NOERROR) return err;
if (cov->tsdim != cov->xdimprev || cov->tsdim !=cov->xdimown) return ERRORDIM;
if (key != NULL) {
// note: key calls first function to simulate the points
// if this is true then AverageIntern/RandomCoinIntern calls Average/RandomCoin
cov_model *intern = cov;
while (intern != NULL && intern->nr!=AVERAGE_INTERN && isAnyDollar(intern))
intern = intern->key != NULL ? intern->key : intern->sub[0];
if (intern == NULL) {
BUG;
}
if ((err = cov->role == ROLE_BASE || cov->role == ROLE_GAUSS
// internal coin
? CHECK(key, dim, dim, ProcessType, XONLY, //err
CoordinateSystemOf(cov->isoown),
SUBMODEL_DEP,
ROLE_POISSON_GAUSS) // info wird von do_mpp verwendet
: cov->role == ROLE_POISSON_GAUSS
? CHECK(key, dim, dim, PointShapeType, XONLY, //err
CoordinateSystemOf(cov->isoown),
SUBMODEL_DEP, ROLE_POISSON_GAUSS)
: ERRORFAILED
) != NOERROR) {
return err;
}
} else { // key == NULL
if (next == pdf) {
if ((err = CHECK(next, dim, dim, PosDefType, XONLY, SYMMETRIC,
SCALAR, ROLE_POISSON_GAUSS)) != NOERROR) {
return err;
}
if ((pdf->pref[Average] == PREF_NONE) +
(pdf->pref[RandomCoin] == PREF_NONE) !=1){
//APMI(pdf);
//assert(false);
return ERRORPREFNONE;
}
} else { // shape != NULL
if ((err = CHECK(next, dim, dim, ShapeType, XONLY,
CoordinateSystemOf(cov->isoown),
SCALAR, ROLE_POISSON)) != NOERROR) { //POISS_GAUSS ??
return err;
}
// APMI(cov);
}
}
setbackward(cov, key != NULL ? key : next);
if ((err = kappaBoxCoxParam(cov, GAUSS_BOXCOX)) != NOERROR) return err;
if ((err = checkkappas(cov, false)) != NOERROR) return err;
return NOERROR;
}
void range_randomcoin(cov_model VARIABLE_IS_NOT_USED *cov, range_type *range) {
GAUSS_COMMON_RANGE;
range->min[RANDOMCOIN_INTENSITY] = RF_NEGINF;
range->max[RANDOMCOIN_INTENSITY] = RF_INF;
range->pmin[RANDOMCOIN_INTENSITY] = -10;
range->pmax[RANDOMCOIN_INTENSITY] = 10;
range->openmin[RANDOMCOIN_INTENSITY] = true;
range->openmax[RANDOMCOIN_INTENSITY] = true;
range->min[RANDOMCOIN_METHOD] = 0;
range->max[RANDOMCOIN_METHOD] = 1;
range->pmin[RANDOMCOIN_METHOD] = 0;
range->pmax[RANDOMCOIN_METHOD] = 1;
range->openmin[RANDOMCOIN_METHOD] = false;
range->openmax[RANDOMCOIN_METHOD] = false;
}
int struct_randomcoin(cov_model *cov, cov_model **newmodel){
cov_model *tmp_shape = NULL,
*pdf = cov->sub[COIN_COV],
*shape = cov->sub[COIN_SHAPE];
location_type *loc = Loc(cov);
int err,
dim = cov->tsdim; // taken[MAX DIM],
ROLE_ASSERT(ROLE_POISSON_GAUSS);
if (loc->Time || (loc->grid && loc->caniso != NULL)) {
TransformLoc(cov, true, GRIDEXPAND_AVOID, false);
SetLoc2NewLoc(pdf == NULL ? shape : pdf, PLoc(cov));
}
if (cov->key != NULL) COV_DELETE(&(cov->key));
ASSERT_NEWMODEL_NULL;
if (pdf != NULL) {
if ((err = CHECK(pdf, dim, dim, PosDefType, XONLY, SYMMETRIC,
SCALAR, ROLE_POISSON_GAUSS)) != NOERROR) {
//sicherstellen, dass pdf von der richtigen Form, insb. role_poisson_gauss
return err;
}
if (pdf->pref[Average] == PREF_NONE && pdf->pref[RandomCoin]==PREF_NONE) {
// if (pdf->nr > LASTDOLLAR) AERR(ERRORPREFNONE);
return ERRORPREFNONE;
}
if ((err = STRUCT(pdf, &tmp_shape)) != NOERROR)
goto ErrorHandling; //&(cov->key)
if (tmp_shape == NULL)
GERR("no structural information for random coins given");
tmp_shape->calling = cov;
if ((err = CHECK(tmp_shape, dim, dim, ShapeType, XONLY,
CoordinateSystemOf(cov->isoown),
SCALAR, ROLE_POISSON_GAUSS)) != NOERROR) {
goto ErrorHandling;
}
} else {
tmp_shape = shape;
// if ((err = covCpy(&(cov->key), shape)) > NOERROR) {
// return err;
// }
}
// if ((err = STRUCT(cov, NULL)) != NOERROR) return err; ?????
SERR("Sorry, 'random coin' does not work currently.");
assert(false); // todo
switch(P0INT(RANDOMCOIN_METHOD)) {
case 0: {
// Alternativ?!
// if ((err = addPointShape(&(cov->key), shape, NULL, tsdim, vdim)) != NOERROR)
if ((err = struct_ppp_pts(&(cov->key), tmp_shape, cov,
cov->tsdim, cov->vdim[0]))
!= NOERROR) goto ErrorHandling;
break;
}
case 1: {
if ((err = covCpy(&(cov->key), shape)) != NOERROR) goto ErrorHandling;
addModel(&(cov->key), RANDOMSIGN, cov);
addModel(&(cov->key), MCMC_PGS, cov);
cov_model *key = cov->key;
if ((err = covCpy(key->sub + PGS_LOC, shape)) != NOERROR)
goto ErrorHandling;
addModel(key, PGS_LOC, POW);
kdefault(key->sub[PGS_LOC], POW_ALPHA, GLOBAL.mpp.shape_power);
addModel(key, PGS_LOC, SCATTER);
PARAMALLOC(key, SCATTER_MAX, dim, 1);
for (int i=0; i<dim; i++) PARAM(key, SCATTER_MAX)[i] = NA_INTEGER;
PARAMALLOC(key, SCATTER_STEP, dim, 1);
for (int i=0; i<dim; i++) PARAM(key, SCATTER_STEP)[i] = RF_NA;
addModel(&(cov->key), MCMC, cov);
break;
}
default: BUG;
}
if ((err = CHECK(cov->key, dim, dim, PointShapeType, XONLY,
CoordinateSystemOf(cov->isoown),
cov->vdim[0], cov->role)) != NOERROR) goto ErrorHandling;
ErrorHandling:
if (pdf != NULL && tmp_shape != NULL) COV_DELETE(&tmp_shape);
return err;
}
int init_randomcoin(cov_model *cov, gen_storage *S) {
cov_model
*covshape = cov->sub[ cov->sub[COIN_SHAPE] != NULL ? COIN_SHAPE : COIN_COV],
*key = cov->key,
*sub = key == NULL ? covshape : key;
char name[] = "Poisson-Gauss";
location_type *loc = Loc(cov);
//mpp_storage *s = &(S->mpp);
int
err = NOERROR;
sprintf(ERROR_LOC, "%s process: ", name);
if (cov->role != ROLE_POISSON_GAUSS) {
// PMI(cov->calling, "init_randomcoin");
ILLEGAL_ROLE;
}
cov->method = covshape->pref[Average] == PREF_NONE ? RandomCoin : Average;
if (cov->method == Average && loc->caniso != NULL) {
bool diag, quasidiag, semiseparatelast, separatelast;
int idx[MAXMPPDIM];
assert(loc->timespacedim <= MAXMPPDIM);
analyse_matrix(loc->caniso, loc->cani_nrow, loc->cani_ncol, &diag,
&quasidiag, idx, &semiseparatelast, &separatelast);
if (!separatelast) SERR("not a model where time is separated");
}
if ((err = init_mpp(cov, S)) != NOERROR) {
//PMI(cov); AERR(err);
return err;
}
sub->Spgs->intensity =
sub->Spgs->totalmass * P0(RANDOMCOIN_INTENSITY);
sub->Spgs->log_density = log(P0(RANDOMCOIN_INTENSITY));
// PMI(cov, "randmcoin");
assert(sub->mpp.moments >= 2);
if (!R_FINITE(sub->Spgs->totalmass) || !R_FINITE(sub->mpp.mM[2]))
SERR("Moments of submodels not known");
sub->role = ROLE_POISSON_GAUSS;
return NOERROR;
}
void do_randomcoin(cov_model *cov, gen_storage *s) {
assert(cov->simu.active);
SAVE_GAUSS_TRAFO;
double *res = cov->rf;
dompp(cov, cov->Sgen==NULL ? s : cov->Sgen);// letzteres falls shape gegeben
BOXCOX_INVERSE;
}
