https://github.com/cran/RandomFields
Tip revision: 270c4ba0da8194f2f81416d85bae1ccbd3ed4933 authored by Martin Schlather on 08 August 1977, 00:00:00 UTC
version 1.0.4
version 1.0.4
Tip revision: 270c4ba
RFcircembed.cc
/*
Authors
Martin Schlather, Martin.Schlather@uni-bayreuth.de
Simulation of a random field by circulant embedding
(see Wood and Chan, or Dietrich and Newsam for the theory )
Copyright (C) 2001 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>
bool CIRCEMBED_FORCE=false; /* after CIRCEMBED_TRIALS, the result will be
accepted, although not within TOL */
int CIRCEMBED_TRIALS=3;/* how often grid is enlarged to hopefully get a
positive definite matrix */
int CIRCEMBED_MMIN = 0; /* cf TBMCE_MTOT, minimum M in each direction */
#ifdef doubleReal
Real CIRCEMBED_TOL_RE=-1e-7;/* <0; real part must be greater than TOL_RE to
be considered as zero */
Real CIRCEMBED_TOL_IM=1e-3;/* >0; maximum module of imaginary part that is
considered as zero */
#else
Real CIRCEMBED_TOL_RE=-1e-5;
Real CIRCEMBED_TOL_IM=1e-3;
#endif
typedef struct CE_storage {
int m[MAXDIM],halfm[MAXDIM],cumm[MAXDIM],nn[MAXDIM];
long totallength;
covfct cov;
Real param[TOTAL_PARAM];
double *c,*d,*local;
FFT_storage FFT;
} CE_storage;
void FFT_destruct(FFT_storage *FFT)
{
if (FFT->work!=NULL) {free(FFT->work); FFT->work=NULL;}
if (FFT->iwork!=NULL) {free(FFT->iwork); FFT->iwork=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);
if (x->local!=NULL) free(x->local);
FFT_destruct(&(x->FFT));
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, Real *tolRe, Real *tolIm,
int *trials, int *mmin)
{
switch(*action) {
case 0 :
CIRCEMBED_TRIALS=*trials;
if (CIRCEMBED_TRIALS<1) {
CIRCEMBED_TRIALS=1;
if (GENERAL_PRINTLEVEL>0) PRINTF("\nWARNING! CIRCEMBED_TRIALS had been less than 1\n");
}
CIRCEMBED_FORCE=(bool)*force;
CIRCEMBED_TOL_RE=*tolRe;
if (CIRCEMBED_TOL_RE>0) {
CIRCEMBED_TOL_RE=0;
if (GENERAL_PRINTLEVEL>0) PRINTF("\nWARNING! CIRCEMBED_TOL_RE had been positive.\n");
}
CIRCEMBED_TOL_IM=*tolIm;
if (CIRCEMBED_TOL_IM<0) {
CIRCEMBED_TOL_IM=0;
if (GENERAL_PRINTLEVEL>0) PRINTF("\nWARNING! CIRCEMBED_TOL_IM had been neagtive.\n");
}
CIRCEMBED_MMIN=*mmin;
if (CIRCEMBED_MMIN>0) {
CIRCEMBED_MMIN =
1 << (1+(int) (log(((double) CIRCEMBED_MMIN)-0.1)/log(2.0)));
if ((CIRCEMBED_MMIN!=*mmin) && (GENERAL_PRINTLEVEL>0))
PRINTF("\nWARNING! CIRCEMBED_MMIN set to %d.\n",CIRCEMBED_MMIN);
}
break;
case 1 :
*force=CIRCEMBED_FORCE;
*tolRe=CIRCEMBED_TOL_RE;
*tolIm=CIRCEMBED_TOL_IM;
*trials=CIRCEMBED_TRIALS;
*mmin=CIRCEMBED_MMIN;
if (GetNotPrint) break;
case 2 :
PRINTF("\nCIRC_EMBED\n==========\nforce=%d,\ntol_Re=%e\ntol_Im=%e\ntrials=%d\nmmin=%d\n",
CIRCEMBED_FORCE,CIRCEMBED_TOL_RE,CIRCEMBED_TOL_IM,
CIRCEMBED_TRIALS,CIRCEMBED_MMIN);
break;
default : PRINTF(" unknown action\n");
}
}
int fastfourier(double *data, int *m, int dim, bool first, FFT_storage *S)
#ifdef RF_GSL
// fftw_plan fftw_create_plan(int n, fftw_direction dir,int flags);
// fftwnd_plan fftwnd_create_plan(int rank, const int *n,
// fftw_direction dir, int flags);
// void fftwnd(fftwnd_plan plan, int howmany, fftw_complex *in,int istride,
// int idist, fftw_complex *out, int ostride, int odist);
{
assert(false);
}
#else
/* this function is taken from the fft function by Robert Gentleman
and Ross Ihaka, in R */
{
int inv, nseg, n,nspn,i,maxf,maxp,error;
if (first) {
int maxmaxf,maxmaxp;
inv = -2;
nseg = maxmaxf = maxmaxp = 1;
/* do whole loop just for error checking and maxmax[fp] .. */
for (i = 0; i <dim; i++) {
if (m[i] > 1) {
fft_factor(m[i], &maxf, &maxp);
if (maxf == 0) {error=ERRORFOURIER; goto ErrorHandling;}
if (maxf > maxmaxf) maxmaxf = maxf;
if (maxp > maxmaxp) maxmaxp = maxp;
nseg *= m[i];
}
}
if ((S->work = (double*)malloc(4 * maxmaxf * sizeof(double)))==NULL) {
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
if ((S->iwork = (int*)malloc(maxmaxp * sizeof(int)))==NULL) {
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
S->nseg = nseg;
// nseg = LENGTH(z); see loop above
} else {
inv = 2;
}
n = 1;
nspn = 1;
nseg = S->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, S->work, S->iwork);
}
}
return 0;
ErrorHandling:
FFT_destruct(S);
return error;
}
#endif
int internal_init_circulantembedding(key_type * key)
{
Real *c,h,invn2[MAXDIM],*param;
int *m, *halfm, *cumm, *nn, error,trials,index[MAXDIM],dim;
int maxm,MAXM[]={67108864,8192,256,64}; /* the maximum grid size:
MAXM[dim-1]^dim */
long mtot,i,j,k,twoi,twoRealmtot;
bool positivedefinite;
covfct cov;
c=NULL;
if(key->cov->check!=NULL) {
if(error=key->cov->check(key)) goto ErrorHandling;
}
dim=key->dim;
maxm=MAXM[dim -1];
nn = ((CE_storage*)key->S)->nn;
cov = ((CE_storage*)key->S)->cov;
param = ((CE_storage*)key->S)->param;
m = ((CE_storage*)key->S)->m;
cumm = ((CE_storage*)key->S)->cumm;
halfm = ((CE_storage*)key->S)->halfm;
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 */
for (i=0;i<dim;i++){
m[i] = 1 << (1 + (int) ceil(log(((Real) nn[i])-EPSILON1000) * INVLOG2));
if (m[i]<CIRCEMBED_MMIN) {m[i]=CIRCEMBED_MMIN;} // maybe CIRCEMBED_MMIN[i] ...
halfm[i] = m[i]/2;
index[i]=1-halfm[i]; //
if (m[i]>maxm) {error=ERRORFAILED; goto ErrorHandling;}
invn2[i] = key->x[i][XSTEP] * key->x[i][XSTEP]; /* These are the nominators in (3.1).
But here a rectangle nn[0] * step x ... x nn[dim-1] * step
is used instead of the [0,1]^d cube.
"*step" is already squared as finally the Euclidean distance
has to calculated.
Here, the notation is more general than needed (for extension
to non-isotropic random fields */
}
positivedefinite = false;
/* Eq. (3.12) shows that only j\in I(m) [cf. (3.2)] is needed,
so only the first to 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]<=m[l]/2
index[l]=h[l]-m[l] if m[l]/2+1<=h[l]<=m[l]-1
Then h[l]=(index[l]+m[l]) mod m[l] !!
*/
/* find Fourier transform of matrix; the latter might be enlarged
automatically */
trials=0;
while (!positivedefinite && (trials<CIRCEMBED_TRIALS)){
trials++;
cumm[0]=1; for(i=0;i<dim-1;i++){cumm[i+1]=cumm[i] * m[i];}
mtot=cumm[dim-1] * m[dim-1];
if (GENERAL_PRINTLEVEL>=2) PRINTF("mtot=%d\n ",mtot);
twoRealmtot = 2*sizeof(Real) * mtot;
//// for the following, see the paper by Wood and Chan!
// meaning of following variable c, see eq. (3.8)
if ((c =(Real*) malloc(twoRealmtot)) ==0){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;
}
for (i=0;i<mtot;i++){ // fill in c(h) column-wise
register int dummy;
j=0;
for (k=0;k<dim;k++) {j+=cumm[k] * ((index[k]+m[k]) % m[k]);}
// here specialisation to rotationinvariant cov.-functions:
h=0.0; // Euclidean distance for the right hand side of eq (3.8)
for (k=0;k<dim;k++) {h+=invn2[k] * (Real) (index[k] * index[k]); }
c[dummy=2*j]=cov(sqrt(h),param); c[dummy+1]=0.0; /* cf. eq. 3.8; first Real
part, then imaginary
part */
k=0;while((k<dim)&&(++index[k]>halfm[k])) {index[k]=1-halfm[k]; k++;}
}
if (GENERAL_PRINTLEVEL>=10) { // checks -- to be deleted
int M=10000;
for (i=0;i<mtot;i++) {
if ((fabs(c[2*i])>M) || (fabs(c[2*i+1])>M)) {
PRINTF("CY(%d %f %f %f %f)\n",i,c[2*i],c[2*i+1],
fabs(c[2*i])>M,fabs(c[2*i+1])>M); /* exit(1); */}
}
}
if (error=fastfourier(c,m,dim,TRUE,&(((CE_storage*)key->S)->FFT)))
goto ErrorHandling;
if (GENERAL_PRINTLEVEL>=10) { // checks -- to be deleted
int M=10000;
for (i=0;i<mtot;i++) {
if ((fabs(c[2*i])>M) || (fabs(c[2*i+1])>M)) {
PRINTF("CX(%d %f %f %f %f)\n",i,c[2*i],c[2*i+1],
fabs(c[2*i])>M,fabs(c[2*i+1])>M); /* exit(1); */
}
}
}
// check if positive definite. If not: enlarge and restart
if (!CIRCEMBED_FORCE || (trials<CIRCEMBED_TRIALS)) {
i=0; twoi=0;
while ((i<mtot)&&(positivedefinite=(c[twoi]>=CIRCEMBED_TOL_RE&&
fabs(c[twoi+1])<CIRCEMBED_TOL_IM)))
{i++; twoi+=2;}
if (!positivedefinite) {
if (GENERAL_PRINTLEVEL>=2)
PRINTF(" nonpos %d %f %f \n",i,c[twoi],c[twoi+1]);
free(c); c=NULL;
FFT_destruct(&(((CE_storage*)key->S)->FFT));
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 */
m[i] <<= 1;
halfm[i]=m[i]/2;
index[i]=1-halfm[i];
if (m[i]>maxm) {error=ERRORFAILED;goto ErrorHandling;}
}
}
} else {if (GENERAL_PRINTLEVEL>=2) PRINTF("forced\n");}
}
if (positivedefinite || CIRCEMBED_FORCE) {
// correct theoretically impossible values, that are still within
// tolerance CIRCEMBED_TOL_RE/CIRCEMBED_TOL_IM
for(i=0,twoi=0;i<mtot;i++) {
c[twoi] = (c[twoi]>0.0 ? c[twoi] : 0.0);
c[twoi+1] = 0.0;
twoi+=2;
}
} else {error=ERRORFAILED;goto ErrorHandling;}
if (GENERAL_STORING) {
if ((((CE_storage*)key->S)->d=(double *)malloc(twoRealmtot))==0){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;} //d
}
((CE_storage*)key->S)->c=c;
if (GENERAL_PRINTLEVEL>=10) {
for (i=0;i<2*mtot;i++) {PRINTF("%f ",c[i]);} PRINTF("\n");
}
return 0;
ErrorHandling:
if (c!=NULL) {free(c);}
assert(key->destruct!=NULL);
key->destruct(&key->S);
key->destruct=NULL;
return error;
}
int init_circulantembedding(key_type * key)
{
int error,i,d;
assert(key->active==false);
assert((key->covnr>=0) &&(key->covnr<currentNrCov));
assert(key->cov->cov!=NULL); // I do not know an example yet,
// where cov is unknown, but should be simulated
assert(key->param[VARIANCE]>=0.0);
assert(key->param[SILL]>=key->param[VARIANCE]);
assert( fabs(key->param[SCALE]*key->param[INVSCALE]-1.0) < EPSILON);
assert((key->S==NULL) && (key->destruct==NULL));
if ((key->S=malloc(sizeof(CE_storage)))==0){
error=ERRORMEMORYALLOCATION; goto ErrorHandling;
}
((CE_storage*)key->S)->c =NULL;
((CE_storage*)key->S)->d =NULL;
((CE_storage*)key->S)->local=NULL;
FFT_NULL(&(((CE_storage*)key->S)->FFT));
key->destruct = CE_destruct;
// this circulant embedding procedure is programmed only for regular grids
// and isotropic covariance functions
if (!key->grid) {error=ERRORMETHODNOTALLOWED;goto ErrorHandling;}
for (d=1; d<key->dim; d++){
if ((key->length[d]!=1) && // length.==1 iff grid is of lower dimension
((key->x[d][XSTEP]!=key->x[0][XSTEP]) ||
(key->length[d]!=key->length[0])))
// [d[2] is the step length
// of the grid, check which dimension
{
error=ERRORNOTPROGRAMMED; goto ErrorHandling;;
}
}
for (i=0; i<TOTAL_PARAM; i++) {((CE_storage*)key->S )->param[i]=key->param[i];}
((CE_storage*)key->S )->cov = key->cov->cov;
// this is still desastrous!--- in do_circ, an enlarged res should be used then
// so write an internal do_circ function
for (d=0; d<key->dim; d++) {
((CE_storage*)key->S)->nn[d]=key->length[d];
}
((CE_storage*)key->S)->totallength=key->totallength;
if (error=internal_init_circulantembedding(key)) goto ErrorHandling;
return 0;
ErrorHandling:
if (key->destruct!=NULL) key->destruct(&key->S);
key->destruct=NULL;
key->active = false;
return error;
}
void do_circulantembedding(key_type *key, Real *res END_WITH_GSLRNG)
/*
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 GENERAL_STORUNG==false when calling init_circulantembedding
and GENERAL_STORUNG==true when calling do_circulantembedding, the
complete programme will fail, since the initialization depends on
the value of GENERAL_STORUNG
*/
{
int i,*halfm, Ntot, k, kk,*m, *nn,
*cumm,index[MAXDIM],halfmtot,lastnonzeroORmiddle,modtot;
Real UU,VV,XX,YY,invsqrtmtot;
Real *c,*d;
bool zeroORmiddle,first;
int error,dim;
long mtot,twoRealmtot;
#ifdef RF_GSL
assert(RANDOM!=NULL);
#endif
assert(key->active);
c=((CE_storage*)key->S)->c;
m= ((CE_storage*)key->S)->m;
dim=key->dim;
halfm=((CE_storage*)key->S)->halfm;
cumm=((CE_storage*)key->S)->cumm; ////
mtot=cumm[dim-1] * m[dim-1];
twoRealmtot = 2*sizeof(Real) * mtot;
invsqrtmtot = 1/sqrt((Real)mtot);
halfmtot = mtot /2;
modtot = cumm[dim-1] * (halfm[dim-1]+1);
Ntot=((CE_storage*)key->S)->totallength;
nn = ((CE_storage*)key->S)->nn;
if (((CE_storage*)key->S)->d==NULL) {d=c;} /* overwrite the intermediate
result directly (algorithm
allows for that) */
else{d=((CE_storage*)key->S)->d;}
if (GENERAL_PRINTLEVEL>=5) 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;} // superfluent?
//// for details see Wood and Chan:
for (i=0; i<modtot; i++) { /* about half of the loop is clearly unnecessary,
as pairs of r.v. are created, see inequality * below;
therefore modtot is used instead of mtot;
as the procedure is complicated, it is easier to let the loop run
over a bit more values of i than using exactly mtot/2 values */
zeroORmiddle = true;
lastnonzeroORmiddle=0; /* theoretically, it is undefined, if
index[k] = (0 or halfm[k]) for all k;
but this would lead an error in the inequality *
below; in the case all the indices are 0 or
halfm[k] then inequality * below is obviously
satisfied with lastnonzeroORmiddle=0
*/
kk=0; /* kk takes the role of k[l] in Prop. 3; but here kk is already
transformed from matrix indices to a vector index */
for (k=0; k<dim; k++){
if (index[k]==0 || index[k]==halfm[k]) {kk += cumm[k] * index[k];}
else {
zeroORmiddle=false;
kk += cumm[k] * (m[k]-index[k]);
lastnonzeroORmiddle=k;
}
}
if (index[lastnonzeroORmiddle]<=halfm[lastnonzeroORmiddle]){ /* * */
/* By Prop. 3, pairs of r.v. exist which are (neg.) dependent;
so, these pairs are created together what is done in
the following. The above inequality guaratees, that
those pairs are created only once ---
Assume m=4 for all k and dimension=3; then the above inequality
becomes clear considering the following pair of tripels whose
corresponding random variables are simulated simultaneously.
[(0, 0, 0); (0, 0, 2)];
[(0, 1, 0); (0, 3, 0)];
[(0, 3, 0); ...] is ruled out !
[(1, 1, 1); (3, 3, 3)];
[(1, 3, 1); (3, 1, 3)];
[(3, 1, 3); ...] is ruled out!
[(3, 3, 3); ...] is ruled out!
[(0, 0, 2); ...] is NOT ruled out by *, but by the following
inequality (i<halfmtot)
*/
if (zeroORmiddle){ /* step 2 of Prop. 3 */
if (i<halfmtot) {
int twoi;
/* create two Gaussian random variables XX & YY with variance 1/2 */
UU= sqrt(-log(1.0 - UNIFORM_RANDOM));
VV = TWOPI * UNIFORM_RANDOM;
XX = UU * sin(VV);
YY = UU * cos(VV);
kk = (i + halfmtot)<<1;
twoi = i << 1;
d[twoi] = sqrt(c[twoi]) * XX; d[twoi+1]=0.0;
d[kk] = sqrt(c[kk]) * YY; d[kk+1]=0.0;
}
} else { /* step 3 of Prop. 3 */
int twoi, twok;
UU= sqrt(-log(1.0 - UNIFORM_RANDOM));
VV = TWOPI * UNIFORM_RANDOM;
XX = UU * sin(VV);
YY = UU * cos(VV);
twoi = i << 1;
twok = kk<< 1;
d[twoi+1]=d[2*i]=sqrt(c[twoi]); d[twoi]*=XX; d[twoi+1]*=YY;
d[twok+1]=d[2*kk]=sqrt(c[twok]); d[twok]*=XX; d[twok+1]*= -YY;
}
}
k=0; while((k<dim) && (++index[k]>=m[k])) {index[k++]=0;}
}
fastfourier(d,m,dim,FALSE,&(((CE_storage*)key->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;}
for (i=0;i<Ntot;i++){
kk=0;for (k=0; k<dim; k++){kk+=cumm[k] * index[k];}
res[i]=d[2*kk]*invsqrtmtot;
if ((fabs(d[2*kk+1])>CIRCEMBED_TOL_IM) &&
((GENERAL_PRINTLEVEL>=2 && first) || GENERAL_PRINTLEVEL>=3)){
PRINTF("IMAGINARY PART <> 0, %f\n",d[2*kk+1]); first=false;
}
k=0; while((k<dim) && (++index[k]>=nn[k])) {index[k++]=0; }
}
if (!(key->active=(GENERAL_STORING && (((CE_storage*)key->S)->d!=NULL)))) {
assert(key->destruct!=NULL);
key->destruct(&key->S);
key->destruct=NULL;
}
}
// see /home/martin/article/C/RFce_local.cc
int init_circ_embed_local(key_type * key) {assert(false); /* not public yet */}
// make sure that internal_init is called with a quadratic scheme
// make sure that the relevant part is cut out correctly
// hinr?! nicht doppeltes Feld generieren, sondern 1faches?
// da sowieso nur 1/3 verwandt wird??
// nein. an sich schon doppeltes Feld; aber bei kompakten Traeger [0,1], muss
// maxabstand der Punkte nur 1/2 betragen (durch Verdoppelung der Matrix
// wird range erreicht oder so aehnlich).
// --> beruecksichtigung bei lokaler simuation
//
// for future development : make sure that the directions are set up correctly!
void do_circ_embed_local(key_type *key, Real *res END_WITH_GSLRNG) {
assert(false); /* not public yet */
}