https://github.com/cran/RandomFields
Tip revision: c72cbea0dd8d0cd806e6d383a207338bd91db917 authored by Martin Schlather on 05 August 2005, 00:00:00 UTC
version 1.3.1
version 1.3.1
Tip revision: c72cbea
RFcircembed.cc
/*
Authors
Yindeng, Jiang, jiangyindeng@gmail.com
Martin Schlather, schlath@hsu-hh.de
Simulation of a random field by circulant embedding
(see Wood and Chan, or Dietrich and Newsam for the theory )
Copyright (C) 2001 -- 2003 Martin Schlather
Copyright (C) 2004 -- 2005 Yindeng Jiang & Martin Schlather
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR 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 <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "RFsimu.h"
#include <assert.h>
#include <R_ext/Applic.h>
ce_param CIRCEMBED={false, true, TRIVIALSTRATEGY, -1e-7, 1e-3, 3, 20000000,
0, 0, 0, 0};
ce_param LOCAL_CE={false, true, TRIVIALSTRATEGY, -1e-9, 1e-7, 1, 20000000,
0, 0, 0, 0};
void FFT_destruct(FFT_storage *FFT)
{
if (FFT->iwork!=NULL) {free(FFT->iwork); FFT->iwork=NULL;}
if (FFT->work!=NULL) {free(FFT->work); FFT->work=NULL;} //?
}
void FFT_NULL(FFT_storage *FFT)
{
FFT->work = NULL;
FFT->iwork = NULL;
}
void CE_destruct(void **S)
{
if (*S!=NULL) {
CE_storage *x;
x = *((CE_storage**)S);
if (x->c!=NULL) free(x->c);
if (x->d!=NULL) free(x->d);
FFT_destruct(&(x->FFT));
free(*S);
*S = NULL;
}
}
void localCE_destruct(void **S)
{
if (*S!=NULL) {
localCE_storage* x;
x = *((localCE_storage**) S);
DeleteKeyNotTrend(&(x->key));
free(*S);
*S = NULL;
}
}
/*********************************************************************/
/* CIRCULANT EMBEDDING METHOD (1994) ALGORITHM */
/* (it will always be refered to the paper of Wood & Chan 1994) */
/*********************************************************************/
void SetParamCircEmbed( int *action, int *force, double *tolRe, double *tolIm,
int *trials, double *mmin, int *useprimes, int *strategy,
double *maxmem)
{
SetParamCE(action, force, tolRe, tolIm, trials, mmin, useprimes, strategy,
maxmem, &CIRCEMBED, "CIRCEMBED");
}
void SetParamLocal( int *action, int *force, double *tolRe, double *tolIm,
int *trials, double *mmin, int *useprimes, int *strategy,
double *maxmem)
{
SetParamCE(action, force, tolRe, tolIm, trials, mmin, useprimes, strategy,
maxmem, &LOCAL_CE, "LOCAL");
}
int fastfourier(double *data, int *m, int dim, bool first, bool inverse,
FFT_storage *FFT)
/* this function is taken from the fft function by Robert Gentleman
and Ross Ihaka, in R */
{
long int inv, nseg, n,nspn,i;
int maxf, maxp, Xerror;
if (first) {
int maxmaxf,maxmaxp;
nseg = maxmaxf = maxmaxp = 1;
/* do whole loop just for Xerror checking and maxmax[fp] .. */
for (i = 0; i<dim; i++) {
if (m[i] > 1) {
fft_factor(m[i], &maxf, &maxp);
if (maxf == 0) {Xerror=ERRORFOURIER; goto ErrorHandling;}
if (maxf > maxmaxf) maxmaxf = maxf;
if (maxp > maxmaxp) maxmaxp = maxp;
nseg *= m[i];
}
}
assert(FFT->work==NULL);
if ((FFT->work = (double*) malloc(4 * maxmaxf * sizeof(double)))==NULL) {
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
assert(FFT->iwork==NULL);
if ((FFT->iwork = (int*) malloc(maxmaxp * sizeof(int)))==NULL) {
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
FFT->nseg = nseg;
// nseg = LENGTH(z); see loop above
}
inv = (inverse) ? 2 : -2;
n = 1;
nspn = 1;
nseg = FFT->nseg;
for (i = 0; i < dim; i++) {
if (m[i] > 1) {
nspn *= n;
n = m[i];
nseg /= n;
fft_factor(n, &maxf, &maxp);
fft_work(&(data[0]), &(data[1]), nseg, n, nspn, inv, FFT->work,FFT->iwork);
}
}
return NOERROR;
ErrorHandling:
FFT_destruct(FFT);
return Xerror;
}
int fastfourier(double *data, int *m, int dim, bool first, FFT_storage *FFT){
return fastfourier(data, m, dim, first, !first, FFT);
}
int init_circ_embed(key_type *key, int m)
{
methodvalue_type *meth;
int Xerror=NOERROR, d, actcov;
long twoRealmtot;
double steps[MAXDIM], *c;
CE_storage *s;
ce_param* cepar;
bool // Critical[MAXDIM],
critical;
c = NULL;
meth = &(key->meth[m]);
cepar = &CIRCEMBED;
if (!key->grid) {Xerror=ERRORMETHODNOTALLOWED;goto ErrorHandling;}
SET_DESTRUCT(CE_destruct, m);
if ((meth->S = malloc(sizeof(CE_storage)))==0){
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = (CE_storage*) meth->S;
s->c =NULL;
s->d =NULL;
FFT_NULL(&(s->FFT));
if ((Xerror = FirstCheck_Cov(key, m, true)) != NOERROR) goto ErrorHandling;
actcov = meth->actcov;
critical = false;
// for (k=0; k<dim; k++) Critical[k]=false; // critical odd not used yet
// if (cov->type==ANISOTROPIC) { // if Hypermodels are
// // included it becomes much more difficult!
// if (!cov->even) {
// for(k=0; k<dim; k++) {
// Critical[k] |= cov->odd[k];
// }
// }
// }
for (d=0; d<key->timespacedim; d++) {
s->nn[d]=key->length[d];
steps[d]=key->x[d][XSTEP];
}
int *mm, *cumm, *halfm, dim;
double hx[MAXDIM], totalm;
int trials, index[MAXDIM], dummy;
long mtot,i,k,twoi;
bool positivedefinite, cur_crit;
mtot=-1;
mm = s->m;
halfm= s->halfm;
cumm = s->cumm;
dim = key->timespacedim;
c=NULL;
if (GENERAL_PRINTLEVEL>=5) PRINTF("calculating the Fourier transform\n");
/* cumm[i+1]=\prod_{j=0}^i m[j]
cumm is used for fast transforming the matrix indices into an
index of the vector (the way the matrix is stored) corresponding
to the matrix */
/* calculate the dimensions of the matrix C, eq. (2.2) in W&C */
// ("CE missing strategy for matrix entension in case of anisotropic fields %d\n
totalm = 1.0;
for (i=0;i<dim;i++){ // size of matrix at the beginning
mm[i] = s->nn[i];
if (cepar->mmin[i]>0.0) {
if (mm[i] > ((int) ceil(cepar->mmin[i])) / 2) {
sprintf(ERRORSTRING_OK, "Minimum size in direction %d is %d",
(int) i, mm[i]);
sprintf(ERRORSTRING_WRONG,"%d", (int) ceil(cepar->mmin[i]));
return ERRORMSG;
}
mm[i] = ((int) ceil(cepar->mmin[i])) / 2;
} else if (cepar->mmin[i] < 0.0) {
assert(cepar->mmin[i] <= -1.0);
mm[i] = (int) ceil((double) mm[i] * -cepar->mmin[i]);
}
if (cepar->useprimes) {
// Note! algorithm fails if mm[i] is not a multiple of 2
// 2 may not be put into NiceFFTNumber since it does not
// guarantee that the result is even even if the input is even !
mm[i] = 2 * NiceFFTNumber((unsigned long) mm[i]);
} else {
mm[i] = (1 << 1 + (int) ceil(log((double) mm[i]) * INVLOG2 - EPSILON1000));
}
totalm *= (double) mm[i];
}
if (totalm > cepar->maxmem) {
sprintf(ERRORSTRING_OK, "%f", cepar->maxmem);
sprintf(ERRORSTRING_WRONG,"%f", totalm);
return ERRORMAXMEMORY;
}
positivedefinite = false;
/* Eq. (3.12) shows that only j\in I(m) [cf. (3.2)] is needed,
so only the first two rows of (3.9) (without the taking the
modulus of h in the first row)
The following variable `index' corresponds to h(l) in the following
way: index[l]=h[l] if 0<=h[l]<=mm[l]/2
index[l]=h[l]-mm[l] if mm[l]/2+1<=h[l]<=mm[l]-1
Then h[l]=(index[l]+mm[l]) mod mm[l] !!
*/
/* The algorithm below:
while (!positivedefinite && (trials<cepar->trials)){
trials++;
calculate the covariance values "c" according to the given "m"
fastfourier(c)
if (!cepar->force || (trials<cepar->trials)) {
check if positive definite
if (!positivedefinite && (trials<cepar->trials)) {
enlarge "m"
}
} else
print "forced"
}
*/
trials=0;
while (!positivedefinite && (trials<cepar->trials)){
trials++;
cumm[0]=1;
for(i=0;i<dim;i++){
halfm[i]=mm[i]/2;
index[i]=1-halfm[i];
cumm[i+1]=cumm[i] * mm[i]; // only relevant up to i<dim-1 !!
}
mtot=cumm[dim-1] * mm[dim-1];
if (GENERAL_PRINTLEVEL>=2) {
for (i=0;i<dim;i++) PRINTF("mm[%d]=%d, ",i,mm[i]);
PRINTF("mtot=%d\n ",mtot);
}
twoRealmtot = 2 * sizeof(double) * mtot;
// for the following, see the paper by Wood and Chan!
// meaning of following variable c, see eq. (3.8)
if ((c = (double*) malloc(twoRealmtot)) == 0) {
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
for(i=0; i<dim; i++){index[i]=0;}
for (i=0; i<mtot; i++){
cur_crit = false;
for (k=0; k<dim; k++) {
hx[k] = steps[k] *
(double) ((index[k] <= halfm[k]) ? index[k] : index[k] - mm[k]);
cur_crit |= (index[k]==halfm[k]);
}
dummy=i << 1;
c[dummy] = (critical && cur_crit) ? 0.0 :
key->covFct(hx, dim, key->cov, meth->covlist, actcov, key->anisotropy);
// statt critical hier besser andere CovFct addieren!
c[dummy+1] = 0.0;
k=0; while( (k<dim) && (++(index[k]) >= mm[k])) {
index[k]=0;
k++;
}
assert( (k<dim) || (i==mtot-1));
}
if (GENERAL_PRINTLEVEL>6) PRINTF("FFT...");
if ((Xerror=fastfourier(c, mm, dim, true, &(s->FFT)))!=0)
goto ErrorHandling;
if (GENERAL_PRINTLEVEL>6) PRINTF("finished\n");
// check if positive definite. If not: enlarge and restart
if (!cepar->force || (trials<cepar->trials)) {
long int mtot2;
mtot2 = mtot * 2;
twoi=0;
// 16.9. < cepar.tol.im changed to <=
while ((twoi<mtot2) && (positivedefinite=(c[twoi]>=cepar->tol_re) &&
(fabs(c[twoi+1])<=cepar->tol_im)))
{twoi+=2;}
if ( !positivedefinite) {
if (GENERAL_PRINTLEVEL>=2)
// 1.1.71: %f changed to %e because c[twoi+1] is usually very small
PRINTF("non-positive eigenvalue: c[%d])=%e + %e i.\n",
(int) (twoi/2), c[twoi], c[twoi+1]);
if (GENERAL_PRINTLEVEL>=4) { // just for printing the smallest
// eigenvalue (min(c))
long int index=twoi, sum=0;
double smallest=c[twoi];
char percent[]="%";
for (twoi=0; twoi<mtot2; twoi += 2) {
if (c[twoi] < 0) {
sum++;
if (c[twoi]<smallest) { index=twoi; smallest=c[index]; }
}
}
PRINTF(" The smallest eigenvalue is c[%d] = %e.\n",
index / 2, smallest);
PRINTF(" The percentage of negative eigenvalues is %5.3f%s.\n",
100.0 * (double) sum / (double) mtot, percent);
}
}
if (!positivedefinite && (trials<cepar->trials)) {
FFT_destruct(&(s->FFT));
free(c); c=NULL;
totalm = 1.0;
switch (cepar->strategy) {
case 0 :
for (i=0;i<dim;i++) { /* enlarge uniformly in each direction, maybe
this should be modified if the grid has
different sizes in different directions */
mm[i] <<= 1;
totalm *= (double) mm[i];
}
break;
case 1 :
double cc, maxcc, hx[MAXDIM];
int maxi;
maxi = -1;
maxcc = RF_NEGINF;
for (i=0; i<dim; i++) hx[i] = 0.0;
for (i=0; i<dim; i++) {
hx[i] = steps[i] * mm[i];
cc = fabs(key->covFct(hx, dim, key->cov, meth->covlist, actcov,
key->anisotropy));
if (GENERAL_PRINTLEVEL>2) PRINTF("%d cc=%e (%e)",i,cc,hx[i]);
if (cc>maxcc) {
maxcc = cc;
maxi = i;
}
hx[i] = 0.0;
}
assert(maxi>=0);
mm[maxi] <<= 1;
for (i=0;i<dim; i++) totalm *= (double) mm[i];
break;
default:
assert(false);
}
if (totalm>cepar->maxmem) {
sprintf(ERRORSTRING_OK, "%f", cepar->maxmem);
sprintf(ERRORSTRING_WRONG,"%f", totalm);
Xerror=ERRORMAXMEMORY;
goto ErrorHandling;
}
// assert(false);
}
} else {if (GENERAL_PRINTLEVEL>=2) PRINTF("forced\n");}
}
assert(mtot>0);
if (positivedefinite || cepar->force) {
// correct theoretically impossible values, that are still within
// tolerance CIRCEMBED.tol_re/CIRCEMBED.tol_im
double r, imag;
r = imag = 0.0;
for(i=0,twoi=0;i<mtot;i++) {
if (c[twoi] > 0.0) {
c[twoi] = sqrt(c[twoi]);
} else {
if (c[twoi] < r) r = c[twoi];
c[twoi] = 0.0;
}
{
register double a;
if ((a=fabs(c[twoi+1])) > imag) imag = a;
}
c[twoi+1] = 0.0;
twoi+=2;
}
if (GENERAL_PRINTLEVEL>1) {
if (r < -GENERAL_PRECISION || imag > GENERAL_PRECISION) {
PRINTF("using approximating circulant embedding:\n");
if (r < -GENERAL_PRECISION)
PRINTF("\tsmallest real part has been %e \n", r);
if (imag > GENERAL_PRECISION)
PRINTF("\tlargest modulus of the imaginary part has been %e \n",imag);
}
}
} else {Xerror=ERRORFAILED;goto ErrorHandling;}
if (GENERAL_PRINTLEVEL>=10) {
for (i=0;i<2*mtot;i++) {PRINTF("%f ",c[i]);} PRINTF("\n");
}
// s->c = c;
// return NOERROR;
if (key->storing) {
if ((s->d=(double *) malloc(twoRealmtot))==0){
Xerror=ERRORMEMORYALLOCATION;goto ErrorHandling;} //d
}
s->c = c;
return NOERROR;
ErrorHandling:
if (c!=NULL) {free(c);}
return Xerror;
}
void covcpy(covinfo_type *dest, covinfo_type *source){
assert(source->x==NULL);
memcpy(dest, source, sizeof(covinfo_type));
}
double GetScaledDiameter(key_type *key, covinfo_type *kc)
{
// SCALE and ANISO is considered as space trafo and envolved here
double diameter;
diameter = 0.0;
if (key->anisotropy) { // code can be reduced to the anisotropic case
// with the obvious loss of time in the isotropic case...
int i, j, k, dim, ncorner_dim;
double distsq;
double sx[ZWEIHOCHMAXDIM * MAXDIM];
dim = kc->truetimespacedim;
ncorner_dim = (1 << key->timespacedim) * dim;
GetCornersOfGrid(key, dim, kc->aniso, sx);
for (i=1; i<ncorner_dim; i+=dim) {
for (j=0; j<i; j+=dim) {
distsq = 0.0;
for (k=0; k<dim; k++) {
register double dummy;
dummy = sx[i + k] - sx[j + k];
distsq += dummy * dummy;
}
if (distsq > diameter) diameter = distsq;
}
}
} else { // see above
int d;
for (d=0; d<key->timespacedim; d++) {
register double dummy;
dummy = key->x[d][XSTEP] * (double) (key->length[d] - 1);
diameter += dummy * dummy;
}
diameter *= kc->aniso[0] * kc->aniso[0];
}
// printf("diameter=%f", sqrt(diameter));
return sqrt(diameter);
}
int init_circ_embed_local(key_type *key, int m)
{
methodvalue_type *meth;
int Xerror=NOERROR, v, store_msg[MAXDIM], msg, simuactcov, instance,
actcov, nc, store_nr;
localCE_storage *s;
cov_fct *hyper;
covinfo_type *sc;
bool selectlocal[Forbidden + 1];
param_type store_param;
double rawRmax;
if (key->covFct != CovFct) { Xerror=ERRORNOTPROGRAMMED; goto ErrorHandling;}
SET_DESTRUCT(localCE_destruct, m);
if (!key->grid) { Xerror=ERRORMETHODNOTALLOWED; goto ErrorHandling;}
meth = &(key->meth[m]);
assert(meth->S==NULL);
if ((meth->S=malloc(sizeof(localCE_storage)))==0){
Xerror=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
s = (localCE_storage*) meth->S;
KEY_NULL(&(s->key));
// prepared for more sophisticated selections if init_circ_embed_local
// and init_circ_embed are merged
for (v=0; v < (int) Forbidden; v++) selectlocal[v]=false;
selectlocal[meth->unimeth] = true;
if ((Xerror = FirstCheck_Cov(key, m, false)) != NOERROR) goto ErrorHandling;
actcov = meth->actcov;
if (2 * actcov > MAXCOV) {Xerror=ERRORNCOVOUTOFRANGE; goto ErrorHandling;}
rawRmax = 0.0;
for (simuactcov=v=0; v<actcov; v++) {
covinfo_type *kc;
store_nr = -1;
kc = &(key->cov[meth->covlist[v]]);
// printf("** %d\n", kc->nr);
sc = &(s->key.cov[simuactcov]);
sc->method = CircEmbed;
sc->param[VARIANCE] = 1.0;
sc->param[HYPERNR] = 1; // as only addition between models allowed
sc->param[DIAMETER] = GetScaledDiameter(key, kc); // SCALE is considered
// as space trafo and envolved here
if (GENERAL_PRINTLEVEL>7) PRINTF("diameter %f\n", sc->param[DIAMETER]);
sc->op = 2;
covcpy(sc + 1, kc);
instance = -1;
store_param[LOCAL_R] = R_PosInf;
store_msg[simuactcov] = MSGLOCAL_JUSTTRY; // lowest preference
for (nc=0; nc < nlocal; nc++) {// in case more cutoff methods are found
instance = 0;
sc->nr = LocalCovList[nc];
hyper = &(CovList[sc->nr]);
if (GENERAL_PRINTLEVEL>7) PRINTF("%d %s:\n", nc, hyper->name);
if (!selectlocal[hyper->localtype]) continue;
// for future use if more hyper functions for cutoff and intrinsic
// are available
// if (PREFERENCE[hyper->localtype] >= 0 &&
// PREFERENCE[hyper->localtype]!=sc->nr) continue;
while ((msg = hyper->getparam(sc + 1, sc->param, instance++))
!= MSGLOCAL_ENDOFLIST) {
// printf("%d %s %d \n", instance, hyper->name, msg);
if ((Xerror=
hyper->checkNinit(sc, COVLISTALL, key->ncov - meth->covlist[v],
CircEmbed)) == NOERROR
&& (
(int) msg < (int) store_msg[simuactcov]
||
((int) msg == (int) store_msg[simuactcov] &&
sc->param[LOCAL_R] < store_param[LOCAL_R])
)
)
{
memcpy(store_param, sc->param, sizeof(param_type));
store_msg[simuactcov] = msg;
store_nr = sc->nr;
} // if
// printf("%d %d %f \n", Xerror, store_msg[simuactcov], store_param[LOCAL_R] );
if (GENERAL_PRINTLEVEL>7) {
PRINTF("v=%d nc=%d inst=%d err=%d hyp.kappa=%f, #=%d, diam=%f\n r=%f curmin_r=%f\n",
v, nc, instance, Xerror, sc->param[HYPERKAPPAI],
(int) sc->param[HYPERNR],
sc->param[DIAMETER], sc->param[LOCAL_R],
store_param[LOCAL_R]);
if (hyper->implemented[CircEmbedIntrinsic] == HYPERIMPLEMENTED)
PRINTF("userr=%f, a0=%f, a2=%f, b=%f\n",
sc->param[INTRINSIC_RAWR],
sc->param[INTRINSIC_A0], sc->param[INTRINSIC_A2],
sc->param[INTRINSIC_B]
);
else {
// printf("%s %d\n",
// hyper->name, hyper->implemented[CircEmbedCutoff]);
assert(hyper->implemented[CircEmbedCutoff] == HYPERIMPLEMENTED);
PRINTF("a=%f, sqrt.a=%f, b=%f\n",
sc->param[CUTOFF_A],
sc->param[CUTOFF_ASQRTR], sc->param[CUTOFF_B]
);
}
}
} // while
if (rawRmax < sc->param[LOCAL_R] / sc->param[DIAMETER])
rawRmax = sc->param[LOCAL_R] /sc->param[DIAMETER];
} // nc
if (!R_FINITE(store_param[LOCAL_R])) {
PRINTF("v=%d nc=%d, inst=%d err=%d hyp.kappa=%f, #=%d, diam=%f\nr=%f curmin_r=%f\n",
v, nc, instance, Xerror, sc->param[HYPERKAPPAI],
(int) sc->param[HYPERNR],
sc->param[DIAMETER], sc->param[LOCAL_R],
store_param[LOCAL_R]);
assert(Xerror!=NOERROR);
goto ErrorHandling;
}
sc->nr = store_nr;
assert(sc->nr >= 0);
memcpy(sc->param, store_param, sizeof(param_type));
simuactcov += (short int) sc->param[HYPERNR] + 1;
sc += (int) sc->param[HYPERNR];
} // v
// prepare for call of internal_InitSimulateRF
int covnr[MAXCOV], op[MAXCOV], CEMethod[MAXCOV], cum_nParam,
aniso, n_aniso, i;
double ParamList[MAXCOV * TOTAL_PARAM];
char errorloc_save[nErrorLoc];
if (key->anisotropy) {
n_aniso = key->timespacedim * key->timespacedim ;
aniso = ANISO;
} else {
n_aniso = 1;
aniso = SCALE;
}
cum_nParam = 0;
for (v=0; v<simuactcov; v++) {
sc = &(s->key.cov[v]);
covnr[v] = sc->nr;
op[v] = sc->op;
CEMethod[v] = (int) CircEmbed;
ParamList[cum_nParam++] = sc->param[VARIANCE];
for (i=0; i<CovList[sc->nr].kappas; i++)
ParamList[cum_nParam++] = sc->param[KAPPA + i];
for (i=0; i<n_aniso; i++)
ParamList[cum_nParam++] = sc->param[aniso + i];
}
//for (i=0; i<cum_nParam; printf("%d %f \n", i, ParamList[i++]));
// assert(false);
strcpy(errorloc_save, ERROR_LOC);
sprintf(ERROR_LOC, "%s%s ", errorloc_save, "cutoff: ");
ce_param ce_save;
memcpy(&ce_save, &CIRCEMBED, sizeof(ce_param));
memcpy(&CIRCEMBED, &LOCAL_CE, sizeof(ce_param));
rawRmax *= sqrt(key->timespacedim);
for (i=0; i<key->timespacedim; i++) {
if (CIRCEMBED.mmin[i]==0.0) CIRCEMBED.mmin[i] = - rawRmax;
}
instance = 0;
for (;;) {
bool anychangings;
Xerror = internal_InitSimulateRF(key->x[0], key->T, key->spatialdim,
3, key->grid, key->Time,
covnr, ParamList, cum_nParam,
simuactcov, key->anisotropy, op,
CEMethod, DISTR_GAUSS, &(s->key),
GENERAL_STORING,
0 /* natural scaling */, CovFct);
anychangings = false;
if (Xerror == NOERROR || Xerror == USEOLDINIT) {
break;
} else {
for (v=0; v<simuactcov; v++, instance++) {
sc = &(s->key.cov[v]);
hyper = &(CovList[sc->nr]);
if ((hyper->localtype == CircEmbedCutoff ||
hyper->localtype == CircEmbedIntrinsic)) {
if (store_msg[v] != MSGLOCAL_OK) {
bool improved;
memcpy(store_param, sc->param, sizeof(param_type));
improved = hyper->alternative(sc, instance) &&
(hyper->checkNinit(sc, COVLISTALL, (int) sc->param[HYPERNR],
CircEmbed) != NOERROR);
if (improved) memcpy(sc->param, store_param, sizeof(param_type));
anychangings |= improved;
}
} else {
assert((v % 2) == 1); // if the method is generalised to
// multiplicative submodels this assert reminds that
// this part must be changed.
}
}
if (!anychangings) break;
}
}
memcpy(&CIRCEMBED, &ce_save, sizeof(ce_param));
strcpy(ERROR_LOC, errorloc_save);
if (Xerror != NOERROR) goto ErrorHandling;
s->factor = 0.0;
for (v=0; v<simuactcov; v++) {
sc = &(s->key.cov[v]);
hyper = &(CovList[sc->nr]);
if (hyper->localtype == CircEmbedIntrinsic) {
static double zero=0.0;
// these steps are necessary, see the two cases in 3dBrownian!
s->factor += sc->param[INTRINSIC_A2] *
key->covFct(&zero, 1, sc, COVLISTALL, (int) sc->param[HYPERNR],
false) / (sc->param[DIAMETER] * sc->param[DIAMETER]);
}
}
// 2.0 : see Stein (2002)
s->factor = sqrt(2.0 * s->factor); // standard deviation of the Gaussian
// variables in do_...
return NOERROR;
ErrorHandling:
return Xerror;
}
void do_circ_embed_local(key_type *key, int m, double *res )
{
double x[MAXDIM], dx[MAXDIM];
long index[MAXDIM], k, r;
localCE_storage *s;
assert(key->active);
s = (localCE_storage*) key->meth[m].S;
assert(s->factor >= 0.0);
internal_DoSimulateRF(&(s->key), 1, res);
if (s->factor>0) {
for (k=0; k<key->timespacedim; k++) {
index[k]=0;
dx[k]= GAUSS_RANDOM(s->factor * key->x[k][XSTEP]);
x[k]= 0.0;
}
for(r=0;;) {
for (k=0; k<key->timespacedim; k++) res[r] += x[k];
r++;
k=0;
while( (k<key->timespacedim) && (++index[k]>=key->length[k])) {
index[k]=0;
x[k] = 0.0;
k++;
}
if (k>=key->timespacedim) break;
x[k] += dx[k];
}
}
}
void do_circ_embed(key_type *key, int m, double *res ){
double *d;
CE_storage *s;
s = (CE_storage*)key->meth[m].S;
assert(key->active);
if (s->d==NULL) { /* overwrite the intermediate
result directly (algorithm allows for that) */
d=s->c;
assert(!key->storing);
}
else {
assert(key->storing);
d=s->d;
}
/*
implemented here only for rotationsinvariant covariance functions
for arbitrary dimensions;
(so it works only for even covariance functions in the sense of
Wood and Chan,p. 415, although they have suggested a more general
algorithm;)
Warning! If key->storing==false when calling init_circ_embed
and key->storing==true when calling do_circ_embed, the
complete programme will fail, since the initialization depends on
the value of key->storing
*/
int i, j, k, HalfMp1[MAXDIM], HalfMaM[2][MAXDIM], index[MAXDIM], dim,
*nn, *mm, *cumm, *halfm;
double XX,YY,invsqrtmtot, *c;
bool first, free[MAXDIM+1], noexception;
long mtot;
dim = key->timespacedim;
nn = s->nn;
mm = s->m;
cumm = s->cumm;
halfm = s->halfm;
c = s->c;
mtot=cumm[dim-1] * mm[dim-1];
for (i=0; i<dim; i++) {
HalfMp1[i] = ((mm[i] % 2)==1) ? -1 : halfm[i];
HalfMaM[0][i] = halfm[i];
HalfMaM[1][i] = mm[i] - 1;
}
//if there are doubts that the algorithm below reaches all elements
//of the matrix, uncomment the lines containing the variable xx
//
//bool *xx; xx = (bool*) malloc(sizeof(bool) * mtot);
//for (i=0; i<mtot;i++) xx[i]=true;
invsqrtmtot = 1/sqrt((double) mtot);
if (GENERAL_PRINTLEVEL>=10) PRINTF("Creating Gaussian variables... \n");
/* now the Gaussian r.v. have to defined and multiplied with sqrt(FFT(c))*/
for (i=0; i<dim; i++) {
index[i]=0;
free[i] = false;
}
free[MAXDIM] = false;
// if MaxStableRF calls Cutofff, c and look uninitialised
for(;;) {
i = j = 0;
noexception = false;
for (k=0; k<dim; k++) {
i += cumm[k] * index[k];
if ((index[k]==0) || (index[k]==HalfMp1[k])) {
j += cumm[k] * index[k];
} else {
noexception = true; // not case 2 in prop 3 of W&C
j += cumm[k] * (mm[k] - index[k]);
}
}
//A 0 0 0 0 0 0.000000 0.000000
//B 0 0 0 0 0 -26.340334 0.000000
// if (index[1] ==0) printf("A %d %d %d %d %d %f %f\n",
// index[0], index[1], i, j, noexception, d[i], d[i+1]);
if (GENERAL_PRINTLEVEL>=10) PRINTF("cumm...");
i <<= 1; // since we have to index imaginary numbers
j <<= 1;
if (noexception) { // case 3 in prop 3 of W&C
XX = GAUSS_RANDOM(INVSQRTTWO);
YY = GAUSS_RANDOM(INVSQRTTWO);
d[i] = d[i+1] = c[i]; d[i] *= XX; d[i+1] *= YY;
d[j] = d[j+1] = c[j]; d[j] *= XX; d[j+1] *= -YY;
} else { // case 2 in prop 3 of W&C
d[i] = c[i] * GAUSS_RANDOM(1.0);
d[i+1] = 0.0;
}
// if (i==224|| i==226 || i==228) {
//printf("B %d %d %d %d %d %f %f\n",
// index[0], index[1], i, j, noexception, d[i], d[i+1]);
// assert(false);
// }
if (GENERAL_PRINTLEVEL>=10) PRINTF("k=%d ", k);
// if (index[1] == 0) printf("B %d %d %d %d %d %f %f\n",
// index[0], index[1], i, j, noexception, d[i], d[i+1]);
/*
this is the difficult part.
We have to run over roughly half the points, but we should
not run over variables twice (time lost)
Due to case 2, we must include halfm.
idea is:
for (i1=0 to halfm[dim-1])
if (i1==0) or (i1==halfm[dim-1]) then endfor2=halfm[dim-2]
else endfor2=mm[dim-2]
for (i2=0 to endfor2)
if ((i1==0) or (i1==halfm[dim-1])) and
((i2==0) or (i2==halfm[dim-2]))
then endfor3=halfm[dim-3]
else endfor3=mm[dim-3]
for (i3=0 to endfor3)
....
i.e. the first one that is not 0 or halfm (regarded from dim-1 to 0)
runs over 0..halfm, all the others over 0..m
this is realised in the following allowing for arbitrary value of dim
free==true <=> endfor==mm[]
*/
k=0;
if (++index[k]>HalfMaM[free[k]][k]) {
// in case k increases the number of indices that run over 0..m increases
free[k] = true;
index[k]= 0;
k++;
while((k<dim) && (++index[k]>HalfMaM[free[k]][k])) {
free[k] = true;
index[k]= 0;
k++;
}
if (k>=dim) break;
// except the very last (new) number is halfm and the next index is
// restricted to 0..halfm
// then k decreases as long as the index[k] is 0 or halfm
if (!free[k] && (index[k]==halfm[k])){//index restricted to 0..halfm?
// first: index[k] is halfm? (test on ==0 is superfluent)
k--;
while ( (k>=0) && ((index[k]==0) || (index[k]==halfm[k]))) {
// second and following: index[k] is 0 or halfm?
free[k] = false;
k--;
}
}
}
}
// printf("%f\n", d[0]);
// printf("%f\n", d[1]);
// printf("%f\n", d[2]);
// printf("%f\n", d[3]);
// double zz=0.0;
// for (i =0; i<2 * mtot; i++) {
// if (true || i<=100) printf("%d %d\n", i, 2 * mtot);
// zz += d[i];
// assert(i!=230);
// // if (i>100); break;
// }
// printf("%f\n", zz);
fastfourier(d, mm, dim, false, &(s->FFT));
/* now we correct the result of the fastfourier transformation
by the factor 1/sqrt(mtot) and read the relevant matrix out of
the large vector c */
first = true;
for(i=0; i<dim; i++){index[i]=0;}
int totpts = key->totalpoints;
for (i=0; i<totpts; i++){
j=0; for (k=0; k<dim; k++) {j+=cumm[k] * index[k];}
res[i] += d[2*j]*invsqrtmtot;
if ((fabs(d[2*j+1])>CIRCEMBED.tol_im) &&
((GENERAL_PRINTLEVEL>=2 && first) || GENERAL_PRINTLEVEL>=6)){
PRINTF("IMAGINARY PART <> 0, %e\n",d[2*j+1]);
first=false;
}
k=0; while((k<dim) && (++index[k]>=nn[k])) {index[k++]=0;}
}
}
