/* Authors Martin Schlather, schlather@math.uni-mannheim.de Definition of correlation functions that are gaussian scale mixtures Copyright (C) 2017 -- 2018 Olga Moreva (bivariate models) & 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 PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include "def.h" #include #include #include #include #include #include "extern.h" #include "questions.h" #include "primitive.h" #include "operator.h" #include "AutoRandomFields.h" #include "shape.h" #include "Processes.h" //#include "xport_import.h" //#define LOG2 M_LN2 #define i11 0 #define i21 1 #define i22 2 #define epsilon 1e-15 /* Cauchy models */ #define CAUCHY_GAMMA 0 void Cauchy(double *x, model *cov, double *v){ double gamma = P0(CAUCHY_GAMMA); *v = POW(1.0 + *x * *x, -gamma); } void logCauchy(double *x, model *cov, double *v, double *Sign){ double gamma = P0(CAUCHY_GAMMA); *v = LOG(1.0 + *x * *x) * -gamma; *Sign = 1.0; } void TBM2Cauchy(double *x, model *cov, double *v){ double gamma = P0(CAUCHY_GAMMA), y2, lpy2; y2 = *x * *x; lpy2 = 1.0 + y2; switch ((int) (gamma * 2.0 + 0.001)) {// ueber check sichergestellt case 1 : *v = 1.0 / lpy2; break; case 3 : *v = (1.0 - y2)/ (lpy2 * lpy2); break; case 5 : *v = (1.0-y2*(2.0+0.333333333333333333333*y2))/(lpy2*lpy2*lpy2); break; case 7 : lpy2 *= lpy2; *v = (1.0- y2*(3.0+y2*(1.0+0.2*y2)))/(lpy2 * lpy2); break; default : ERR("TBM2 for cauchy only possible for alpha=0.5 + k; k=0, 1, 2, 3 "); } } void DCauchy(double *x, model *cov, double *v){ double y=*x, gamma = P0(CAUCHY_GAMMA); *v = (-2.0 * gamma * y) * POW(1.0 + y * y, -gamma - 1.0); } void DDCauchy(double *x, model *cov, double *v){ double ha = *x * *x, gamma = P0(CAUCHY_GAMMA); *v = 2.0 * gamma * ((2.0 * gamma + 1.0) * ha - 1.0) * POW(1.0 + ha, -gamma - 2.0); } void InverseCauchy(double*x, model *cov, double *v){ double gamma = P0(CAUCHY_GAMMA); if (*x == 0.0) *v = RF_INF; else *v = SQRT(POW(*x, -1.0 / gamma) - 1.0); } int checkCauchy(model VARIABLE_IS_NOT_USED *cov){ RETURN_NOERROR; } void rangeCauchy(model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[CAUCHY_GAMMA] = 0.0; range->max[CAUCHY_GAMMA] = RF_INF; range->pmin[CAUCHY_GAMMA] = 0.09; range->pmax[CAUCHY_GAMMA] = 10.0; range->openmin[CAUCHY_GAMMA] = true; range->openmax[CAUCHY_GAMMA] = true; } void coinitCauchy(model VARIABLE_IS_NOT_USED *cov, localinfotype *li) { li->instances = 1; li->value[0] = 1.0; // q[CUTOFF_A] li->msg[0] = MSGLOCAL_JUSTTRY; } void DrawMixCauchy(model VARIABLE_IS_NOT_USED *cov, double *random) { //better GR 3.381.4 ?? !!!! // split into parts <1 and >1 ?> *random = -LOG(1.0 -UNIFORM_RANDOM); } double LogMixDensCauchy(double VARIABLE_IS_NOT_USED *x, double logV, model *cov) { double gamma = P0(CAUCHY_GAMMA); // da dort 1/w vorgezogen return 0.5 * (gamma - 1.0) * logV - 0.5 * lgammafn(gamma); } /** another Cauchy model */ #define CTBM_ALPHA 0 #define CTBM_BETA 1 #define CTBM_GAMMA 2 void Cauchytbm(double *x, model *cov, double *v){ double ha, alpha = P0(CTBM_ALPHA), beta = P0(CTBM_BETA), gamma = P0(CTBM_GAMMA), y=*x; if (y==0) { *v = 1.0; } else { ha = POW(y, alpha); *v = (1.0 + (1.0 - beta / gamma) * ha) * POW(1.0 + ha, -beta / alpha - 1.0); } } void DCauchytbm(double *x, model *cov, double *v){ double y= *x, ha, alpha = P0(CTBM_ALPHA), beta = P0(CTBM_BETA), gamma = P0(CTBM_GAMMA); if (y == 0.0) *v = 0.0; // WRONG VALUE, but multiplied else { // by zero anyway ha = POW(y, alpha - 1.0); *v = beta * ha * (-1.0 - alpha / gamma + ha * y * (beta / gamma - 1.0)) * POW(1.0 + ha * y, -beta /alpha - 2.0); } } void rangeCauchytbm(model *cov, range_type *range){ range->min[CTBM_ALPHA] = 0.0; range->max[CTBM_ALPHA] = 2.0; range->pmin[CTBM_ALPHA] = 0.05; range->pmax[CTBM_ALPHA] = 2.0; range->openmin[CTBM_ALPHA] = true; range->openmax[CTBM_ALPHA] = false; range->min[CTBM_BETA] = 0.0; range->max[CTBM_BETA] = RF_INF; range->pmin[CTBM_BETA] = 0.05; range->pmax[CTBM_BETA] = 10.0; range->openmin[CTBM_BETA] = true; range->openmax[CTBM_BETA] = true; range->min[CTBM_GAMMA] = (double) OWNLOGDIM(0); range->max[CTBM_GAMMA] = RF_INF; range->pmin[CTBM_GAMMA] = range->min[CTBM_GAMMA]; range->pmax[CTBM_GAMMA] = range->pmin[CTBM_GAMMA] + 10.0; range->openmin[CTBM_GAMMA] = false; range->openmax[CTBM_GAMMA] = true; } // spatially constant (matrix); #define CONSTANT_M 0 void kappaconstant(int i, model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc) { *nr = *nc = i == CONSTANT_M ? 0 : -1; } void constant(double VARIABLE_IS_NOT_USED *x, model *cov, double *v){ int vdimSq = VDIM0 * VDIM1; MEMCOPY(v, P(CONSTANT_M), vdimSq * sizeof(double)); } void nonstatconstant(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED *y, model *cov, double *v){ int vdimSq = VDIM0 * VDIM1; MEMCOPY(v, P(CONSTANT_M), vdimSq * sizeof(double)); } int checkconstant(model *cov) { int err; // anzahl listen elemente if (GLOBAL.internal.warn_constant) { GLOBAL.internal.warn_constant = false; warning("NOTE that the definition of 'RMconstant' has changed in version 3.0.61. Maybe 'RMfixcov' is the model your are looking for"); } VDIM0 = VDIM1 = cov->nrow[CONSTANT_M]; if (VDIM0 != VDIM1) RETURN_ERR(ERROR_MATRIX_SQUARE); if (equalsVariogram(OWNTYPE(0))) { // ACHTUNG wirklich genau VariogramType SERR("strange call"); } if (cov->q != NULL) { return (int) cov->q[0]; } else QALLOC(1); cov->q[0] = NOERROR; // check whether positive eigenvalue if (!Ext_is_positive_definite(P(CONSTANT_M), VDIM0)) { cov->monotone = MONOTONE; cov->ptwise_definite = pt_indef; if (isnowPosDef(cov)) return cov->q[0] = ERROR_MATRIX_POSDEF; } else { cov->monotone = COMPLETELY_MON; cov->ptwise_definite = pt_posdef; } cov->matrix_indep_of_x = true; int vdim = VDIM0, vdimP1 = vdim + 1; double *p = P(CONSTANT_M); for (int i=0; impp.maxheights[i] = p[i * vdimP1]; err = checkkappas(cov); RETURN_ERR(err); } void rangeconstant(model VARIABLE_IS_NOT_USED *cov, range_type *range){ // auch in rangec verwendet! int dim = VDIM0; if (isnowPosDef(cov)) { if (isnowTcf(cov)) { range->min[CONSTANT_M] = (int) (dim == 1); range->max[CONSTANT_M] = 1; range->pmin[CONSTANT_M] = range->min[CONSTANT_M]; range->pmax[CONSTANT_M] = 1; range->openmin[CONSTANT_M] = false; range->openmax[CONSTANT_M] = false; } else { range->min[CONSTANT_M] = dim == 1 ? 0 : RF_NEGINF; range->max[CONSTANT_M] = RF_INF; range->pmin[CONSTANT_M] = dim == 1 ? 0 : -1e10; range->pmax[CONSTANT_M] = 1e10; range->openmin[CONSTANT_M] = dim != 1; range->openmax[CONSTANT_M] = true; } } else { range->min[CONSTANT_M] = RF_NEGINF; range->max[CONSTANT_M] = RF_INF; range->pmin[CONSTANT_M] = - 1e10; range->pmax[CONSTANT_M] = 1e10; range->openmin[CONSTANT_M] = true; range->openmax[CONSTANT_M] = true; } } /* exponential model */ void exponential(double *x, model VARIABLE_IS_NOT_USED *cov, double *v){ *v = EXP(- *x); // printf(": %10g %10g\n", *x, *v); } void logexponential(double *x, model VARIABLE_IS_NOT_USED *cov, double *v, double *Sign){ *v = - *x; *Sign = 1.0; } void TBM2exponential(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { double y = *x; *v = 1.0; if (y!=0.0) *v = 1.0 - PIHALF * y * Ext_I0mL0(y); } void Dexponential(double *x, model VARIABLE_IS_NOT_USED *cov, double *v){ *v = - EXP(- *x); } void DDexponential(double *x, model VARIABLE_IS_NOT_USED *cov, double *v){ *v = EXP(-*x); } void Inverseexponential(double *x, model VARIABLE_IS_NOT_USED *cov, double *v){ *v = (*x == 0.0) ? RF_INF : -LOG(*x); } void nonstatLogInvExp(double *x, model *cov, double *left, double *right){ double z = *x <= 0.0 ? - *x : 0.0; int d, dim = PREVLOGDIM(0); for (d=0; dspec); if (PREVLOGDIM(0) <= 2) RETURN_NOERROR; cs->density = densityexponential; return search_metropolis(cov, s); } else if (hasSmithFrame(cov)) { //Inverseexponential(&GLOBAL.mpp. about_ zero, NULL, &(cov->mpp.refradius)); // R *= GLOBAL.mpp.radius_natscale_factor; /* if (cov->mpp.moments >= 1) { int xi, d; double i[3], dimfak, Rfactor, sum, dimHalf = 0.5 * D; dimfak = gammafn(D); for (xi=1; xi<=2; xi++) { R = xi * cov->mpp.refradius; // // http://de.wikipedia.org/wiki/Kugel i[xi] = dimfak * 2 * POW(M_PI / (double) (xi*xi), dimHalf) / gammafn(dimHalf); if (R < RF_INF) { // Gradstein 3.351 fuer n=d-1 //printf("here\n"); for (sum = 1.0, factor=1.0, d=1; dmpp.mM[1] = cov->mpp.mMplus[1] = i[1]; if (cov->mpp.moments >= 2) { cov->mpp.mM[2] = cov->mpp.mMplus[2] = i[2]; } } */ assert(cov->mpp.maxheights[0] == 1.0); if (cov->mpp.moments >= 1) { cov->mpp.mM[1]= cov->mpp.mMplus[1] = SurfaceSphere(dim - 1, 1.0) * gammafn(D); } } else if (hasRandomFrame(cov)) { RETURN_NOERROR; } else ILLEGAL_FRAME; RETURN_NOERROR; } void do_exp(model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) { } void spectralexponential(model *cov, gen_storage *S, double *e) { spectral_storage *s = &(S->Sspectral); int dim = PREVLOGDIM(0); if (dim <= 2) { double A = 1.0 - UNIFORM_RANDOM; E12(s, dim, SQRT(1.0 / (A * A) - 1.0), e); // printf("dim = %d %f %f %f\n", dim, A, e[0]); } else { metropolis(cov, S, e); } } int checkexponential(model *cov) { int dim = OWNLOGDIM(0); if (dim > 2) cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = PREF_NONE; if (dim != 2) cov->pref[Hyperplane] = PREF_NONE; RETURN_NOERROR; } int hyperexponential(double radius, double *center, double *rx, model VARIABLE_IS_NOT_USED *cov, bool simulate, double ** Hx, double ** Hy, double ** Hr) { // lenx : half the length of the rectangle // center : center of the rectangle // simulate=false: estimated number of lines returned; // simulate=true: number of simulated lines returned; // hx, hy : direction of line // hr : distance of the line from the origin // rectangular area where center gives the center // // the function expects scale = 1; double lambda, phi, lx, ly, *hx, *hy, *hr; long i, p, // q=0; // q = RF_NA; assert(false); int k, Err; assert(OWNLOGDIM(0)==2); // we should be in two dimensions // first, we simulate the lines for a rectangle with center (0,0) // and half the side length equal to lenx lx = rx[0]; ly = rx[1]; lambda = TWOPI * radius * 0.5; /* total, integrated, intensity */ // 0.5 in order to get scale 1 if (!simulate) return lambda < 999999 ? (int) lambda : 999999 ; assert(*Hx==NULL); assert(*Hy==NULL); assert(*Hr==NULL); p = (long) rpois(lambda); if ((hx=*Hx=(double *) MALLOC(sizeof(double) * (p + 8 * sizeof(int))))==NULL|| (hy=*Hy=(double *) MALLOC(sizeof(double) * (p + 8 *sizeof(int))))==NULL|| (hr=*Hr=(double *) MALLOC(sizeof(double) * (p + 8 * sizeof(int))))==NULL){ Err=ERRORMEMORYALLOCATION; goto ErrorHandling; } /* creating the lines; some of the lines are not relevant as they do not intersect the rectangle investigated --> k!=4 (it is checked if all the corners of the rectangle are on one side (?) ) */ for(i=0; iinstances = 1; li->value[0] = 1.0; // q[CUTOFF_A] li->msg[0] = MSGLOCAL_OK; } void ieinitExp(model VARIABLE_IS_NOT_USED *cov, localinfotype *li) { li->instances = 1; li->value[0] = 1.0; li->msg[0] = MSGLOCAL_OK; } void DrawMixExp(model VARIABLE_IS_NOT_USED *cov, double *random) { // GR 3.325: int_-infty^infty EXP(x^2/2 - b/x^2) = EXP(-\sqrt b) double x = GAUSS_RANDOM(1.0); *random = 1.0 / (x * x); } double LogMixDensExp(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED logV, model VARIABLE_IS_NOT_USED *cov) { // todo: check use of mixdens --- likely to programme it now completely differently return 0.0; } /* Gausian model */ void Gauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { *v = EXP(- *x * *x); //printf("x =%10g v=%10g\n", x[0], *v); } void logGauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v, double *Sign) { *v = - *x * *x; *Sign = 1.0; } void DGauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { double y = *x; *v = -2.0 * y * EXP(- y * y); } void DDGauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { double y = *x * *x; *v = (4.0 * y - 2.0)* EXP(- y); } void D3Gauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { double y = *x * *x; *v = *x * (12 - 8 * y) * EXP(- y); } void D4Gauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { double y = *x * *x; *v = ((16 * y - 48) * y + 12) * EXP(- y); } void InverseGauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *v) { *v = SQRT(-LOG(*x)); } void nonstatLogInvGauss(double *x, model VARIABLE_IS_NOT_USED *cov, double *left, double *right) { double z = 0.0; if (*x < 0.0) z =SQRT(- *x); int d, dim = PREVLOGDIM(0); for (d=0; dframe) { case PoissonGaussType :// optimierte density fuer den Gauss-Fall double invscale; addModel(newmodel, GAUSS, cov); addModel(newmodel, DOLLAR); kdefault(*newmodel, DSCALE, INVSQRTTWO); addModel(newmodel, TRUNCSUPPORT); InverseGauss(&GLOBAL.mpp.about_zero, cov, &invscale); kdefault(*newmodel, TRUNC_RADIUS, invscale); break; default : if (hasSmithFrame(cov)) { // optimierte density fuer den Gauss-Fall // crash(); addModel(newmodel, GAUSS_DISTR, cov); // to kdefault(*newmodel, GAUSS_DISTR_MEAN, 0.0); kdefault(*newmodel, GAUSS_DISTR_SD, INVSQRTTWO); RETURN_NOERROR; } ILLEGAL_FRAME_STRUCT; } RETURN_NOERROR; } double IntUdeU2_intern(int d, double R, double expMR2) { // int_0^R u^{d-1} EXP(-u^2) \D u if (d==0) return (pnorm(R, 0.0, INVSQRTTWO, true, false) - 0.5) * SQRTPI; else if (d == 1) return 0.5 * (1.0 - expMR2); return 0.5 * (expMR2 + (d - 1.0) * IntUdeU2_intern(d - 2, R, expMR2)); } double IntUdeU2(int d, double R) { // int_0^R u^{d-1} EXP(-u^2) \D u return IntUdeU2_intern(d, R, EXP(-R*R)); } int initGauss(model *cov, gen_storage *s) { // if (hasAnyFrame(cov)) RETURN_NOERROR; if (HAS_SPECTRAL_FRAME(cov)) { spec_properties *cs = &(s->spec); if (OWNLOGDIM(0) <= 2) RETURN_NOERROR; cs->density = densityGauss; return search_metropolis(cov, s); } else if (hasSmithFrame(cov)) { // int_{b(0,R) e^{-a r^2} dr = d b_d int_0^R r^{d-1} e^{-a r^2} dr // where a = 2.0 * xi / sigma^2 // here : 2.0 is a factor so that the mpp function leads to the // gaussian covariance model EXP(-x^2) // xi : 1 : integral () // 2 : integral ()^2 double R = RF_INF; int dim = OWNLOGDIM(0); assert(COVNR == GAUSS); assert(cov->mpp.maxheights[0] = 1.0); if (cov->mpp.moments >= 1) { cov->mpp.mM[1] = cov->mpp.mMplus[1] = SurfaceSphere(dim - 1, 1.0) * IntUdeU2(dim - 1, R); int i; for (i=2; i<=cov->mpp.moments; i++) { cov->mpp.mM[i] = cov->mpp.mM[1] * POW((double) i, -0.5 * dim); } } // cov->mpp.maxheights[0] = 1.0; } else if (hasRandomFrame(cov) || hasAnyPoissonFrame(cov)) { RETURN_NOERROR; } else ILLEGAL_FRAME; RETURN_NOERROR; } void do_Gauss(model VARIABLE_IS_NOT_USED *cov, gen_storage VARIABLE_IS_NOT_USED *s) { } void spectralGauss(model *cov, gen_storage *S, double *e) { spectral_storage *s = &(S->Sspectral); int dim = PREVLOGDIM(0); if (dim <= 2) { E12(s, dim, 2.0 * SQRT(-LOG(1.0 - UNIFORM_RANDOM)), e); } else { metropolis(cov, S, e); } } void DrawMixGauss(model VARIABLE_IS_NOT_USED *cov, double VARIABLE_IS_NOT_USED *random) { *random = 1.0; } double LogMixDensGauss(double VARIABLE_IS_NOT_USED *x, double VARIABLE_IS_NOT_USED logV, model VARIABLE_IS_NOT_USED *cov) { return 0.0; } /* void densGauss(double *x, model *cov, double *v) { int factor[MAXMPPDIM+1] = {0, 1 / M_SQRT_PI, INVPI, INVPI / M_SQRT_PI, INVPI * INVPI}, dim = cov->tsdim; *v = factor[dim] * EXP(- *x * *x); } */ /* void getMassGauss(double *a, model *cov, double *kappas, double *m) { int i, j, k, kold, dim = cov->tsdim; double val[MAXMPPDIM + 1], sqrt2pi = SQRT2 * SQRTPI, factor[MAXMPPDIM+1] = {1, SQRT(2) / M_SQRT_PI, 2 * INVPI, 2 * SQRT(2) * INVPI / M_SQRT_PI, 4 * INVPI * INVPI}; val[0] = 1.0; for (i=1; i<=dim; i++) val[i] = (2.0 * pnorm(SQRT2 * a[i], 0.0, 1.0, false, false) - 1.0) * M_SQRT_PI; for (k=kold=i=0; i 1.0) ? 0.0 : (alpha < 1.0) ? RF_NEGINF : -beta); } else { ha = POW(y, alpha - 1.0); *v = -beta * ha * POW(1.0 + ha * y, -beta / alpha - 1.0); } } void DDgeneralisedCauchy(double *x, model *cov, double *v){ double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), ha, y=*x; if (y ==0.0) { *v = ((alpha==1.0) ? beta * (beta + 1.0) : (alpha==2.0) ? -beta : (alpha < 1) ? RF_INF : RF_NEGINF); } else { ha = POW(y, alpha); *v = beta * ha / (y * y) * (1.0 - alpha + (1.0 + beta) * ha) * POW(1.0 + ha, -beta / alpha - 2.0); } } void D3generalisedCauchy(double *x, model *cov, double *v){ double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), ha, y=*x; if (y ==0.0) { *v = ((alpha==2.0) ? 0 : (alpha == 1)? -beta*(beta+1)*(beta+2) : (alpha < 1)? RF_NEGINF : RF_INF); } else { ha=POW(y, alpha); *v = -beta * ha / (y * y* y) * ( (beta + 1)*(beta + 2)*ha*ha - (3*beta + alpha + 4)*(alpha - 1)*ha + (alpha-1)*(alpha -2) ) * POW(1.0 + ha, -beta / alpha - 3.0); } } void D4generalisedCauchy(double *x, model *cov, double *v){ double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), y=*x; if (y ==0.0) { *v = ((alpha==2.0) ? 3*beta*(beta + 2) : (alpha == 1)? beta*(beta+1)*(beta+2)*(beta + 3) : (alpha < 1)? RF_INF : RF_NEGINF); } else { int ha=POW(y, alpha), haSq = ha * ha; *v = beta * ha / (y * y* y *y) * ( -(alpha-1) * (alpha-2)*(alpha-3) - (alpha - 1) * ( 4*alpha*beta + alpha * (alpha + 7) + 6*beta*beta + 22*beta + 18 ) * haSq + (alpha - 1) * (alpha * (4 * alpha + 7 * beta + 4) - 11*beta - 18) * ha + (beta + 1) * (beta +2) * (beta + 3) * ha * haSq) * POW(1.0 + ha, -beta / alpha - 4.0); } } void loggeneralisedCauchy(double *x, model *cov, double *v, double *Sign){ double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA); *v = LOG(1.0 + POW(*x, alpha)) * -beta/alpha; *Sign = 1.0; } void InversegeneralisedCauchy(double *x, model *cov, double *v) { double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA); *v = RF_INF; if (*x != 0.0) *v = POW(POW(*x, -alpha / beta) - 1.0, 1.0 / alpha); // MLE works much better with 0.01 then with 0.05 } int checkgeneralisedCauchy(model *cov){ if (OWNLOGDIM(0) > 2) cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = PREF_NONE; cov->monotone = P0(GENC_ALPHA) <= 1.0 ? COMPLETELY_MON : NORMAL_MIXTURE; RETURN_NOERROR; } void rangegeneralisedCauchy(model *cov, range_type *range) { bool tcf = isnowTcf(cov) || equalsSphericalIsotropic(OWNISO(0)); range->min[GENC_ALPHA] = 0.0; range->max[GENC_ALPHA] = tcf ? 1.0 : 2.0; range->pmin[GENC_ALPHA] = 0.05; range->pmax[GENC_ALPHA] = range->max[GENC_ALPHA]; range->openmin[GENC_ALPHA] = true; range->openmax[GENC_ALPHA] = false; range->min[GENC_BETA] = 0.0; range->max[GENC_BETA] = RF_INF; range->pmin[GENC_BETA] = 0.05; range->pmax[GENC_BETA] = 10.0; range->openmin[GENC_BETA] = true; range->openmax[GENC_BETA] = true; } void coinitgenCauchy(model *cov, localinfotype *li) { double thres[2] = {0.5, 1.0}, alpha=P0(GENC_ALPHA); if (alpha <= thres[0]) { li->instances = 2; li->value[0] = 0.5; li->value[1] = 1.0; li->msg[0] = li->msg[1] = MSGLOCAL_OK; } else { li->instances = 1; li->value[0] = 1.0; // q[CUTOFF_A] li->msg[0] = (alpha <= thres[1]) ? MSGLOCAL_OK : MSGLOCAL_JUSTTRY; } } void ieinitgenCauchy(model *cov, localinfotype *li) { li->instances = 1; if (P0(GENC_ALPHA) <= 1.0) { li->value[0] = 1.0; li->msg[0] = MSGLOCAL_OK; } else { li->value[0] = 1.5; li->msg[0] = MSGLOCAL_JUSTTRY; } } /*---------------------------------*/ // ************************* bivariate Cauchy void kappa_biCauchy(int i, model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){ *nc = 1; *nr = i == BICauchyalpha || i == BICauchybeta|| i == BICauchyscale ? 3 : 1; } int checkbiCauchy(model *cov) { if (OWNLOGDIM(0) > 2) cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = PREF_NONE; RETURN_NOERROR; } void rangebiCauchy(model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[BICauchyalpha] = 0.0; range->max[BICauchyalpha] = 1; range->pmin[BICauchyalpha] = 0.05; range->pmax[BICauchyalpha] = 1; range->openmin[BICauchyalpha] = true; range->openmax[BICauchyalpha] = true; range->min[BICauchybeta] = 0.0; range->max[BICauchybeta] = RF_INF; range->pmin[BICauchybeta] = 0.05; range->pmax[BICauchybeta] = 10.0; range->openmin[BICauchybeta] = true; range->openmax[BICauchybeta] = true; range->min[BICauchyscale] = 0.0; range->max[BICauchyscale] = RF_INF; range->pmin[BICauchyscale] = 1e-2; range->pmax[BICauchyscale] = 100.0; range->openmin[BICauchyscale] = true; range->openmax[BICauchyscale] = true; // to do: check rhos range->min[BICauchyrho] = - 1; range->max[BICauchyrho] = 1; range->pmin[BICauchyrho] = - 1; range->pmax[BICauchyrho] = 1; range->openmin[BICauchyrho] = false; range->openmax[BICauchyrho] = false; } void coinitbiCauchy(model *cov, localinfotype *li) { double thres = 1, *alpha = P(BICauchyalpha); if ( (alpha[0] <= thres) && (alpha[1] <= thres) && (alpha[2] <= thres) ) { li->instances = 1; li->value[0] = CUTOFF_THIRD_CONDITION; // q[CUTOFF_A] li->msg[0] =MSGLOCAL_OK; } else { li->instances = 1; li->value[0] = CUTOFF_THIRD_CONDITION ; // q[CUTOFF_A] li->msg[0] = MSGLOCAL_JUSTTRY; } } void biCauchy (double *x, model *cov, double *v) { int i; double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), y = *x, z; assert(BICauchyalpha == GENC_ALPHA); assert(BICauchybeta == GENC_BETA); for (i=0; i<3; i++) { z = y/P(BICauchyscale)[i]; P(GENC_ALPHA)[0] = P(BICauchyalpha)[i]; P(GENC_BETA)[0] = P(BICauchybeta)[i]; generalisedCauchy(&z, cov, v + i); } P(BICauchyalpha)[0] = alpha; P(BICauchybeta)[0] = beta; v[3] = v[2]; v[1] *= P0(BICauchyrho); v[2] = v[1]; } void DbiCauchy(double *x, model *cov, double *v) { int i; double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), y = *x, z; assert(BICauchyalpha == GENC_ALPHA); assert(BICauchybeta == GENC_BETA); for (i=0; i<3; i++) { z = y/P(BICauchyscale)[i]; P(GENC_ALPHA)[0] = P(BICauchyalpha)[i]; P(GENC_BETA)[0] = P(BICauchybeta)[i]; DgeneralisedCauchy(&z, cov, v + i); v[i] /=P(BICauchyscale)[i]; } P(BICauchyalpha)[0] = alpha; P(BICauchybeta)[0] = beta; v[3] = v[2]; v[1] *= P0(BICauchyrho); v[2] = v[1]; } void DDbiCauchy(double *x, model *cov, double *v) { int i; double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), y = *x, z; assert(BICauchyalpha == GENC_ALPHA); assert(BICauchybeta == GENC_BETA); for (i=0; i<3; i++) { z = y/P(BICauchyscale)[i]; P(GENC_ALPHA)[0] = P(BICauchyalpha)[i]; P(GENC_BETA)[0] = P(BICauchybeta)[i]; DDgeneralisedCauchy(&z, cov, v + i); v[i] /=(P(BICauchyscale)[i]*P(BICauchyscale)[i]); } P(BICauchyalpha)[0] = alpha; P(BICauchybeta)[0] = beta; v[3] = v[2]; v[1] *= P0(BICauchyrho); v[2] = v[1]; } void D3biCauchy(double *x, model *cov, double *v) { int i; double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), y = *x, z; assert(BICauchyalpha == GENC_ALPHA); assert(BICauchybeta == GENC_BETA); for (i=0; i<3; i++) { z = y/P(BICauchyscale)[i]; P(GENC_ALPHA)[0] = P(BICauchyalpha)[i]; P(GENC_BETA)[0] = P(BICauchybeta)[i]; D3generalisedCauchy(&z, cov, v + i); v[i] /=(P(BICauchyscale)[i]*P(BICauchyscale)[i]*P(BICauchyscale)[i]); } P(BICauchyalpha)[0] = alpha; P(BICauchybeta)[0] = beta; v[3] = v[2]; v[1] *= P0(BICauchyrho); v[2] = v[1]; } void D4biCauchy(double *x, model *cov, double *v) { int i; double alpha = P0(GENC_ALPHA), beta=P0(GENC_BETA), y = *x, z; assert(BICauchyalpha == GENC_ALPHA); assert(BICauchybeta == GENC_BETA); for (i=0; i<3; i++) { z = y/P(BICauchyscale)[i]; P(GENC_ALPHA)[0] = P(BICauchyalpha)[i]; P(GENC_BETA)[0] = P(BICauchybeta)[i]; D4generalisedCauchy(&z, cov, v + i); double dummy = P(BICauchyscale)[i]*P(BICauchyscale)[i]; v[i] /=(dummy*dummy); } P(BICauchyalpha)[0] = alpha; P(BICauchybeta)[0] = beta; v[3] = v[2]; v[1] *= P0(BICauchyrho); v[2] = v[1]; } /* epsC -> generalised Cauchy : leading 1 is now an eps */ #define EPS_ALPHA 0 #define EPS_BETA 1 #define EPS_EPS 2 void epsC(double *x, model *cov, double *v){ double alpha = P0(EPS_ALPHA), beta=P0(EPS_BETA), eps=P0(EPS_EPS); *v = POW(eps + POW(*x, alpha), -beta/alpha); } void logepsC(double *x, model *cov, double *v, double *Sign){ double alpha = P0(EPS_ALPHA), beta=P0(EPS_BETA), eps=P0(EPS_EPS); *v = LOG(eps + POW(*x, alpha)) * -beta/alpha; *Sign = 1.0; } void DepsC(double *x, model *cov, double *v){ double ha, y=*x, alpha = P0(EPS_ALPHA), beta=P0(EPS_BETA), eps=P0(EPS_EPS); if (y ==0.0) { *v = (eps == 0.0) ? RF_NEGINF : ((alpha > 1.0) ? 0.0 : (alpha < 1.0) ? RF_NEGINF : -beta); } else { ha = POW(y, alpha - 1.0); *v = -beta * ha * POW(eps + ha * y, -beta / alpha - 1.0); } } void DDepsC(double *x, model *cov, double *v){ double ha, y=*x, alpha = P0(EPS_ALPHA), beta=P0(EPS_BETA), eps=P0(EPS_EPS); if (y ==0.0) { *v = (eps == 0.0) ? RF_INF : ((alpha==2.0) ? beta * (beta + 1.0) : RF_INF); } else { ha = POW(y, alpha); *v = beta * ha / (y * y) * ( (1.0 - alpha) * eps + (1.0 + beta) * ha) * POW(eps + ha, -beta / alpha - 2.0); } } void InverseepsC(double *x, model *cov, double *v){ double alpha = P0(EPS_ALPHA), beta=P0(EPS_BETA), eps=P0(EPS_EPS); *v = RF_INF; if (*x != 0.0) *v = POW(POW(*x, -alpha / beta) - eps, 1.0 / alpha); } int checkepsC(model *cov){ double eps=P0(EPS_ALPHA); int i, err; if (OWNLOGDIM(0) > 2) cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = PREF_NONE; if ((err = checkkappas(cov, false)) != NOERROR) RETURN_ERR(err); kdefault(cov, EPS_ALPHA, 1.0); kdefault(cov, EPS_BETA, 1.0); kdefault(cov, EPS_EPS, 0.0); if (ISNAN(eps) || eps == 0.0) { // OWNDOM(0)=GENERALISEDCOVARIANCE; // later for (i=CircEmbed; ipref[i] = PREF_NONE; } RETURN_NOERROR; } void rangeepsC(model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[EPS_ALPHA] = 0.0; range->max[EPS_ALPHA] = 2.0; range->pmin[EPS_ALPHA] = 0.05; range->pmax[EPS_ALPHA] = 2.0; range->openmin[EPS_ALPHA] = true; range->openmax[EPS_ALPHA] = false; range->min[EPS_BETA] = 0.0; range->max[EPS_BETA] = RF_INF; range->pmin[EPS_BETA] = 0.05; range->pmax[EPS_BETA] = 10.0; range->openmin[EPS_BETA] = true; range->openmax[EPS_BETA] = true; range->min[EPS_EPS] = 0.0; range->max[EPS_EPS] = RF_INF; range->pmin[EPS_EPS] = 0.0; range->pmax[EPS_EPS] = 10.0; range->openmin[EPS_EPS] = true; // false for generlised covariance range->openmax[EPS_EPS] = true; } /* hyperbolic */ #define BOLIC_NU 0 #define BOLIC_XI 1 #define BOLIC_DELTA 2 void hyperbolic(double *x, model *cov, double *v){ double Sign; loghyperbolic(x, cov, v, &Sign); *v = EXP(*v); } void loghyperbolic(double *x, model *cov, double *v, double *Sign){ double nu = P0(BOLIC_NU), xi=P0(BOLIC_XI), delta=P0(BOLIC_DELTA), y=*x; *Sign = 1.0; if (y==0.0) { *v = 0.0; return; } else if (y==RF_INF) { *v = RF_NEGINF; *Sign = 0.0; return; } if (delta==0) { // whittle matern if (nu > 80) warning("extremely imprecise results due to nu>80"); *v = logWM(y * xi, nu, cov->q[WM_LOGGAMMA], 0.0); } else if (xi==0) { //cauchy => NU2 < 0 ! y /= delta; /* note change in Sign as NU2<0 */ *v = LOG(1.0 + y * y) * 0.5 * nu; } else { y=SQRT(delta * delta + y * y); double bk[MATERN_NU_THRES + 1L], xiy = xi * y, nuThres = nu <= MATERN_NU_THRES ? nu : MATERN_NU_THRES, factor = 1 / SQRT(nuThres); if (xiy <= LOW_MATERN) {*v = 1.0; return;} if (xiy == RF_INF) {*v = 0.0; return;} // also for cosine i.e. nu=-1/2 *v = QVALUE3 + nu * LOG(y) + LOG(bessel_k_ex(xiy, nu, 2.0, bk)) - xiy; if (nu > MATERN_NU_THRES) { // factor!=0.0 && // printf("UU \n"); double w, g = MATERN_NU_THRES / nu; y = 0.5 * xiy * factor; Gauss(&y, NULL, &w); *v = *v * g + (1.0 - g) * w; } } } void Dhyperbolic(double *x, model *cov, double *v) { double nu = P0(BOLIC_NU), xi=P0(BOLIC_XI), delta=P0(BOLIC_DELTA); double y = *x; if (y==0.0) { *v = 1.0; return; } if (delta==0) { // whittle matern *v = xi * Ext_DWM(y * xi, nu, cov->q[WM_LOGGAMMA], 0.0); *v *= xi; } else if (xi==0) { //cauchy double ha; y /= delta; ha = y * y; /* note change in Sign as NU2<0 */ *v = nu * FABS(y) * POW(1.0 + ha, 0.5 * nu - 1.0) / delta; } else { double nuThres = nu <= MATERN_NU_THRES ? nu : MATERN_NU_THRES, s=SQRT(delta * delta + y * y), xi_s = xi * s, logs = LOG(s), bk[MATERN_NU_THRES + 1L]; *v = - y * xi * EXP(QVALUE3 + (nuThres-1.0) * logs + LOG(bessel_k_ex(xi_s, nuThres-1.0, 2.0, bk)) - xi_s); if (nu > MATERN_NU_THRES) { double w, g = MATERN_NU_THRES / nu, factor = 1.0 / SQRT(nuThres), scale = 0.5 * factor; y = xi_s * scale; DGauss(&y, NULL, &w); w *= scale; *v = *v * g + (1.0 - g) * w; } } } int inithyperbolic(model *cov, gen_storage VARIABLE_IS_NOT_USED *s) { double nu = P0(BOLIC_NU), nuThres = nu <= MATERN_NU_THRES ? nu : MATERN_NU_THRES, xi=P0(BOLIC_XI), delta=P0(BOLIC_DELTA), xd = xi * delta, // xidelta bk[MATERN_NU_THRES + 1L]; QVALUE3 = xd - LOG(bessel_k_ex(xd, nuThres, 2.0, bk)) - nuThres * LOG(delta); if (nu > MATERN_NU_THRES) { // factor!=0.0 && // printf("UU \n"); double w, factor = 1 / SQRT(nuThres), g = MATERN_NU_THRES / nu, y = 0.5 * xd * factor; Gauss(&y, NULL, &w); QVALUE3 = QVALUE3 * g + (1.0 - g) * w; } if (!ISNA(delta) && delta == 0.0 && !ISNA(nu)) { assert(cov->q + WM_LOGGAMMA == &(QVALUE)); assert(cov->q + WM_GAMMA == &(QVALUE2)); cov->q[WM_LOGGAMMA] = lgammafn(nuThres); // QVALUE cov->q[WM_GAMMA] = gammafn(nuThres); // QVALUE2 } RETURN_NOERROR; } int checkhyperbolic(model *cov){ double nu = P0(BOLIC_NU), xi=P0(BOLIC_XI), delta=P0(BOLIC_DELTA); int i; for (i=0; i<= Nothing; i++) cov->pref[i] *= ( ((bool) ISNAN(nu)) || nu < WhittleUpperNu[i]); if (nu>0) { if ((delta<0) || (xi<=0)) { SERR3("xi>0 and delta>=0 if nu>0. Got nu=%10g and delta=%10g for xi=%10g", nu, delta, xi); } } else if (nu<0) { if ((delta<=0) || (xi<0)) { SERR3("xi>=0 and delta>0 if nu<0. Got nu=%10g and delta=%10g for xi=%10g", nu, delta, xi); } } else { // nu==0.0 if ((delta<=0) || (xi<=0)) { SERR3("xi>0 and delta>0 if nu=0. Got nu=%10g and delta=%10g for xi=%10g", nu, delta, xi); } } if (cov->q == NULL) { EXTRA_Q; inithyperbolic(cov, NULL); } RETURN_NOERROR; } void rangehyperbolic(model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[BOLIC_NU] = RF_NEGINF; range->max[BOLIC_NU] = RF_INF; range->pmin[BOLIC_NU] = -20.0; range->pmax[BOLIC_NU] = 20.0; range->openmin[BOLIC_NU] = true; range->openmax[BOLIC_NU] = true; int i; for (i=1; i<=2; i++) { range->min[i] = 0.0; range->max[i] = RF_INF; range->pmin[i] = 0.000001; range->pmax[i] = 10.0; range->openmin[i] = false; range->openmax[i] = true; } } /* stable model */ #define STABLE_ALPHA 0 void stable(double *x, model *cov, double *v){ double y = *x, alpha = P0(STABLE_ALPHA); *v = 1.0; if (y!=0.0) *v = EXP(-POW(y, alpha)); } void logstable(double *x, model *cov, double *v, double *Sign){ double y = *x, alpha = P0(STABLE_ALPHA); *v = 0.0; if (y!=0.0) *v= -POW(y, alpha); *Sign = 1.0; } void Dstable(double *x, model *cov, double *v){ double z, y = *x, alpha = P0(STABLE_ALPHA); if (y == 0.0) { *v = (alpha > 1.0) ? 0.0 : (alpha < 1.0) ? RF_NEGINF : -1.0; } else { z = POW(y, alpha - 1.0); *v = -alpha * z * EXP(-z * y); } } /* stable: second derivative at t=1 */ void DDstable(double *x, model *cov, double *v) { double z, y = *x, alpha = P0(STABLE_ALPHA), xalpha; if (y == 0.0) { // olga please check *v = (alpha == 1.0) ? 1.0 : (alpha == 2.0) ? -2.0 : //alpha * (1 - alpha) alpha < 1.0 ? RF_INF : RF_NEGINF; } else { z = POW(y, alpha - 2.0); xalpha = z * y * y; *v = alpha * (1.0 - alpha + alpha * xalpha) * z * EXP(-xalpha); } } void D3stable(double *x, model *cov, double *v) { double z, y = *x, alpha = P0(STABLE_ALPHA), xalpha; if (y == 0.0) { // olga, please check; please note that RF_INF has been incorrect (reserved for integer variables only) *v = (alpha == 1.0) ? -1.0 : (alpha == 2.0) ? 0 : alpha < 1.0 ? RF_NEGINF : RF_INF; } else { z = POW(y, alpha - 3.0); xalpha = z * y * y * y; *v = -alpha * ( 3*alpha*(xalpha -1) + alpha*alpha* (xalpha*xalpha - 3*xalpha + 1) + 2) * z * EXP(-xalpha); } } void D4stable(double *x, model *cov, double *v) { double z, y = *x, alpha = P0(STABLE_ALPHA), xalpha; if (y == 0.0) { *v = (alpha == 1.0) ? 1.0 : (alpha == 2.0) ? 0 : alpha < 1.0 ? RF_INF : RF_NEGINF; } else { z = POW(y, alpha - 4.0); xalpha = z * y * y * y * y; *v = alpha*( alpha*alpha*alpha*(7*xalpha -6*xalpha*xalpha + xalpha*xalpha*xalpha - 1 ) + +6*alpha*alpha*(-3*xalpha + xalpha*xalpha +1 )+ 11*alpha*(xalpha - 1 ) + 6 ) * z * EXP(-xalpha); } } void D5stable(double *x, model *cov, double *v) { double z, y = *x, alpha = P0(STABLE_ALPHA), xalpha; if (y == 0.0) { *v = (alpha == 1.0) ? -1.0 : (alpha == 2.0) ? 0 : alpha < 1.0 ? RF_NEGINF : RF_INF;; } else { z = POW(y, alpha - 5.0); xalpha = z * y * y * y * y * y; *v = -alpha*( POW(alpha, 4)*(-15*xalpha +25*xalpha*xalpha - 10*POW(xalpha, 3) + POW(xalpha, 4)+ 1 ) + 10*alpha*alpha*alpha*(7*xalpha -6*xalpha*xalpha + POW(xalpha, 3) - 1 )+ 35*alpha*alpha*(-3*xalpha + xalpha*xalpha + 1 )+ 50*alpha*(xalpha -1 ) + 24 )* z * EXP(-xalpha); } } void Inversestable(double *x, model *cov, double *v){ double y = *x, alpha = P0(STABLE_ALPHA); *v = y>1.0 ? 0.0 : y == 0.0 ? RF_INF : POW( - LOG(y), 1.0 / alpha); } void nonstatLogInversestable(double *x, model *cov, double *left, double *right){ double alpha = P0(STABLE_ALPHA), z = *x <= 0.0 ? POW( - *x, 1.0 / alpha) : 0.0; int d, dim = OWNLOGDIM(0); for (d=0; d 2) cov->pref[CircEmbedCutoff] = cov->pref[CircEmbedIntrinsic] = PREF_NONE; if (alpha == 2.0) cov->pref[CircEmbed] = 2; cov->monotone = alpha <= 1.0 ? COMPLETELY_MON : NORMAL_MIXTURE; RETURN_NOERROR; } void rangestable(model *cov, range_type *range){ bool tcf = isnowTcf(cov) || equalsSphericalIsotropic(OWNISO(0)); range->min[STABLE_ALPHA] = 0.0; range->max[STABLE_ALPHA] = tcf ? 1.0 : 2.0; range->pmin[STABLE_ALPHA] = 0.06; range->pmax[STABLE_ALPHA] = range->max[STABLE_ALPHA]; range->openmin[STABLE_ALPHA] = true; range->openmax[STABLE_ALPHA] = false; } void coinitstable(model *cov, localinfotype *li) { coinitgenCauchy(cov, li); } void ieinitstable(model *cov, localinfotype *li) { ieinitgenCauchy(cov, li); } /* DOUBLEISOTROPIC stable model for testing purposes only */ void stableX(double *x, model *cov, double *v){ double y, alpha = P0(STABLE_ALPHA); y = x[0] * x[0] + x[1] * x[1]; *v = 1.0; if (y!=0.0) *v = EXP(-POW(y, 0.5 * alpha)); } void DstableX(double *x, model *cov, double *v){ double z, y, alpha = P0(STABLE_ALPHA); y = x[0] * x[0] + x[1] * x[1]; if (y == 0.0) { *v = ((alpha > 1.0) ? 0.0 : (alpha < 1.0) ? RF_INF : 1.0); } else { z = POW(y, 0.5 * alpha - 1.0); *v = -alpha * x[0] * z * EXP(- z * y); } } /* END DOUBLEISOTROPIC stable model for testing purposes only */ // ************************* bivariate powered exponential or bivariate stable void kappa_biStable(int i, model VARIABLE_IS_NOT_USED *cov, int *nr, int *nc){ *nc = 1; *nr = ( i == BIStablealpha || i == BIStablescale ) ? 3 : (i == BIStablecdiag || i == BIStablealphadiag ) ? 2 : (i == BIStablerho || i == BIStablebetared || i == BIStablerhored ) ? 1 : -1; } int checkbiStable(model *cov) { int err; gen_storage s; gen_NULL(&s); s.check = true; if ((err = checkkappas(cov, false)) != NOERROR) RETURN_ERR(err); // bistable_storage *S = cov->Sbistable; if (cov->Sbistable == NULL) { ONCE_NEW_STORAGE(bistable); bistable_storage *S = cov->Sbistable; S->alphadiag_given = !PisNULL(BIStablealphadiag); S->rhored_given = !PisNULL(BIStablerhored); // printf("check: alphadiag_given = %d\n", S->alphadiag_given); } if ((err=initbiStable(cov, &s)) != NOERROR) RETURN_ERR(err); VDIM0 = VDIM1 = 2; if ((err = checkkappas(cov)) != NOERROR) RETURN_ERR(err); //PMI(cov); RETURN_NOERROR; } // alphadiag always out sortsofparam sortof_bistable(model *cov, int k, int VARIABLE_IS_NOT_USED *row, int VARIABLE_IS_NOT_USED *col, sort_origin origin) { bistable_storage *S = cov->Sbistable; if (S == NULL) return UNKNOWNPARAM; switch(k) { case BIStablealphadiag : case BIStablebetared : return S->alphadiag_given || origin != original ? ANYPARAM : IGNOREPARAM; case BIStablescale: return SCALEPARAM; case BIStablecdiag : return VARPARAM; case BIStablerho : return S->rhored_given || origin == mle_conform ? IGNOREPARAM : ONLYRETURN; case BIStablerhored : return S->rhored_given || origin != original ? ANYPARAM : IGNOREPARAM; case BIStablealpha : return S->alphadiag_given || origin == mle_conform ? IGNOREPARAM : ONLYRETURN; default : BUG; } } sortsofparam sortof_bistable_INisOUT(model *cov, int k, int VARIABLE_IS_NOT_USED *row, int VARIABLE_IS_NOT_USED *col) { bistable_storage *S = cov->Sbistable; if (S == NULL) return UNKNOWNPARAM; switch(k) { case BIStablealphadiag : return !S->alphadiag_given ? ONLYMLE : ANYPARAM; case BIStablescale: return SCALEPARAM; case BIStablecdiag : return VARPARAM; case BIStablerho : return S->rhored_given ? ONLYRETURN : IGNOREPARAM; case BIStablerhored : return S->rhored_given ? ONLYMLE : ANYPARAM; case BIStablebetared : return !S->alphadiag_given ? ONLYMLE : ANYPARAM; case BIStablealpha : return !S->alphadiag_given ? ONLYRETURN : IGNOREPARAM; default : BUG; } } void rangebiStable(model VARIABLE_IS_NOT_USED *cov, range_type *range){ range->min[BIStablealpha] = 0.0; range->max[BIStablealpha] = 1; range->pmin[BIStablealpha] = 0.06; range->pmax[BIStablealpha] = 1; range->openmin[BIStablealpha] = true; range->openmax[BIStablealpha] = false; range->min[BIStablealphadiag] = 0.0; range->max[BIStablealphadiag] = 1; range->pmin[BIStablealphadiag] = 0.06; range->pmax[BIStablealphadiag] = 1; range->openmin[BIStablealphadiag] = true; range->openmax[BIStablealphadiag] = false; range->min[BIStablebetared] = 0.0; range->max[BIStablebetared] = 1.0; range->pmin[BIStablebetared] = 0.0; range->pmax[BIStablebetared] =1.0; range->openmin[BIStablebetared] = false; range->openmax[BIStablebetared] = false; range->min[BIStablescale] = 0.0; range->max[BIStablescale] = RF_INF; range->pmin[BIStablescale] = 1e-2; range->pmax[BIStablescale] = 100.0; range->openmin[BIStablescale] = true; range->openmax[BIStablescale] = true; // *c = P0(BIc]; // to do: check rhos // int dim = cov->tsdim; // double *scale = P(BIStablescale); range->min[BIStablerho] = - 1.0; range->max[BIStablerho] = 1.0; range->pmin[BIStablerho] = - 1.0; range->pmax[BIStablerho] = 1.0; range->openmin[BIStablerho] = false; range->openmax[BIStablerho] = false; range->min[BIStablerhored] = - 1.0; range->max[BIStablerhored] = 1.0; range->pmin[BIStablerhored] = - 1.0; range->pmax[BIStablerhored] = 1.0; range->openmin[BIStablerhored] = false; range->openmax[BIStablerhored] = false; range->min[BIStablecdiag] = 0; range->max[BIStablecdiag] = RF_INF; range->pmin[BIStablecdiag] = 0.001; range->pmax[BIStablecdiag] = 1000; range->openmin[BIStablecdiag] = false; range->openmax[BIStablecdiag] = true; } void biStablePolynome(double r, double alpha, double a, int dim, double *v ) { double x = POW(a*r, alpha); if (dim == 1) { *v = alpha*x - alpha + 1; } if ( dim == 2 || dim == 3 ) { *v = alpha*alpha*x*x - 3*alpha*alpha*x + 4*alpha*x + alpha*alpha - 4*alpha + 3; } } void biStableUnderInfLog(double r, double *alpha, double *a, int dim, double *res ) { double x = POW(a[i11]*r, alpha[i11]), y = POW(a[i21]*r, alpha[i21]), z = POW(a[i22]*r, alpha[i22]); double p1, p2, p3; biStablePolynome(r, alpha[i11], a[i11], dim, &p1 ); biStablePolynome(r, alpha[i21], a[i21], dim, &p2 ); biStablePolynome(r, alpha[i22], a[i22], dim, &p3 ); if (r == 0) { *res = 0; } else { *res = (alpha[i11] + alpha[i22] - 2*alpha[i21])*LOG(r) + (2*y - x - z) + LOG(p1*p3/(p2*p2)); } } void biStableInterval(double *alpha, double *a, int dim, double *left, double *right ) { double middle = 1, fmiddle, fright, fleft, scalesratio1, scalesratio2, firstmin = 1.0/epsilon; *left = *right = middle; // check first if POW(a[0]/a[1], alpha[0] > 11 or // check first if POW(a[2]/a[1], alpha[2] > 11 scalesratio1 = POW(a[0]/a[1], alpha[0]); scalesratio2 = POW(a[2]/a[1], alpha[2]); if ( (scalesratio1 >= 11) || (scalesratio2 >= 11) ) { biStableUnderInfLog( 1.0/(POW(2, 1/alpha[1])*a[1] ), alpha, a, dim, &firstmin); if (EXP(firstmin) < epsilon ) { *left = 0; *right = 0; return; } } biStableUnderInfLog(middle, alpha, a, dim, &fmiddle); // if the value at 1/(POW(2, 1/alpha[1])*a[1] ) is smaller that the value at 1 // start searching around 1/(POW(2, 1/alpha[1])*a[1] ) if ( fmiddle > firstmin ) { middle = 1.0/(POW(2, 1/alpha[1])*a[1] ); *left = *right = middle; fmiddle = firstmin; } fright = fleft = fmiddle; while (( fmiddle >= MIN(fleft, fright)) && (EXP(MIN(fleft, MIN(fright, fmiddle))) > epsilon ) ) { if ( fleft <= fmiddle ) { middle = *left; fmiddle = fleft; *left = *left/2; } if (fright <= fmiddle) { middle = *right; fmiddle = fright; *right = *right*2; } biStableUnderInfLog(*right, alpha, a, dim, &fright ); biStableUnderInfLog(*left, alpha, a, dim, &fleft ); } if (EXP(MIN(fleft, MIN(fright, fmiddle))) <= epsilon ) { *left = 0; *right = 0; } } /* * Golden search algorithm from Numerical Recipes in C, * book by William H. Press * * (c, c) - interval for searching * *alpha, *a - array of parameters of biStable model * *rhomax - maximum allowable rho for *alpha, *a, the result of * the searching * dim - dimension of the model * */ #define GOLDENR 0.61803399 #define GOLDENC (1.0 -GOLDENR) #define GOLDENTOL 1e-9 #define GOLDENSHIFT2(a, b, c) (a) = (b); (b) = (c); #define GOLDENSHIFT3(a, b, c, d) (a) = (b); (b) = (c); (c) = (d); void biStableMinRho(model *cov, double *alpha, double *a, double ax, double cx, double *rhomax) { double bx = ax + (cx - ax)*GOLDENC, f1, f2, x0, x1, x2, x3, dummy, logconst = LOG(alpha[i11]*alpha[i22]/(alpha[i21]*alpha[i21])* POW(a[i11], alpha[i11])*POW(a[i22], alpha[i22])/ POW(a[i21], 2*alpha[i21])); x0 = ax; x3 = cx; if ( FABS(cx - bx) > FABS(bx - ax) ) { x1 = bx; x2 = bx + GOLDENC*(cx - bx); } else { x2 = bx; x1 = bx - GOLDENC*(bx - ax); } int dim = OWNLOGDIM(0); biStableUnderInfLog(x1, alpha, a, dim, &f1); biStableUnderInfLog(x2, alpha, a, dim, &f2); while ( FABS(x3 - x0) > GOLDENTOL*(FABS(x1) + FABS(x2) ) ) { if (f2 < f1) { GOLDENSHIFT3(x0, x1, x2, GOLDENR*x1 + GOLDENC*x3 ) biStableUnderInfLog(x2, alpha, a, dim, &dummy); GOLDENSHIFT2(f1, f2, dummy ) } else { GOLDENSHIFT3(x3, x2, x1, GOLDENR*x2 + GOLDENC*x0 ) biStableUnderInfLog(x1, alpha, a, dim, &dummy); GOLDENSHIFT2(f2, f1, dummy ) } } if (f1 < f2) { *rhomax = SQRT(EXP(f1+logconst)); } else { *rhomax = SQRT(EXP(f2+logconst)); } *rhomax = MIN(*rhomax, 1); } int initbiStable(model *cov, gen_storage VARIABLE_IS_NOT_USED *s) { double a[3], *cdiag = P(BIStablecdiag), *rho = P(BIStablerho), *rhored = P(BIStablerhored), // *cdiag = P(BIStablecdiag), *alpha = P(BIStablealpha), betared = RF_NAN, *alphadiag, *scale = P(BIStablescale), betaa, betac, //beta_red in [0, 1]; beta_red*betaa + betac in [betac, betaa + betac] rhomax = -2, left = 0, right = 0; int dim = OWNLOGDIM(0); bool notallowed = false; bistable_storage *S = cov->Sbistable; assert(S != NULL); if (PisNULL(BIStablescale)) { PALLOC(BIStablescale, 3, 1); scale = P(BIStablescale); for (int i = 0; i < 3; scale[i++] = 1.0); } // set a = 1/s a[i11] = 1.0/scale[i11]; a[i21] = 1.0/scale[i21]; a[i22] = 1.0/scale[i22]; // default variance is 1 if (PisNULL(BIStablecdiag)) { PALLOC(BIStablecdiag, 2, 1); cdiag = P(BIStablecdiag); cdiag[0] = 1; cdiag[1] = 1; } // alpha diagonal are set => we are in likelohood estimation // we map alphadiag and beta red to alpha[i11], alpha[i21], alpha[i22] // alpha[i21] = betared*betaa + betac // beta = betaa + (2 - betaa)*beta_red if (S->alphadiag_given) { alphadiag = P(BIStablealphadiag); betaa = MAX(alphadiag[0], alphadiag[1]); betac = 2 - betaa; if (!PisNULL(BIStablebetared)) { betared = P0(BIStablebetared); } if (s->check && !PisNULL(BIStablealpha) ) { // if both alphas and alphadiag are given and check is true, // then check if they alphas and alphadiag are equal (or close to equal) alpha = P(BIStablealpha); if (cov->nrow[BIStablealpha] != 3 || cov->ncol[BIStablealpha] != 1) QERRC(BIStablealpha, "must be a 3 x 1 vector"); if (FABS(alpha[i11] - alphadiag[0]) > alpha[i11] * epsilon || FABS(alpha[i22] - alphadiag[1]) > alpha[i22] * epsilon || FABS(alpha[i21] - (betaa + betac*betared) ) > alpha[i21] * epsilon) { // PMI(cov); //printf("\n\n\n --- FABS(alpha[i21] - (betared*betaa + betac) ) = %10g \n", FABS(alpha[i21] - (betaa + betared*betac) ) ); // printf("\n\n\n --- alpha[i21] * epsilon = %10g \n", alpha[i21] * epsilon ); QERRC2(BIStablealpha, "does not match '%.50s' and '%.50s'.", BIStablealphadiag, BIStablebetared); } } // allocate memory for the alphas if (PisNULL(BIStablealpha)) { PALLOC(BIStablealpha, 3, 1); } alpha = P(BIStablealpha); // fill the alphas alpha[i11] = alphadiag[0]; alpha[i22] = alphadiag[1]; alpha[i21] = betaa + betared*betac; // These combinations of smoothness parameters are not allowed // if user mode, return an Error and stop // if likelihood mode, set correlation to 0 //notallowed = ( alpha[i21] < MAX(alpha[i11], alpha[i22]) ) || } else { //alphas are given if ( !PisNULL(BIStablealpha) ) { // all alphas are given, s->check => parameters are set by user, // not likelihood, or they are NA // we map alphas to alphadiag and betared if (PisNULL(BIStablealphadiag)) PALLOC(BIStablealphadiag, 2, 1); alphadiag = P(BIStablealphadiag); if (PisNULL(BIStablebetared)) PALLOC(BIStablebetared, 1, 1); alphadiag[0] = alpha[i11]; // alpha[i21]; alphadiag[1] = alpha[i22]; betaa = MAX(alphadiag[0], alphadiag[1]); betac = 2 - MAX(alphadiag[0], alphadiag[1]); P(BIStablebetared)[0] = (alpha[i21] - betaa)/betac ; betared = P(BIStablebetared)[0]; } if (PisNULL(BIStablealpha)) { // PMI(cov); QERRC2(BIStablealpha, "either '%.50s' or '%.50s' must be set", BIStablealpha, BIStablealphadiag); } } notallowed = ( ( alpha[i11] == alpha[i21] ) && ( alpha[i22] == alpha[i21] ) && ( POW(a[i21], alpha[i11]) < 0.5*POW(a[i11], alpha[i11])+0.5* POW(a[i22], alpha[i11]) ) ) || ( ( alpha[i11] == alpha[i21] ) && ( alpha[i11] > alpha[i22] ) && ( a[i21] <= POW(0.5, 1/alpha[i11])*a[i11] ) ) || ( ( alpha[i22] == alpha[i21] ) && ( alpha[i22] > alpha[i11] ) && ( a[i21] <= POW(0.5, 1/alpha[i22])*a[i22] ) )|| ( alpha[i11] > alpha[i21] ) || ( alpha[i22] > alpha[i21] ) ; if( notallowed && !s->check ) { {rhomax = 0;} } if( notallowed && s->check ) { QERRC(BIStablealpha, "This combination of smoothness parameters is not allowed."); } // alphas are set // find maximum allowable rho if (rhomax != 0) { biStableInterval(alpha, a, dim, &left, &right ); if ( (right == 0) && (left == 0) ) { rhomax = 0; } else { // PMI(cov); // printf("\nDear Olga, please check the PMI output; if the output looks reasonable you should check your code in the next line.\n"); biStableMinRho(cov, alpha, a, left, right, &rhomax); // printf("done\n"); } } // if rho is given, then it is set by user. Check if it is correct or not // ERR("Olga, there must be a call of type \"if (rhored_given)...\""); // if (!PisNULL(BIStablerhored) ){ if ( S->rhored_given ){ //if rhored is given, we are in the likelihood estimation if (PisNULL(BIStablerho)) PALLOC(BIStablerho, 1, 1); rho = P(BIStablerho); *rho = P(BIStablerho)[0] = (*rhored)*rhomax; } else if (!PisNULL(BIStablerho) && s->check ) { if ( PisNULL(BIStablerhored) ) PALLOC(BIStablerhored, 1, 1); rhored = P(BIStablerhored); *rhored = (*rho)/rhomax; } else if (!PisNULL(BIStablerho) && s->check ) { if (FABS(*rho) > rhomax ) { QERRC(BIStablealpha, "The value of cross-correlation parameter rho is too high."); } if ( PisNULL(BIStablerhored) ) PALLOC(BIStablerhored, 1, 1); rhored = P(BIStablerhored); *rhored = P(BIStablerhored)[0]= (*rho)/rhomax; } // printf("end init: alphadiag_given = %d\n", S->alphadiag_given); cov->initialised = true; RETURN_NOERROR; } void coinitbiStable(model *cov, localinfotype *li) { double thres = 1, *alpha = P(BIStablealpha); if ( ( alpha[0] <= thres ) && ( alpha[1] <= thres ) && ( alpha[2] <= thres ) ) { li->instances = 1; li->value[0] = 1; // q[CUTOFF_A] li->msg[0] =MSGLOCAL_OK; } else { li->instances = 1; li->value[0] = CUTOFF_THIRD_CONDITION ; // q[CUTOFF_A] li->msg[0] = MSGLOCAL_JUSTTRY; } } void biStable (double *x, model *cov, double *v) { int i; double alpha = P0(STABLE_ALPHA), *cdiag = P(BIStablecdiag), y = *x, z; assert(BIStablealpha == STABLE_ALPHA); /* biStable_storage *S = cov->SbiStable; assert(S != NULL); */ for (i=0; i<3; i++) { z = y/P(BIStablescale)[i]; P(STABLE_ALPHA)[0] = P(BIStablealpha)[i]; stable(&z, cov, v + i); } P(BIStablealpha)[0] = alpha; v[0] = v[0]*cdiag[0]; v[3] = v[2]*cdiag[1]; v[1] *= P0(BIStablerho)*SQRT(cdiag[0]*cdiag[1]); v[2] = v[1]; // printf("(%4.3f, %4.3f; %4.3e %4.3e %4.3e %4.3e)\t", x[0], x[1], v[0], v[1], v[2], v[3]); // PMI(cov->calling->calling); // if (!R_FINITE(v[0]) || !R_FINITE(v[1]) || !R_FINITE(v[2]) || !R_FINITE(v[3])) { PMI(cov); printf("(%4.3f, %4.3f; %4.3e %4.3e %4.3e %4.3e)\t", x[0], x[1], v[0], v[1], v[2], v[3]); BUG; } } void DbiStable(double *x, model *cov, double *v) { int i; double alpha = P0(STABLE_ALPHA), *cdiag = P(BIStablecdiag), y = *x, z; assert(BIStablealpha == STABLE_ALPHA); /* biStable_storage *S = cov->SbiStable; assert(S != NULL); */ for (i=0; i<3; i++) { z = y/P(BIStablescale)[i]; P(STABLE_ALPHA)[0] = P(BIStablealpha)[i]; Dstable(&z, cov, v + i); v[i] /= P(BIStablescale)[i]; } P(BIStablealpha)[0] = alpha; v[0] = v[0]*cdiag[0]; v[3] = v[2]*cdiag[1]; v[1] *= P0(BIStablerho)*SQRT(cdiag[0]*cdiag[1]); v[2] = v[1]; } void DDbiStable(double *x, model *cov, double *v) { int i; double alpha = P0(STABLE_ALPHA), *cdiag = P(BIStablecdiag), y = *x, z; assert(BIStablealpha == STABLE_ALPHA); /* biStable_storage *S = cov->SbiStable; assert(S != NULL); */ for (i=0; i<3; i++) { z = y/P(BIStablescale)[i]; P(STABLE_ALPHA)[0] = P(BIStablealpha)[i]; DDstable(&z, cov, v + i); v[i] /= P(BIStablescale)[i]*P(BIStablescale)[i]; } P(BIStablealpha)[0] = alpha; v[0] = v[0]*cdiag[0]; v[3] = v[2]*cdiag[1]; v[1] *= P0(BIStablerho)*SQRT(cdiag[0]*cdiag[1]); v[2] = v[1]; } void D3biStable(double *x, model *cov, double *v) { int i; double alpha = P0(STABLE_ALPHA), *cdiag = P(BIStablecdiag), y = *x, z; assert(BIStablealpha == STABLE_ALPHA); /* biStable_storage *S = cov->SbiStable; assert(S != NULL); */ for (i=0; i<3; i++) { z = y/P(BIStablescale)[i]; P(STABLE_ALPHA)[0] = P(BIStablealpha)[i]; D3stable(&z, cov, v + i); v[i] /= P(BIStablescale)[i]*P(BIStablescale)[i]*P(BIStablescale)[i]; } P(BIStablealpha)[0] = alpha; v[0] = v[0]*cdiag[0]; v[3] = v[2]*cdiag[1]; v[1] *= P0(BIStablerho)*SQRT(cdiag[0]*cdiag[1]); v[2] = v[1]; } void D4biStable(double *x, model *cov, double *v) { int i; double alpha = P0(STABLE_ALPHA), *cdiag = P(BIStablecdiag), y = *x, z; // assert(BIStablealpha == STABLE_ALPHA); /* biStable_storage *S = cov->SbiStable; assert(S != NULL); */ // assert(cov->initialised); for (i=0; i<3; i++) { z = y/P(BIStablescale)[i]; P(STABLE_ALPHA)[0] = P(BIStablealpha)[i]; D4stable(&z, cov, v + i); double dummy =P(BIStablescale)[i]*P(BIStablescale)[i]; v[i] /= dummy*dummy; } P(BIStablealpha)[0] = alpha; v[0] = v[0]*cdiag[0]; v[3] = v[2]*cdiag[1]; v[1] *= P0(BIStablerho)*SQRT(cdiag[0]*cdiag[1]); v[2] = v[1]; }