/* Authors Martin Schlather, schlather@math.uni-mannheim.de simulation of max-stable random fields Copyright (C) 2001 -- 2017 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 #include #include "questions.h" #include "Processes.h" #include "shape.h" #include "operator.h" #define POISSON_INTENSITY 0 #define RANDOMCOIN_INTENSITY (COMMON_GAUSS + 1) #define RANDOMCOIN_METHOD (COMMON_GAUSS + 2) /* not used double GetTotal(model* cov) { defn *C = DefList + COVNR; 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; isub[i]); } for (i=0; ikappasub[i] != NULL) v *= GetTotal(cov->kappasub[i]); } } return v; } */ void range_mpp(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(model *cov, gen_storage *S) { // cov ist process model *sub = cov->key; if (sub == NULL) sub = cov->sub[cov->sub[0] == NULL]; if (sub == NULL) SERR("substructure could be detected"); location_type *loc = Loc(cov); //window_info *w = &(S->window); int err; bool poisson = hasPoissonFrame(sub), maxstable = hasMaxStableFrame(sub); pgs_storage *pgs; if (!maxstable && !poisson) ILLEGAL_FRAME; if (!equalsnowPointShape(sub)) SERR1("%.50s is not shape/point function", NICK(sub)); if (loc->distances) RETURN_ERR(ERRORFAILED); if ((err = INIT(sub, maxstable ? 1 : hasPoissonFrame(sub) ? 0 : 2, S)) != NOERROR) RETURN_ERR(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) { SERR("integral of positive part of submodel unkown"); } 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 (hasPoissonGaussFrame(sub)) { if (ISNAN(sub->mpp.mM[2]) || !R_FINITE(sub->mpp.mM[2] || sub->mpp.mM[2] <=0.0)){ SERR("second moment unkown"); } } else assert(hasPoissonFrame(sub)); if ((err = ReturnOwnField(cov)) != NOERROR) { // must be later than INIT ! RETURN_ERR(err); } cov->initialised = cov->simu.active = err == NOERROR; RETURN_ERR(err); } #define SET_SUB { \ sub = cov->key; /* if changed, change do_max_pgs! */ \ if (sub == NULL) sub = cov->sub[cov->sub[0] == NULL]; \ if (sub == NULL) ERR("substructure could not 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(model *cov, gen_storage *s, double *simuxi) { // cov ist process assert(cov->simu.active); location_type *loc = Loc(cov); if (Getcaniso(cov) != NULL) BUG; model *sub = NULL, *key = cov->key; pgs_storage *pgs = NULL; long spatial = Getspatialpoints(cov); int err, *gridlen = NULL, *end = NULL, *start = NULL, *delta = NULL, *nx = NULL, dim = ANYOWNDIM, every = GLOBAL.general.every, // nthreshold = (every>0) ? every : MAXINT, 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, Minimum = GLOBAL.extreme.min_shape_gumbel, // -10 // TO DO; -1e150 *res = cov->rf; coord_type xgr = Getxgr(cov); assert(res != NULL); ext_bool fieldreturn, loggiven; bool maxstable = hasMaxStableFrame(key), simugrid = Getgrid(cov), partialfield = false, // AVERAGE braucht Sonderbehandlung: randomcoin = cov->method == RandomCoin //only for FRAME_ POISSON_ GAUSS ; long zaehler, nn, cumgridlen[MAXMPPDIM +1], Total_n, total = Gettotalpoints(cov), vdim = VDIM0, vdimtot = total * vdim, control = 0; assert(VDIM0 == VDIM1); if (vdim != 1) ERR("Poisson point process based methods only work in the univariate case"); SET_SUB; if (!equalsnowPointShape(sub)) BUG; bool compoundpoisson = hasPoissonFrame(sub), poissongauss = hasPoissonGaussFrame(sub); Total_n = -1; if (!maxstable) Total_n = rpois(pgs->intensity); // printf("pg = %d %10g %d\n", poissongauss, pgs->intensity, maxstable); // A PMI0(sub); 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 (int d=0; dmpp.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 = 1e15; //RF_INF; pgs->globalmin = Minimum; pgs->currentthreshold = loggiven ? pgs->globalmin - threshold : pgs->globalmin / threshold; // printf("pgs->current = %10g loggiven=%d %10g %10g\n", pgs->currentthreshold, loggiven, threshold, pgs->globalmin); // BUG; if (simuxi != NULL) { RFERROR("extremal t not programmed yet"); // to do : RPstudent, extremal t process } } else if (compoundpoisson){ assert(simuxi == NULL); for (int i=0; impp.mM[2]); pgs->currentthreshold = GLOBAL.mpp.about_zero; } for(nn=1; ; nn++) { //printf("%d %d tot=%d\n", (int) nn, (int) control, (int) Total_n); assert(nn < 1000); // if (n % 1000 == 0) // printf("n=%d tot=%d contr=%d: %10g < %10g; res[0]=%10g\n", (int) n, (int) Total_n, (int) control, res[control], threshold, res[0]); // assert(n <= 10); if (!maxstable && (nn > Total_n)) break; // for (d=0; dmin[d] = info->max[d] = 0.0; if (sub->randomkappa) { if ((err = REINIT(sub, sub->mpp.moments, s)) != NOERROR) BUG; DO(sub, s); SET_SUB; } else DO(sub, s); // printf("%10g %d\n", sub->rf[0], sub->fieldreturn); //APMI0(sub); fieldreturn = sub->fieldreturn; // MARCO: hattest Du dies reingenommen? // if (subf rame == BrMethodType) { // assert(sub->SBR != NULL); // n += sub->SBR->hatnumber; // } //printf("log_den %10g\n", pgs->log_density); logdens = pgs->log_density; if (!R_FINITE(logdens)) { BUG; // get_shape = DefList[sub->gatternr].cov; // get_logshape = DefList[sub->gatternr].log; } if (compoundpoisson) { summand = 0.0; } else if (poissongauss) { summand = -0.5 * (logdens + logM2); } else { // max-stable field assert(hasMaxStableFrame(sub)); if (nn > 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(poisson); threshold = logthreshold = (double) (gumbel+LOG(sub->mpp.maxheights[0])); if (!loggiven) threshold = EXP(threshold); // Frechet summand = gumbel - logdens; //@MARCO: summand ist eine rein // technische Variable, die die Anzahl der Berechnungen ein // winziges bisschen reduziert im Falle des Zhou-Ansatzes. //printf("logdens = %10g\n", logdens);BUG; //printf("for thres=%10e %d res=%10e zhou=%10g(%10g) logden=%10g gumbel=%10g loggiven=%d summand=%10g\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=threshold)) { // printf("control: %d %10g %10g\n", (int) control, res[control], threshold); //assert(false); control++; } if (control >= total) { break; } // printf("n=%d every =%d %d %d\n", (int) n, GLOBAL.extreme.check_every, PL, PL_STRUCTURE); if (nn % GLOBAL.extreme.check_every == 0) { double min = res[control]; for (int i=control + 1; i= PL_STRUCTURE) { PRINTF("control: %ld %10g %10g glob=%10g n=%ld every=%d %10g log.max=%10g\n", // %ld OK control, res[control], threshold, min, nn, GLOBAL.extreme.check_every, sub->mpp.maxheights[0], LOG(sub->mpp.maxheights[0])); } pgs->globalmin = min < Minimum ? Minimum : min; } pgs->currentthreshold = loggiven ? pgs->globalmin - gumbel : pgs->globalmin / frechet; // printf("xxx = %10e %10e %10e Min=%10e \n",pgs->globalmin, gumbel, pgs->currentthreshold, Minimum); } // 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); // 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; int d; for (d=0; dsupportmin[d] - 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] - xgr[d][XSTART]) / step )); end[d] = dummy < 0 ? 0 : dummy > gridlen[d] ? gridlen[d] : (int) dummy; partialfield |= end[d] < gridlen[d]; // printf("window [%10g %10g] [%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=%10g %10g, loc.start=%10g %10g\n", d, start[d], end[d], pgs->supportmin[d], pgs->supportmax[d], 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] = xgr[d][XSTART] + (double) start[d] * inc[d]; if (false || zaehler < 0) { PRINTF("d=%d start=%d, end=%d xstart=%10g %10g pgs=[%4.3f, %4.3f] xgr=%10g %10g %10g inc=%3.2f\n", d, start[d], end[d], xstart[d], pgs->xstart[d], pgs->supportmin[d], pgs->supportmax[d], xgr[d][XSTART], xgr[d][XSTEP], xgr[d][XLENGTH], inc[d]); // assert(false); } } if (d < dim) { //printf("cont: d=%d\n", d); continue; } } #define STANDARDINKREMENT_ZAEHLER \ int 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("%10g : %10g fctr=%10g %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]); #define RFMAX(OP) { \ rf = sub->rf; \ for (int i=0; irf; \ for (int i=0; i=0 && zaehler=total) { 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=%10g %10g delta=%10e %d %10e\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)); */ \ res[zaehler] = res[zaehler] < *value ? *value : res[zaehler]; \ } #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=0 && zaehler < total); \ int jj, idx, index=0; \ for (jj=0; jj= pgs->own_grid_len[jj]) break; \ index += pgs->own_grid_cumsum[jj] * idx; \ } \ if (jj >= dim) { OP; } else if (false) printf("%d %10g %10g %10g len=%10g %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; res[zaehler] = res[zaehler] < *value ? *value : res[zaehler];) #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){/* the realized 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) { \ double zw,inct, segt, ct, st, A, cit, sit; \ inct = sub->q[AVERAGE_YFREQ] * xgr[dim][XSTEP]; \ cit = COS(inct); \ sit = SIN(inct); \ G; \ segt = OP1; \ A = OP0; \ ct = A * COS(segt); \ st = A * SIN(segt); \ for (long zt = zaehler; zt < total; 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 (; zaehlerloggiven) { AVEAPPLY(LOGGET, Sign[0] *EXP(value[0] + summand), Sign[1] *EXP(value[1])); } else AVEAPPLY(GET, value[0] * factor, value[1]); } else { // printf("here\n"); assert(hasMaxStableFrame(sub) || hasPoissonFrame(sub)); // reicht das? if (fieldreturn == Huetchenownsize) { // atomdependent rfbased/not; ALL_APPLY(OWNGRIDAPPLY, NULL, OG_MAXAPPLY, NULL, OG_SUMAPPLY, NULL); } else if (fieldreturn == wahr) { // 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 == falsch) { // not field return // printf("kein feld\n"); ALL_APPLY(APPLY, LOGGETPOS, MAXAPPLY, LOGGET, SUMAPPLY, GET); } else BUG; // neither Huetchenownsize, nor true nor false } if (nn > nthreshold) { if (maxstable) { LPRINT("%ld%.50s pos.: value=%1.3f threshold=%1.3f gu=%1.3f %1.3f %d (%ld %d)\n", // %ld ok control, TH(control), (double) res[control], (double) threshold, gumbel, logdens, loggiven, nn, nthreshold); //, deltathresh); nthreshold += deltathresh; } else { PRINTF("n=%d Total=%d\n", (int) nn, (int)(Total_n)); } } R_CheckUserInterrupt(); }// for(n,... -- while control < total // formerly (before 19 Aug 2017: contents of // void finalmaxstable(model *cov, double *res, int n) // had been here --- much better to have it at the very end -- // much less addional evaluations necessary to get c_zhou! // for (k=0; k= PL_DETAILSUSER) { PRINTF("number of shape functions used: %ld\n", nn); } return; } // end dompp void dompp(model *cov, gen_storage *s) { dompp(cov, s, NULL); } // 4805920678128 121 void finalmaxstable(model *cov, double *Res, int n, gen_storage *s) { // cov is process model *key = cov->key, *sub = cov->key;/* if changed, change do_max_pgs! */ if (sub == NULL) sub = cov->sub[cov->sub[0] == NULL]; if (sub == NULL) ERR("substructure could be detected"); long total = Gettotalpoints(cov), vdim = VDIM0, vdimtot = total * vdim; double meansq, var, mean = RF_NA, n_min = RF_INF, eps = GLOBAL.extreme.eps_zhou, *endres = Res + n * vdimtot; ext_bool loggiven = key == NULL ? cov->sub[0]->loggiven : key->loggiven; pgs_storage *pgs = sub->Spgs; assert(pgs != NULL); 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) { int err; if ((err = REINIT(sub, sub->mpp.moments, s)) != NOERROR) BUG; DO(sub, s); } while (true) { // n_min depends on the estimation precision 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 (PL > 5) { PRINTF(" n=%ld, min=%10g %10g mean=%10g zhou=%10g eps=%10g\n", pgs->n_zhou_c, n_min, (double) pgs->sum_zhou_c, mean, pgs->zhou_c, eps); } if (n_min >= GLOBAL.extreme.max_n_zhou || pgs->n_zhou_c >= n_min) break; int endfor=(n_min - pgs->n_zhou_c); endfor = endfor < 10 ? 10 : endfor > 50 ? 50 : endfor; for (int k=0; kmpp.moments, s)) != NOERROR) BUG; DO(sub, s); // SET_SUB missing ?? } } // } else { mean = pgs->zhou_c; } // mean /= M1; //printf(" n=%ld, min=%10g, mean=%10g zhou=%10g eps=%10g %10g max=%10g\n", pgs->n_zhou_c, n_min, (double) pgs->sum_zhou_c, mean, pgs->zhou_c, eps, sub->mpp.maxheights[0]); // printf("loggiven = %d %10g %10g\n", loggiven, mean, pgs->zhou_c); double xi = P0(GEV_XI); if (loggiven && !pgs->logmean) mean =LOG(mean); // sieht noch umstaendlich aus, aber jede variable kann ein anderes // xi irgendwann haben. for (double *res = Res; resalpha != 1.0) { // default = 1 in NULL.cc; s. opitz if (loggiven) for (int k=0; kalpha; else for (int k=0; kalpha); } // printf("xi=%10g %d logmean=%d %10g\n", xi, loggiven, pgs->logmean, mean); if (loggiven) { // 3.6.13 sub->loggiven geaendert zu loggiven for (int k=0; kkappasub[GEV_S] != NULL) { ERR1("'%.50s' 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/xi=%10g\n", ss); if (FABS(ss - 1.0) > 1e-15) for (int k=0; kkappasub[GEV_MU] != NULL) { ERR1("'%.50s' cannot be chosen as a function yet.", KNAME(GEV_MU)); } else { double mu = P0(GEV_MU); // printf("mu=%10g xi=%10g %10g\n", mu, xi, s); // printf("mu=%10g\n", mu); if (FABS(mu) > 1e-15) for (int k=0; kcalling, // cov is poisson process *next = *Cov; ASSERT_ONESYSTEM; int err, dim = XDIM(PREVSYSOF(next), 0), vdim = next->vdim[0]; Types frame = next->frame; assert(VDIM0 == VDIM1); 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(OWNISO(0)), \ vdim, frame)) != NOERROR) RETURN_ERR(err) checkstandardnext; if (!CallingSet(*Cov)) BUG; if (hasPoissonFrame(next)) { addModel(key, PGS_LOC, UNIF); assert(key->nsub == 2); // ??? PARAMALLOC(key->sub[PGS_LOC], UNIF_MIN, dim, 1); PARAMALLOC(key->sub[PGS_LOC], UNIF_MAX, dim, 1); } else { // TREE(cov); // printf("structh %ld %.50s\n", key->sub + PGS_LOC, TYPE_NAMES[cov->frame]); if ((err = STRUCT(key, key->sub + PGS_LOC)) != NOERROR) RETURN_ERR(err); SET_CALLING(key->sub[PGS_LOC], key); } if (!CallingSet(*Cov)) BUG; checkstandardnext; RETURN_NOERROR; } int check_poisson(model *cov) { model *next=cov->sub[0], *key = cov->key, *sub = key == NULL ? next : key; int err, dim = ANYDIM; // taken[MAX DIM], mpp_param *gp = &(GLOBAL.mpp); Types type = sub == key ? PointShapeType : ShapeType; //print("eee\n"); kdefault(cov, POISSON_INTENSITY, gp->intensity[dim]); if ((err = checkkappas(cov, true)) != NOERROR) RETURN_ERR(err); if ((err = CHECK(sub, dim, dim, type, XONLY, CoordinateSystemOf(OWNISO(0)), SUBMODEL_DEP, PoissonType)) != NOERROR) RETURN_ERR(err); setbackward(cov, sub); RETURN_NOERROR; } void range_poisson(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(model *cov, model **newmodel){ model *next=cov->sub[0]; location_type *loc = Loc(cov); // if (newmodel != NULL) crash(); ASSERT_NEWMODEL_NULL; assert(isnowProcess(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 (!equalsnowPointShape(next)) { int err; if ((err = covcpy(&(cov->key), next)) != NOERROR) RETURN_ERR(err); if ((err = addStandardPoisson(&(cov->key))) != NOERROR) RETURN_ERR(err); } RETURN_NOERROR; } int init_poisson(model *cov, gen_storage *S) { // location_type *loc = Loc(cov); // mpp_storage *s = &(S->mpp); model *key=cov->key; int err; if ((err = init_mpp(cov, S)) != NOERROR) { RETURN_ERR(err); } if (!equalsnowPointShape(key)) SERR("no definition of a shape function found"); key->Spgs->intensity = key->Spgs->totalmass * P0(POISSON_INTENSITY); cov->simu.active = cov->initialised = err == NOERROR; RETURN_ERR(err); } //////////////////////////////////////////////////////////////////// // Schlather void extremalgaussian(double *x, model *cov, double *v) { // schlather process model *next = cov->sub[0]; COV(x, next, v); if (hasGaussMethodFrame(next)) *v = 1.0 - SQRT(0.5 * (1.0 - *v)); } //#define MPP_SHAPE 0 //#define MPP_TCF 1 int SetGEVetc(model *cov) { // cov is process int err; extremes_param *ep = &(GLOBAL.extreme); if (cov->sub[MPP_TCF] != NULL && cov->sub[MPP_SHAPE]!=NULL) SERR2("either '%.50s' or '%.50s' must be given", SNAME(MPP_TCF), SNAME(MPP_SHAPE)); if ((err = checkkappas(cov, false)) != NOERROR) RETURN_ERR(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 %10g %10g %10g\n", P0(GEV_XI), P0(GEV_S), P0(GEV_MU)); assert(false); if ((err = checkkappas(cov, true)) != NOERROR) RETURN_ERR(err); RETURN_NOERROR; } int check_schlather(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."); model *key = cov->key, *next = cov->sub[cov->sub[MPP_TCF] != NULL ? MPP_TCF : MPP_SHAPE]; int err, dim = ANYDIM; // taken[MAX DIM], // mpp_param *gp = &(GLOBAL.mpp); double v; bool schlather = DefList[COVNR].Init == init_mpp; // else opitz // printf("A \n"); ASSERT_ONE_SUBMODEL(cov); /// printf("B A \n"); if ((err = SetGEVetc(cov)) != NOERROR) RETURN_ERR(err); // printf("CA \n"); model *sub = cov->key==NULL ? next : key; // printf("GA \n"); if (key == NULL) { // printf("key=NULL\n"); if (equalsBernoulliProcess(sub) && !schlather) SERR1("'%.50s' does not work with Bernoulli processes", NAME(cov)); Types frame = isProcess(sub) || isPointShape(sub) ? SchlatherType : EvaluationType; //isGaussMethod(sub) ? GaussMethodType // : isPointShape(sub) && schlather ? SchlatherType //: equalsBernoulliProcess(sub) && schlather ? NormedProcessType // egal //: isProcess(sub) ? SchlatherType // : EvaluationType; // printf(">>> %.50s %d %.50s\n", TYPE_NAMES[frame], isProcess(sub), NAME(sub)); // if (isProcess(next)) err = CHECK(next, dim, dim, ProcessType, XONLY, CoordinateSystemOf(OWNISO(0)), SCALAR, frame); else err = CHECK(next, dim,dim, PosDefType, XONLY, IsotropicOf(OWNISO(0)), SCALAR, frame); if (err != NOERROR) RETURN_ERR(err); if (sub->vdim[0] != 1) SERR("only univariate processes are allowed"); assert(VDIM0 == VDIM1); setbackward(cov, sub); if (hasEvaluationFrame(next)) { if (next->pref[Nothing] == PREF_NONE) RETURN_ERR(ERRORPREFNONE); COV(ZERO(next), next, &v); if (v != 1.0) SERR2("a correlation function is required as submodel, but '%.50s' has variance %10g.", NICK(next), v); } } else { // key != NULL //printf("key!=NULL\n"); if ( (err = CHECK(key, dim, dim, PointShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), SUBMODEL_DEP, SchlatherType)) != NOERROR) { RETURN_ERR(err); } setbackward(cov, sub); } //print("end check schlather\n"); RETURN_NOERROR; } int struct_schlather(model *cov, model **newmodel){ model *sub = cov->sub[cov->sub[MPP_TCF] != NULL ? MPP_TCF : MPP_SHAPE]; int err, ErrNoInit; bool schlather = DefList[COVNR].Init == init_mpp; // else opitz // assert(isnowSchlather(cov)); 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(err); SET_CALLING(cov->key, cov); } else { if ((err = covcpy(&(cov->key), sub)) != NOERROR) RETURN_ERR(err); } assert(cov->key->calling == cov); if (MODELNR(cov->key) != GAUSSPROC && !equalsBernoulliProcess(cov->key) && MODELNR(cov->key) != BRNORMED ) { if (isnowVariogram(cov->key)) addModel(&(cov->key), GAUSSPROC); else { if (isGaussMethod(cov->key)) { SERR("invalid model specification"); } else SERR2("'%.50s' currently only allowed for gaussian processes %.50s", NICK(cov), schlather ? "and binary gaussian processes" : ""); } } assert(cov->key->calling == cov); Types frame = SchlatherType; // MODELNR(cov->key) == GAUSSPROC ? GaussMethodType //Marco, 29.5.13 // : equalsBernoulliProcess(cov->key) ? NormedProcessType // egal // : SchlatherType; // if ((err = CHECK(cov->key, cov->tsdim, cov->xdimown, ProcessType, // OWNDOM(0), // OWNISO(0), cov->vdim, frame)) != NOERROR) RETURN_ERR(err); if ((err = CHECK_PASSTF(cov->key, ProcessType, VDIM0, frame)) != NOERROR) RETURN_ERR(err); if ((ErrNoInit = STRUCT(cov->key, NULL)) > NOERROR)//err return ErrNoInit; addModel(&(cov->key), STATIONARY_SHAPE); // if ((err = CHECK(cov->key, cov->tsdim, cov->xdimown, PointShapeType, // OWNDOM(0), // OWNISO(0), cov->vdim, SchlatherType)) != NOERROR) //RETURN_ERR(err); if ((err = CHECK_PASSTF(cov->key, PointShapeType, VDIM0, SchlatherType)) != NOERROR) RETURN_ERR(err); return ErrNoInit; } #define OPITZ_ALPHA (LAST_MAXSTABLE + 1) int init_opitzprocess(model *cov, gen_storage *S) { int err; if ( (err = init_mpp(cov, S)) != NOERROR) RETURN_ERR(err); double alpha = P0(OPITZ_ALPHA); 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 * alpha - 0.5) * gammafn(0.5 * alpha + 0.5); pgs->zhou_c = 1.0 / key->mpp.mMplus[1]; pgs->alpha = alpha; cov->simu.active = cov->initialised = true; RETURN_NOERROR; } void range_opitz(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(model *to, 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 0) { for (i=0; isub[i] != NULL) param_set_identical(to->sub[i], from->sub[i], depth - 1); } } } int FillInPts(model *cov, // struct of shape model *shape // shape itself ) { // 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 // // smith, randomcoin assert(cov != NULL); ASSERT_ONESYSTEM; // printf("fx \n"); int err, i, dim = XDIM(SYSOF(shape), 0), nr = COVNR; if (cov->sub[PGS_LOC] != NULL) RETURN_NOERROR; if ((err = covcpy(cov->sub + PGS_FCT, shape)) != NOERROR) RETURN_ERR(err); if (nr == ZHOU) { //printf("here %ld %d %d %d\n", cov->sub[PGS_LOC], ScaleOnly(shape), // !shape->randomkappa, shape->sub[0]->randomkappa); cov->nsub = 2; if (cov->sub[PGS_LOC] != NULL) BUG; bool randomscale = ScaleOnly(shape) && shape->randomkappa && !shape->sub[0]->randomkappa; if ((err = covcpyWithoutRandomParam(cov->sub + PGS_LOC, randomscale ? shape->sub[0] : shape)) != NOERROR) RETURN_ERR(err); if (hasPoissonGaussFrame(shape) //|| hasMaxStableFrame(shape) ) { assert(dim > 0 && dim <= MAXMPPDIM); addModel(cov, PGS_LOC, POW); assert(cov->nsub == 2); kdefault(cov->sub[PGS_LOC], POW_ALPHA, GLOBAL.mpp.shape_power); addModel(cov, PGS_LOC, SCATTER); model *scatter = cov->sub[PGS_LOC]; // PMI0(scatter); printf("sm=%d %d %d\n", SCATTER_MAX, dim,SCATTER_STEP); PARAMALLOC(scatter, SCATTER_MAX, dim, 1); for (i=0; irandomkappa) { addSetDistr(cov->sub + PGS_LOC, cov->sub[PGS_FCT], param_set_identical, true, MAXINT); } addModel(cov, PGS_LOC, RECTANGULAR); if (randomscale) { addModel(cov, PGS_LOC, LOC); addSetDistr(cov->sub + PGS_LOC, shape, ScaleDollarToLoc, true, 0); } } else if (nr == STANDARD_SHAPE) { assert(cov != NULL && shape != NULL); if ((err = STRUCT(shape, cov->sub + PGS_LOC)) != NOERROR) RETURN_ERR(err); SET_CALLING(cov->sub[PGS_LOC], cov); if (cov->sub[PGS_LOC] == NULL) SERR("simple enlarging of the simulation window does not work"); } else if (nr == BALLANI) { RETURN_ERR(ERRORNOTPROGRAMMEDYET); } else BUG; RETURN_NOERROR; } int addPGS(model **Key, // struct of shape model *shape,// shape itself model *pts, int dim, int vdim, Types frame) { // SEE ALSO addPGSLocal in RMS.cc !!!!!! /// random coin und smith // versucht die automatische Anpassung einer PointShapeType-Funktion; // derzeit wird // * PGS und // * STANDARD_SHAPE (weiss nicht mehr wieso -> coins?) // versucht; inklusive Init // // ruft t.w. FillInPts auf #define specific (nPOISSON_SCATTER - 1) bool maxstable = hasMaxStableFrame(shape); int i, err = NOERROR, method = GLOBAL.extreme.scatter_method, pgs[specific] = {maxstable ? ZHOU : BALLANI, STANDARD_SHAPE}; #define infoN 200 char info[nPOISSON_SCATTER - 1][LENERRMSG]; assert(shape != NULL); assert(*Key == NULL); // to do: strokorbball: raeumlich waere Standard/Uniform besser; // praeferenzen programmieren? for (i=0; i 0) { errorMSG(err, info[i-1]); // XERR(err); // eigentlich muss das hier weg } // if (i > 0) XERR(err); assert(i ==0); if (*Key != NULL) COV_DELETE(Key); assert(shape->calling != NULL); addModel(Key, pgs[i], shape->calling); if (pts == NULL) { if ((err = FillInPts(*Key, shape)) != NOERROR) continue; } else { if ((err = covcpy((*Key)->sub + PGS_FCT, shape)) != NOERROR || (err = covcpy((*Key)->sub + PGS_LOC, pts)) != NOERROR ) { model *cov = *Key; RETURN_ERR(err); } // wechselseitiger kreuzverweis Ssetcpy((*Key)->sub[PGS_FCT], (*Key)->sub[PGS_LOC], shape, pts); Ssetcpy((*Key)->sub[PGS_LOC], (*Key)->sub[PGS_FCT], pts, shape); } model *cov = *Key; SET_CALLING(cov, shape->calling); SET_CALLING(cov->sub[PGS_FCT], cov); SET_CALLING(cov->sub[PGS_LOC], cov); assert(cov->sub[PGS_LOC] != NULL && cov->sub[PGS_FCT] != NULL); cov->nsub = 2; // printf(">>>>>>> KEY() start %d %.50s \n", i, NICK((*Key))); if ((err = CHECK(*Key, dim, dim, PointShapeType, XONLY, CoordinateSystemOf(ISO(SYSOF(shape), 0)), vdim, frame)) != NOERROR) { // printf(">>>>>>> KEY() break %d %.50s %d\n", i, Nick(*Key), dim); // XERR(err); continue; } NEW_COV_STORAGE(cov, gen); if ((err = INIT(cov, 1, cov->Sgen)) == NOERROR) break; } // for i_pgs model *cov = *Key; if (err != NOERROR) { SERR("error occured when creating the point-shape fctn"); // errorstring_type xx = "error occured when creating the point-shape fctn"; // errorMSG(err, xx, cov->base, info[i-1]); // SERR2("no method found to simulate the given model:\n\tpgs: %.50s\n\tstandard: %.50s)", info[0], info[1]); } RETURN_NOERROR; } int struct_ppp_pts(model **Key, // = *cov->key model *shape, model *calling, // = process int timespacedim, int vdim, Types frame) { // MAIN FUNCTION FOR CREATING THE POINT-SHAP_OBJECT FOR // SMITH AND RANDOM COIN // Ruft S TRUCT(shape, Key) auf und verarbeitet weiter, je nachdem // was nun in Key steht: pgs, shape, random, nichts. model *dummy = NULL; int err = NOERROR, logdim = LOGDIM(PREVSYSOF(shape), 0), xdim = XDIM(PREVSYSOF(shape), 0); #ifdef RANDOMFIELDS_DEBUGGING model *cov = shape; #endif err = STRUCT(shape, Key); // get the random location corresponding to shape; // either zhou paper or ballani if (err == NOERROR && *Key != NULL) { // indeed there is a corresponding random location function // could be both smith model and random coin SET_CALLING(*Key, calling); Types type = BadType; if (isPointShape(*Key)) type = PointShapeType; else if ((err = CHECK_R(*Key, logdim)) == NOERROR) type = RandomType; if (!equalsRandom(type)) { if ((err = CHECK(*Key, logdim, xdim, type, DOM(PREVSYSOF(shape), 0), ISO(PREVSYSOF(shape), 0), shape->vdim, frame)) != NOERROR) goto ErrorHandling; type = type == BadType ? ShapeType : type; } if (type == RandomType) { // random locations given for the shape fct; // so, it must be of pgs type (or standard): // could be both smith model and random coin #ifdef RANDOMFIELDS_DEBUGGING model *cov = *Key; #endif dummy = *Key; *Key = NULL; if ((err = addPGS(Key, shape, dummy, timespacedim, vdim, frame)) != NOERROR) goto ErrorHandling; if (*Key == NULL) BUG; SET_CALLING(*Key, calling); } else { model *cov = *Key; ASSERT_ONESYSTEM; if (type==PointShapeType) { // happens with struct_strokorbPoly if ((err = FillInPts(*Key, shape)) != NOERROR) goto ErrorHandling; } else { // (type == ShapeType) assert(CALLINGNR != SMITHPROC); // i.e. random coin // i.e. it delivers the (non-random) coin for a given // covariance function dummy = *Key; *Key = NULL; // suche nach geeigneten locationen if ((err = addPGS(Key, dummy, NULL, timespacedim, vdim, frame)) != NOERROR) goto ErrorHandling; //printf("e ddd\n"); } } // ! randomtype } else { // S-TRUCT does not return anything int err2; //assert(false); // XERR(err); if (*Key != NULL) { SET_CALLING(*Key, calling);// make sure that the locations are not deleted COV_DELETE(Key); } if ((err2 = addPGS(Key, shape, NULL, timespacedim, vdim, frame))!=NOERROR){ if (err == NOERROR) err = err2; goto ErrorHandling; } err = NOERROR; } ErrorHandling: if (dummy != NULL) COV_DELETE(&dummy); model *cov = *Key; RETURN_ERR(err); } int check_smith(model *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, dim = ANYDIM; // taken[MAX DIM], ASSERT_ONE_SUBMODEL(cov); if ((err = SetGEVetc(cov)) != NOERROR) RETURN_ERR(err); if (key != NULL) { if ((err = CHECK(key, dim, dim, PointShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), SUBMODEL_DEP, SmithType)) != NOERROR) RETURN_ERR(err); } else { // key == NULL if (next == TCF) { if ((err = CHECK(next, dim, dim, TcfType, XONLY, ISOTROPIC, SCALAR, SmithType)) != NOERROR) { RETURN_ERR(err); } if ((dim == 1 && next->rese_derivs < 1) || (dim >= 2 && dim <= 3 && next->rese_derivs < 2) || dim > 3) { //printf("dim = %d\n", dim); SERR("submodel does not have enough derivatives (programmed)."); } } else { // shape -- geaendert Types frame = SmithType; // isPointShape(sub) || is Shape(sub) ? SmithType //: isGaussMethod(sub) ? GaussMethodType //: equalsBernoulliProcess(sub) ? NormedProcessType // egal //: BadType; if ((err = CHECK(next, dim, dim, ShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), SCALAR, frame)) != NOERROR) { RETURN_ERR(err); } if (next->full_derivs < 0) SERR1("'%.50s' requires an explicit submodel.", NICK(cov)); } } setbackward(cov, next); RETURN_NOERROR; } int struct_smith(model *cov, model **newmodel){ 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, logdim = LOGDIM(PREVSYSOF(sub), 0), xdim = XDIM(PREVSYSOF(sub), 0); assert(hasSmithFrame(sub)); 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, logdim, xdim, ShapeType, DOM(PREVSYSOF(tcf), 0), ISO(PREVSYSOF(tcf), 0), tcf->vdim, SmithType)) != NOERROR) goto ErrorHandling; tmp_shape = tcf_shape; // SET_CALLING_NULL tcf_shape->calling ; ?? } else { tmp_shape = shape; } // if ((err = struct_ppp_pts(&(cov->key), tmp_shape, cov, ANYDIM, VDIM0, SmithType)) != NOERROR) goto ErrorHandling; err = NOERROR; ErrorHandling: if (tcf_shape != NULL && tmp_shape != NULL) COV_DELETE(&tmp_shape); RETURN_ERR(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, model *cov, // = process logDfct logD, double *v) { // DARF JA NICHT DURCH MEHRERE CORES AUFGERUFEN WERDEN!! model *sub = cov; while (isnowProcess(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); } ASSERT_ONESYSTEM; location_type *loc = Loc(cov); int dim = OWNTOTALXDIM; int len = loc->totalpoints; if (data != NULL) { double xi = P0(GEV_XI), mu = P0(GEV_MU), s = P0(GEV_S); if (xi == 0) { for (int i=0; iq[i] =EXP((data[i] - mu) / s); } else { double xiinv = 1.0 / xi, xis = xi / s; for (int i=0; iq[i] =POW(1.0 + xis * (data[i] - mu), xiinv); } } switch(GLOBAL.fit.likelihood) { case LL_AUTO: case LL_COMPOSITE: { double z[2], gamma, gamma0, x[MAXMPPDIM + 1], y[MAXMPPDIM + 1], *xx = loc->x; COV(ZERO(sub), sub, &gamma0); for (int i=0; iq[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, model *cov, double *v) { loglikelihoodMaxstable(data, cov, HueslerReisslogD, v); } void loglikelihoodSchlather(double *data, model *cov, double *v) { loglikelihoodMaxstable(data, cov, schlatherlogD, v); } ////////////////////////////////////////////////////////////////////// // Random Coin ////////////////////////////////////////////////////////////////////// int check_randomcoin(model *cov) { model *pdf = cov->sub[COIN_COV], *shape = cov->sub[COIN_SHAPE], *next = shape != NULL ? shape : pdf, *key = cov->key; int err, dim = ANYDIM; // 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); if ((!hasPoissonGaussFrame(next) && !hasGaussMethodFrame(next) && !hasProcessFrame(next)) || cov->key!=NULL) ILLEGAL_FRAME; kdefault(cov, RANDOMCOIN_INTENSITY, gp->intensity[dim]); kdefault(cov, RANDOMCOIN_METHOD, 0); if ((err = checkkappas(cov, false)) != NOERROR) RETURN_ERR(err); if (key != NULL) { // note: key calls first function to simulate the points // if this is true then AverageIntern/RandomCoinIntern calls Average/RandomCoin model *intern = cov; while (intern != NULL && MODELNR(intern) != AVERAGE_INTERN && isAnyDollar(intern)) intern = intern->key != NULL ? intern->key : intern->sub[0]; if (intern == NULL) { BUG; } if (hasGaussMethodFrame(next)) err = CHECK(key, dim, dim, ProcessType, XONLY, CoordinateSystemOf(OWNISO(0)), SUBMODEL_DEP, PoissonGaussType); else if (hasPoissonGaussFrame(next)) err = CHECK(key, dim, dim, PointShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), SUBMODEL_DEP, PoissonGaussType); else err = ERRORFAILED; if (err != NOERROR) RETURN_ERR(err); } else { // key == NULL if (next == pdf) { // s ymmetric: insbesondere also univariate if ((err = CHECK(next, dim, dim, PosDefType, XONLY, SYMMETRIC, SCALAR, PoissonGaussType)) != NOERROR) { RETURN_ERR(err); } if ((pdf->pref[Average] == PREF_NONE) + (pdf->pref[RandomCoin] == PREF_NONE) !=1){ //assert(false); RETURN_ERR(ERRORPREFNONE); } } else { // shape != NULL if ((err = CHECK(next, dim, dim, ShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), SCALAR, PoissonType)) != NOERROR) { //POISS_GAUSS ?? RETURN_ERR(err); } } } setbackward(cov, key != NULL ? key : next); if ((err = kappaBoxCoxParam(cov, GAUSS_BOXCOX)) != NOERROR) RETURN_ERR(err); if ((err = checkkappas(cov, false)) != NOERROR) RETURN_ERR(err); RETURN_NOERROR; } void range_randomcoin(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(model *cov, model **newmodel){ model *tmp_shape = NULL, *pdf = cov->sub[COIN_COV], *shape = cov->sub[COIN_SHAPE]; location_type *loc = Loc(cov); int err, dim = ANYDIM; // taken[MAX DIM], 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) {// s ymmetric: insbesondere also univariate if ((err = CHECK(pdf, dim, dim, PosDefType, XONLY, SYMMETRIC, SCALAR, PoissonGaussType)) != NOERROR) { //sicherstellen, dass pdf von der richtigen Form, insb. frame_poisson_gauss RETURN_ERR(err); } if (pdf->pref[Average] == PREF_NONE && pdf->pref[RandomCoin]==PREF_NONE) { // if (pdf->nr > LASTDOLLAR) AERR(ERRORPREFNONE); RETURN_ERR(ERRORPREFNONE); } if ((err = STRUCT(pdf, &tmp_shape)) != NOERROR) goto ErrorHandling; //&(cov->key) if (tmp_shape == NULL) GERR("no structural information for random coins given"); SET_CALLING(tmp_shape, cov); if ((err = CHECK(tmp_shape, dim, dim, ShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), SCALAR, PoissonGaussType)) != NOERROR) { goto ErrorHandling; } } else { tmp_shape = shape; // if ((err = covcpy(&(cov->key), shape)) > NOERROR) { // RETURN_ERR(err); // } } // if ((err = STRUCT(cov, NULL)) != NOERROR) RETURN_ERR(err); ????? SERR("Sorry, 'random coin' does not work currently."); assert(false); // todo switch(P0INT(RANDOMCOIN_METHOD)) { case 0: { // Alternativ?! // if ((err = addPGS(&(cov->key), shape, NULL, timespacedim, vdim)) != NOERROR) if ((err = struct_ppp_pts(&(cov->key), tmp_shape, cov, dim, VDIM0, PoissonGaussType)) != 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); 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); model *scatter = key->sub[PGS_LOC]; PARAMALLOC(scatter, SCATTER_MAX, dim, 1); for (int i=0; ikey), MCMC, cov); break; } default: BUG; } if ((err = CHECK(cov->key, dim, dim, PointShapeType, XONLY, CoordinateSystemOf(OWNISO(0)), VDIM0, PoissonGaussType)) != NOERROR) goto ErrorHandling; ErrorHandling: if (pdf != NULL && tmp_shape != NULL) COV_DELETE(&tmp_shape); RETURN_ERR(err); } int init_randomcoin(model *cov, gen_storage *S) { 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; KEY_type *KT = cov->base; SPRINTF(KT->error_loc, "%.50s process", name); assert(hasPoissonGaussFrame(sub)); 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) { RETURN_ERR(err); } sub->Spgs->intensity = sub->Spgs->totalmass * P0(RANDOMCOIN_INTENSITY); sub->Spgs->log_density =LOG(P0(RANDOMCOIN_INTENSITY)); assert(sub->mpp.moments >= 2); if (!R_FINITE(sub->Spgs->totalmass) || !R_FINITE(sub->mpp.mM[2])) SERR("Moments of submodels not known"); RETURN_NOERROR; } void do_randomcoin(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; }