Revision 6aa7805ff1eef8ecb881b5099f3cd30f7c2bd637 authored by Martin Schlather on 08 August 1977, 00:00:00 UTC, committed by Gabor Csardi on 08 August 1977, 00:00:00 UTC
1 parent a3c58bc
RFgetNset.cc
/*
Authors
Martin Schlather, martin.schlather@cu.lu
library for simulation of random fields -- get's, set's, and print's
Copyright (C) 2001 -- 2004 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 <assert.h>
#include <string.h>
#include "RFsimu.h"
#include "MPPFcts.h"
#include <R_ext/Applic.h>
void GetParameterIndices(int *variance, int *scale,
int *kappa, int *lastkappa, int *invscale,
int *aniso)
{
*variance=VARIANCE; *scale=SCALE; *kappa=KAPPA;
*lastkappa=LASTKAPPA; *invscale=INVSCALE; *aniso=ANISO;
}
void GetrfParameters(int *covmaxchar, int *methodmaxchar,
int *distrmaxchar,
int *covnr, int *methodnr, int *distrnr,
int *maxdim, int *maxmodels) {
if (currentNrCov==-1) InitModelList();
*covmaxchar=COVMAXCHAR;
*methodmaxchar=METHODMAXCHAR;
*distrmaxchar=DISTRMAXCHAR;
*covnr=currentNrCov;
*methodnr=(int) Nothing;
*distrnr=DISTRNR;
*maxdim=MAXDIM;
*maxmodels=MAXKEYS;
}
void GetDistrName(int *nr,char **name){
if ((*nr<0) ||(*nr>=DISTRNR)) {strncpy(*name,"",DISTRMAXCHAR); return;}
strncpy(*name, DISTRNAMES[*nr], DISTRMAXCHAR);
}
void GetDistrNr(char **name, int *n, int *nr)
{
unsigned int ln,v;
// == -1 if no matching function is found
// == -2 if multiple matching fnts are found, without one matching exactly
// if more than one match exactly, the last one is taken (enables overwriting
// standard functions)
for (v=0; v<*n; v++) {
nr[v]=0;
ln=strlen(name[v]);
while ( (nr[v] < DISTRNR) && strncmp(name[v], DISTRNAMES[nr[v]], ln)) {
(nr[v])++;
}
if (nr[v]<DISTRNR) {
// a matching function is found. Are there other functions that match?
int j;
bool exactmatching,multiplematching;
exactmatching=(ln==strlen(DISTRNAMES[nr[v]]));
multiplematching=false;
j=nr[v]+1; // if two or more distributions have the same name
// the last one is taken
while (j<DISTRNR) {
while ( (j<DISTRNR) && strncmp(name[v],DISTRNAMES[j],ln)) {j++;}
if (j<DISTRNR) {
if (ln==strlen(DISTRNAMES[j])) {nr[v]=j; exactmatching=true;}
else {multiplematching=true;}
}
j++;
}
if (!exactmatching && multiplematching) {nr[v]=-2;}
} else nr[v]=-1;
}
}
void SetParamCE( int *action, int *force, Real *tolRe, Real *tolIm,
int *trials,
int *mmin, // mmin must be vector of length MAXDIM!!
int *userfft, int *strategy,
ce_param *CE, char *name)
{
int d;
switch(*action) {
case 0 :
CE->force=(bool)*force;
CE->tol_re=*tolRe;
if (CE->tol_re>0) {
CE->tol_re=0;
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! %s.tol_re had been positive.\n",name);
}
CE->tol_im=*tolIm;
if (CE->tol_im<0) {
CE->tol_im=0;
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! %s.tol_im had been neagtive.\n",name);
}
CE->trials=*trials;
if (CE->trials<1) {
CE->trials=1;
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! %s.trials had been less than 1\n",name);
}
for (d=0; d<MAXDIM; d++) {
CE->mmin[d] = mmin[d];
if (CE->mmin[d]>0) {
CE->mmin[d] =
1 << (1+(int) (log(((double) CE->mmin[d])-0.1)/log(2.0)));
if ((CE->mmin[d]!=mmin[d]) && (GENERAL_PRINTLEVEL>0))
PRINTF("\nWARNING! %s.mmin[%d] set to %d.\n",name,d,CE->mmin[d]);
}
}
CE->userfft=(bool)*userfft;
if (*strategy>LASTSTRATEGY) {
if (GENERAL_PRINTLEVEL>0)
PRINTF("\nWARNING! %s_STRATEGY not set\n",name);
} else CE->strategy=(char) *strategy;
break;
case 1 :
*force=CE->force;
*tolRe=CE->tol_re;
*tolIm=CE->tol_im;
*trials=CE->trials;
for (d=0; d<MAXDIM; d++) {mmin[d]=CE->mmin[d];}
*userfft= CE->userfft;
*strategy=(int) CE->strategy;
if (GetNotPrint) break;
case 2 :
PRINTF("\n%s\n==========\nforce=%d\ntol_re=%e\ntol_im=%e\ntrials=%d\nuserfft=%d\nstrategy=%d\n",name,
CE->force,CE->tol_re,CE->tol_im, CE->trials, CE->userfft,
(int) CE->strategy);
for (d=0; d<MAXDIM; d++) PRINTF("%d ", CE->mmin[d]);
PRINTF("\n");
break;
default : PRINTF(" unknown action\n");
}
}
void GetModelName(int *nr,char **name){
if (currentNrCov==-1) InitModelList();
if ((*nr<0) ||(*nr>=currentNrCov)) {strncpy(*name,"",COVMAXCHAR); return;}
strncpy(*name,CovList[*nr].name,COVMAXCHAR);
}
void GetNrParameters(int *covnr, int *n, int* kappas) {
int v;
if (currentNrCov==-1) InitModelList();
for (v=0; v<*n; v++) {
if ((covnr[v]<0) ||(covnr[v]>=currentNrCov)) {kappas[v]=-1;}
else kappas[v] = CovList[covnr[v]].kappas;
}
}
void GetModelNr(char **name, int *n, int *nr)
{
unsigned int ln,v;
// == -1 if no matching function is found
// == -2 if multiple matching fnts are found, without one matching exactly
// if more than one match exactly, the last one is taken (enables overwriting
// standard functions)
if (currentNrCov==-1) InitModelList();
for (v=0; v<*n; v++) {
nr[v]=0;
ln=strlen(name[v]);
while ( (nr[v]<currentNrCov) && strncmp(name[v],CovList[nr[v]].name,ln)) {
(nr[v])++;
}
if (nr[v]<currentNrCov) {
// a matching function is found. Are there other functions that match?
int j;
bool exactmatching,multiplematching;
exactmatching=(ln==strlen(CovList[nr[v]].name));
multiplematching=false;
j=nr[v]+1; // if two or more covariance functions have the same name
// the last one is taken
while (j<currentNrCov) {
while ( (j<currentNrCov) && strncmp(name[v],CovList[j].name,ln)) {j++;}
if (j<currentNrCov) {
if (ln==strlen(CovList[j].name)) {nr[v]=j; exactmatching=true;}
else {multiplematching=true;}
}
j++;
}
if (!exactmatching && multiplematching) {nr[v]=-2;}
} else nr[v]=-1;
}
}
void GetMethodNr(char **name, int *n, int *nr)
{
unsigned int ln, v;
// == -1 if no matching method is found
// == -2 if multiple matching methods are found, without one matching exactly
for (v=0; v<*n; v++) {
nr[v]=0;
ln=strlen(name[v]);
while ( (nr[v]<(int) Forbidden) && strncmp(name[v],METHODNAMES[nr[v]],ln)) {
(nr[v])++;
}
if (nr[v]<(int) Forbidden) {
// a matching method is found. Are there other methods that match?
int j;
bool exactmatching,multiplematching;
exactmatching=(ln==strlen(METHODNAMES[nr[v]]));
multiplematching=false;
j=nr[v]+1; // if two or more methods have the same name the last one is
// taken; stupid here, but nice in GetCovNr
while (j<(int) Forbidden) {
while ( (j<(int) Forbidden) && strncmp(name[v],METHODNAMES[j],ln)) j++;
if (j<(int) Forbidden) {
if (ln==strlen(METHODNAMES[j])) {nr[v]=j; exactmatching=true;}
else {multiplematching=true;}
}
j++;
}
if (!exactmatching && multiplematching) {nr[v]=-2;}
} else nr[v]=-1;
}
}
void GetMethodName(int *nr, char **name)
{
if ((*nr<0) ||(*nr>=(int) Forbidden)) {
strncpy(*name,"",METHODMAXCHAR);
return;
}
strncpy(*name,METHODNAMES[*nr],METHODMAXCHAR);
}
void PrintMethods()
{
int i;
PRINTF("\n\n Methods for generating Gaussian random fields\n =============================================\n\n");
for (i=0;i<(int) Nothing;i++) { PRINTF(" * %s\n",METHODNAMES[i]); }
PRINTF("\n\n Methods for non-Gaussian random fields\n ======================================\n");
for (i=1+(int) Nothing;i<(int) Forbidden; i++) {
PRINTF(" * %s\n",METHODNAMES[i]); }
PRINTF("\n");
}
void PrintCovDetailed(int i)
{
PRINTF("%5d: %s cov=%d %d %d, tbm=%d %d %d, spec=%d,\n abf=%d hyp=%d, oth=%d %d\n",
i,CovList[i].name,CovList[i].cov,CovList[i].cov_loc,
CovList[i].square_factor,
CovList[i].cov_tbm2,CovList[i].cov_tbm3,CovList[i].tbm_method,
CovList[i].spectral,CovList[i].add_mpp_scl,CovList[i].hyperplane,
CovList[i].initother,CovList[i].other);
}
void PrintModelListDetailed()
{
int i;
if (CovList==NULL) {PRINTF("There are no functions available!\n");}
else {
PRINTF("\nList of covariance functions:\n=============================\n");
for (i=0;i<currentNrCov;i++) PrintCovDetailed(i);
}
}
void GetRange(int *nr, int *dim, int *index, Real *range, int *lrange){
// never change Real without crosschecking with fcts in RFCovFcts.cc!
// index is increased by one except index is the largest value possible
int i;
if (currentNrCov==-1) InitModelList();
if ((*nr<0) || (*nr>=currentNrCov) ||
(*lrange != CovList[*nr].kappas * 4) ) {
for (i=0; i<*lrange; i++) range[i]=RF_NAN;
*index = -100;
return;
}
assert(CovList[*nr].range != NULL);
CovList[*nr].range(*dim, index, range);
// index>0 : further parts of the range are missing
// index=-1: no further part of the range
// index=-2: dimension not valid
//assert(false);
}
extern int IncludeModel(char *name, int kappas, checkfct check,
int isotropic, bool variogram, methodfct method,
rangefct range)
{
int i;
if (currentNrCov==-1) InitModelList(); assert(CovList!=NULL);
assert(currentNrCov>=0);
if (currentNrCov>=MAXNRCOVFCTS)
{PRINTF("Error. Model list full.\n");return -1;}
strncpy(CovList[currentNrCov].name, name, COVMAXCHAR-1);
CovList[currentNrCov].name[COVMAXCHAR-1]='\0';
if (strlen(name)>=COVMAXCHAR) {
PRINTF("Warning! Covariance name is truncated to `%s'",
CovList[currentNrCov].name);
}
CovList[currentNrCov].kappas = kappas;
CovList[currentNrCov].check = check;
CovList[currentNrCov].isotropic = isotropic;
CovList[currentNrCov].cov=NULL;
CovList[currentNrCov].naturalscale=NULL;
CovList[currentNrCov].cov_loc=NULL;
CovList[currentNrCov].cov_tbm2=NULL;
CovList[currentNrCov].cov_tbm3=NULL;
CovList[currentNrCov].ableitung=NULL,
CovList[currentNrCov].spectral=NULL;
CovList[currentNrCov].add_mpp_scl=NULL;
CovList[currentNrCov].add_mpp_rnd=NULL;
CovList[currentNrCov].hyperplane=NULL;
CovList[currentNrCov].initother=NULL;
CovList[currentNrCov].other=NULL;
CovList[currentNrCov].square_factor=NULL;
CovList[currentNrCov].tbm_method = Forbidden;
CovList[currentNrCov].isotropic = isotropic;
CovList[currentNrCov].variogram = variogram;
CovList[currentNrCov].even = true;
CovList[currentNrCov].range= range;
assert((CovList[currentNrCov].method=method)!=NULL);
for (i=0; i<MAXDIM; i++) {
CovList[currentNrCov].odd[i] = false;
}
if (kappas>0) assert(check!=NULL);
currentNrCov++;
return(currentNrCov-1);
}
extern void addodd(int nr, int dim) {
assert((dim>=0) && (dim<MAXDIM));
assert((nr>=0) && (nr<currentNrCov));
assert(CovList[nr].isotropic==ANISOTROPIC);
PRINTF("addodd not implemented yet. Sorry.");
assert(false);
CovList[nr].even = false;
CovList[nr].odd[dim] = true;
}
//extern void addSimu(int nr, SimulationType r1, SimulationType r2,
// SimulationType r3)
// // see RFsimu.h for comments
//{
// int i;
// assert((nr>=0) && (nr<currentNrCov));
// CovList[nr].first[0]=r1;
// CovList[nr].first[1]=r2;
// for (i=2; i<MAXDIM; i++) CovList[nr].first[i]=r3;
//}
extern void addCov(int nr, covfct cov, natscalefct naturalscale,
covfct cov_loc, scalefct square_factor)
{
assert((nr>=0) && (nr<currentNrCov));
CovList[nr].cov=cov;
CovList[nr].naturalscale=naturalscale;
CovList[nr].cov_loc=cov_loc;
CovList[nr].square_factor=square_factor;
if ((square_factor!=NULL) || (cov_loc!=NULL) )
assert( (square_factor!=NULL) && (cov_loc!=NULL) );
}
extern void addTBM(int nr, isofct cov_tbm2, isofct cov_tbm3,
isofct ableitung,
SimulationType tbm_method,
randommeasure spectral) {
// must be called always AFTER addCov !!!!
assert((nr>=0) && (nr<currentNrCov));
CovList[nr].cov_tbm2=cov_tbm2;
CovList[nr].cov_tbm3=cov_tbm3;
CovList[nr].ableitung=ableitung,
CovList[nr].tbm_method = tbm_method;
CovList[nr].spectral=spectral;
if ( (ableitung!=NULL) || (cov_tbm3!=NULL) || (tbm_method!=Nothing))
assert((ableitung!=NULL) && (cov_tbm3!=NULL) && (CovList[nr].cov!=NULL)
&& (tbm_method!=Nothing));
}
extern void addOther(int nr, initmppfct add_mpp_scl, MPPRandom add_mpp_rnd,
generalSimuInit initother, generalSimuMethod other){
assert((nr>=0) && (nr<currentNrCov));
CovList[nr].add_mpp_scl=add_mpp_scl;
CovList[nr].add_mpp_rnd=add_mpp_rnd;
CovList[nr].initother=initother;
CovList[nr].other=other;
CovList[nr].hyperplane=NULL;
if ((add_mpp_scl!=NULL) || (add_mpp_rnd!=NULL))
assert((add_mpp_scl!=NULL) && (add_mpp_rnd!=NULL));
if ((other!=NULL) || (initother!=NULL))
assert((other!=NULL) && (initother!=NULL));
}
void PrintModelList()
{
int i;
char percent[]="%";
char empty[]="";
char header[]="Circ local TBM2 TBM3 sp dir add hyp oth\n";
char coded[3][2]={"-","X","+"};
char firstcolumn[20],line[80];
if (currentNrCov==-1) {
InitModelList();
if (GENERAL_PRINTLEVEL>5)
PRINTF("List of covariance functions initiated.\n");
}
if (CovList==NULL) {PRINTF("There are no functions available!\n");}
else {
sprintf(firstcolumn,"%s%ds ",percent,COVMAXCHAR);
sprintf(line,"%s3s %s3s %s3s %s3s %s2s %s2s %s2s %s2s %s2s\n",
percent,percent,percent,percent,percent,percent,percent,percent,
percent);
PRINTF("\n\n");
PRINTF(firstcolumn,empty); PRINTF(" List of models\n");
PRINTF(firstcolumn,empty); PRINTF(" ==============\n");
PRINTF(firstcolumn,empty); PRINTF("[See also PrintMethodList()]\n\n");
PRINTF(firstcolumn,empty);PRINTF(header,empty);
for (i=0;i<currentNrCov;i++) {
PRINTF(firstcolumn,CovList[i].name);
if (strcmp(CovList[i].name,"nugget")==0) {
PRINTF(line,
coded[2], coded[false], coded[2], coded[2], coded[2],
coded[2], coded[2], coded[false], coded[false]);
} else {
PRINTF(line,
coded[(CovList[i].cov!=NULL) && (CovList[i].cov_loc==NULL)],
coded[CovList[i].cov_loc!=NULL],
coded[CovList[i].cov_tbm2!=NULL],
coded[CovList[i].cov_tbm3!=NULL],
coded[CovList[i].spectral!=NULL],
coded[(CovList[i].cov!=NULL) && (CovList[i].cov_loc==NULL)],
// nugget is left out
coded[CovList[i].add_mpp_scl!=NULL],
coded[CovList[i].hyperplane!=NULL],
coded[CovList[i].initother!=NULL]
);
}
}
PRINTF(firstcolumn,empty);PRINTF(header,empty);
}
}
void GetModelList(int* idx) {
int i, j, methods;
if (currentNrCov==-1) {
InitModelList();
if (GENERAL_PRINTLEVEL>5)
PRINTF("List of covariance functions initiated.\n");
}
if (CovList==NULL) return;
assert(Nothing==10);
methods = Nothing - 2; // nugget and nothing itself
for (j=i=0; i<currentNrCov; i++) {
idx[j++] = (CovList[i].cov!=NULL) && (CovList[i].cov_loc==NULL);
idx[j++] = CovList[i].cov_loc!=NULL;
idx[j++] = CovList[i].cov_tbm2!=NULL;
idx[j++] = CovList[i].cov_tbm3!=NULL;
idx[j++] = CovList[i].spectral!=NULL;
idx[j++] = (CovList[i].cov!=NULL) && (CovList[i].cov_loc==NULL);
// nugget is left out;
idx[j++] = CovList[i].add_mpp_scl!=NULL;
idx[j++] = CovList[i].hyperplane!=NULL;
idx[j++] = CovList[i].initother!=NULL;
}
return;
}
void GetNaturalScaling(int *covnr, Real *q, int *naturalscaling,
Real *natscale, int *error)
{
// q: kappas only
// values of naturalscaling:
//#define NATSCALE_EXACT 1
//#define NATSCALE_APPROX 2
//#define NATSCALE_MLE 3 /* check mleRF when changing !! */
// +10 if numerical is allowed
static int oldcovnr = -99;
static Real oldp[TOTAL_PARAM], p[TOTAL_PARAM];
static Real OldNatScale;
int TEN, endfor, k;
bool numeric;
*error = 0;
TEN = 10;
endfor = KAPPA + CovList[*covnr].kappas;
for (k=KAPPA; k<endfor; k++) p[k]=q[k-KAPPA];
if (*naturalscaling) {
if (CovList[*covnr].isotropic!=FULLISOTROPIC) {
*error=ERRORANISOTROPIC; return;}
if (!(numeric=*naturalscaling>TEN)) {TEN=0;}
*natscale=0.0;
if (CovList[*covnr].naturalscale!=NULL) {
*natscale = CovList[*covnr].naturalscale(p,*naturalscaling-TEN);
if (*natscale!=0.0) return;
}
if (numeric) {
Real x,newx,yold,y,newy;
covfct cov;
int parami, wave,i;
if ((cov=CovList[*covnr].cov)==NULL) {*error=ERRORNOTDEFINED;return;}
if ((cov=CovList[*covnr].cov)==nugget) {*error=ERRORRESCALING; return;}
// already calculated ?
parami=KAPPA; // do not compare mean,variance, etc.
if (oldcovnr==*covnr) {
for (; parami<endfor; parami++) {
if (oldp[parami]!=p[parami]) {
break;
}
}
if (parami==endfor) {
*natscale=OldNatScale;
return;
}
}
// the other parameters need not to be copied as checked that
// they are identical to previous ones
for (;parami<endfor;parami++) {oldp[parami]=p[parami];}
oldp[VARIANCE] = oldp[SCALE] = oldp[INVSCALE] = 1.0;
oldcovnr = -99; /* if error occurs, the next call will realise that
the previous result cannot be used */
/* **************************************
Now, find a quick and good solution for NewInvScale --
*/
wave = 0;
x = 1.0;
if ( (yold=cov(&x,oldp,1)) > 0.05) {
Real leftx;
x *= 2.0;
while ( (y=cov(&x,oldp,1)) > 0.05) {
if (yold<y){ wave++; if (wave>10) {*error=ERRORWAVING; return;} }
yold = y;
x *= 2.0;
if (x>1E30) {*error=ERRORRESCALING; return;} // or use a separate ERROR
}
leftx = x * 0.5;
for (i=0; i<3 /* good choice?? */ ;i++) {
if (y==yold) {*error=ERRORWAVING; return;} // should never appear
newx = x + (x-leftx)/(y-yold)*(0.05-y);
if ( (newy=cov(&newx,oldp,1)) >0.05) {
leftx = newx;
yold = newy;
} else {
x = newx;
y = newy;
}
}
if (y==yold) {*error=ERRORWAVING; return;} // should never appear
*natscale = 1.0 / ( x + (x-leftx)/(y-yold)*(0.05-y) );
} else {
Real rightx;
x *= 0.5;
while ( (y=cov(&x,oldp,1)) < 0.05) {
if (yold>y){ wave++; if (wave>10) {*error=ERRORWAVING; return;} }
yold = y;
x *= 0.5;
if (x<1E-30) {*error=ERRORRESCALING; return;} //or use a separate ERROR
}
rightx = x * 2.0;
for (i=0; i<3 /* good choice?? */ ;i++) {
if (y==yold) {*error=ERRORWAVING; return;} // should never appear
newx = x + (x-rightx)/(y-yold)*(0.05-y);
if ( (newy=cov(&newx,oldp,1)) <0.05) {
rightx = newx;
yold = newy;
} else {
x = newx;
y = newy;
}
}
if (y==yold) {*error=ERRORWAVING; return;} // should never appear
*natscale = 1.0 / ( x + (x-rightx)/(y-yold)*(0.05-y) );
}
oldcovnr=*covnr; /* result is reliable -> stored for next call */
OldNatScale = *natscale;
} else {*error=ERRORRESCALING; }
} else {
*natscale = 1.0;
}
}
void SetParam(int *action, int *storing,int *printlevel,int *naturalscaling,
char **pch) {
switch (*action) {
case 0 :
GENERAL_STORING= (bool) *storing;
GENERAL_PRINTLEVEL=*printlevel;
GENERAL_NATURALSCALING=*naturalscaling;
if ((strlen(*pch)>1) && (GENERAL_PRINTLEVEL>0))
PRINTF("\n`pch' has more than one character -- first character taken only\n");
strncpy(GENERAL_PCH, *pch, 1);
break;
case 1 :
*storing = GENERAL_STORING;
*printlevel = GENERAL_PRINTLEVEL;
*naturalscaling = GENERAL_NATURALSCALING;
strcpy(*pch, GENERAL_PCH);
if (GetNotPrint) break;
case 2 :
PRINTF("\nGeneral Parameters\n==================\nstoring=%d\nprint level=%d\nnatural scaling=%d\npch=%s\n",
GENERAL_STORING,GENERAL_PRINTLEVEL,GENERAL_NATURALSCALING,*pch);
break;
default : PRINTF(" unknown action\n");
}
}
void StoreTrend(int *keyNr, int *modus, char **trend,
int *lt, Real *lambda, int *ll, int *error)
{
unsigned long len;
key_type *key;
if (currentNrCov==-1) InitModelList(); assert(CovList!=NULL);
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {*error=ERRORREGISTER;goto ErrorHandling;}
key = &(KEY[*keyNr]);
if (key->LinearTrend!=NULL) free(key->LinearTrend);
key->LinearTrend = NULL;
if (key->TrendFunction!=NULL) free(key->TrendFunction);
key->TrendFunction = NULL;
key->lTrendFct = 0;
key->lLinTrend = 0;
switch(*modus){
case 0 : break;
case 3 : case 1 :
key->LinearTrend = (Real *) malloc(len = sizeof(Real) * *ll);
memcpy(key->LinearTrend, lambda, len);
if (*modus==1) break;
key->lLinTrend = *ll;
case 2 :
key->TrendFunction = (char *) malloc(len = sizeof(char) * *lt);
memcpy(key->TrendFunction, *trend, len);
key->lTrendFct = *lt;
break;
default: key->TrendModus = -1; *error=ERRORUNSPECIFIED; goto ErrorHandling;
}
key->TrendModus = *modus;
*error = 0;
return;
ErrorHandling:
return;
}
void GetTrueDim(bool anisotropy, int timespacedim, Real* param,
int *TrueDim, bool *no_last_comp,
int *start_aniso, int *index_dim)
// output: start_aniso (where the non-zero positions start, anisotropy)
// index_dim (dimension nr of the non-zero positions, aniso)
// no_last_comp(last dim of aniso matrix. is identically zero, i.e,
// effectively no time component)
// TrueDim The reduced dimension, including the time direction
{
int i,j,k;
unsigned long endfor, startfor;
if (anisotropy) {
/* check whether the dimension can be reduced */
/* TODO:this algorithm only detects zonal anisotropies -- could be improved*/
for (j=-1, i=0; i<timespacedim; i++) {
startfor = i * timespacedim + ANISO;
endfor= startfor + timespacedim;
for (k=startfor; k<endfor; k++) {
if (param[k]!=0.0) break;
}
if (k<endfor) {
start_aniso[++j]=startfor; /* where the non-zero positions*/
/* of the anisotropy matrix start */
index_dim[j] = i;
}
}
*no_last_comp = k>=endfor; // last i is time, if k=endfor then all
// components zero
*TrueDim = j+1;
} else {
*no_last_comp = true;
*TrueDim = timespacedim; // == spatialdim
}
}
void GetTrendLengths(int *keyNr, int *modus, int *lt, int *ll, int *lx)
{
key_type *key;
if (currentNrCov==-1) {*modus=-1;goto ErrorHandling;}
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {*modus=-2; goto ErrorHandling;}
key = &(KEY[*keyNr]);
if (!key->active) {*modus=-3; goto ErrorHandling;}
*modus = key->TrendModus;
*lt = key->lTrendFct;
*ll = key->lLinTrend;
*lx = key->totalpoints;
return;
ErrorHandling:
return;
}
void GetTrend(int *keyNr, char **trend, Real *lambda, Real *x, int *error)
{
key_type *key;
if (currentNrCov==-1) {*error=-1;goto ErrorHandling;}
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {*error=-2; goto ErrorHandling;}
key = &(KEY[*keyNr]);
if (!key->active) {*error=-3; goto ErrorHandling;}
switch(key->TrendModus){
case 0 : break;
case 3 : case 1 :
memcpy(lambda, key->LinearTrend, key->lLinTrend);
if (key->TrendModus==1) break;
case 2:
strcpy(*trend, key->TrendFunction);
break;
default: *error=-4; goto ErrorHandling;
}
*error = 0;
return;
ErrorHandling:
return;
}
void GetCoordinates(int *keyNr, Real *x, int *error)
{
key_type *key;
if (currentNrCov==-1) {*error=-1;goto ErrorHandling;}
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {*error=-2; goto ErrorHandling;}
key = &(KEY[*keyNr]);
if (!key->active) {*error=-3; goto ErrorHandling;}
if (key->TrendModus==0) {*error=-4; goto ErrorHandling;}
*error = ERRORNOTPROGRAMMED; // TODO !!
return;
ErrorHandling:
return;
}
void GetKeyInfo(int *keyNr, int *total, int *lengths, int *spatialdim,
int *timespacedim, int *grid, int *distr, int *maxdim)
{
int d;
key_type *key;
if (*maxdim<MAXDIM) { *total=-3; *maxdim=MAXDIM; return;}
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {*total=-1; return;}
key = &(KEY[*keyNr]);
if (!key->active) {*total=-2;}
else {
*total=key->totalpoints;
*spatialdim = key->spatialdim;
*timespacedim = key->timespacedim;
for (d=0; d<MAXDIM; d++) lengths[d]=key->length[d];
*grid = (int) key->grid;
*distr = (int) key->distribution;
*maxdim = MAXDIM;
}
}
void GetCornersOfGrid(key_type *key, int Stimespacedim, int* start_aniso,
Real *param, Real *sxx){
// returning sequence (#=2^MAXDIM) vectors of dimension s->simutimespacedim
Real sx[ZWEIHOCHMAXDIM * MAXDIM];
int index,i,k,g,endforM1;
long endfor,indexx;
char j[MAXDIM];
endfor = 1 << key->timespacedim;
endforM1 = endfor -1;
for (i=0; i<key->timespacedim; i++) j[i]=0;
for (index=i=0; i<endfor; i++) {
for (k=0; k<key->timespacedim; k++) {
sx[index] = key->x[k][XSTART];
if (j[k])
sx[index] += key->x[k][XSTEP] * (double) (key->length[k]-1);
index++;
}
k=0; (j[k])++;
if (i!=endforM1)
while(j[k]>=2) { assert(j[k]==2); j[k]=0; k++; (j[k])++;}
}
for (index=indexx=i=0; i<endfor;
i++, index+=key->timespacedim, indexx+=Stimespacedim) {
for (k=0; k<Stimespacedim; k++) {
register Real dummy;
dummy = 0.0;
for (g=0; g<key->timespacedim; g++)
dummy += param[start_aniso[k]+g] * sx[index+g];
sxx[indexx + k] = dummy;
}
}
}
void GetCornersOfElement(key_type *key, int Stimespacedim, int* start_aniso,
Real *param, Real *sxx){
Real sx[ZWEIHOCHMAXDIM * MAXDIM];
int index,i,k,g,endforM1;
long endfor,indexx;
char j[MAXDIM];
endfor = 1 << key->timespacedim;
endforM1 = endfor -1;
for (i=0; i<key->timespacedim; i++) j[i]=0;
for (index=i=0; i<endfor; i++) {
for (k=0; k<key->timespacedim; k++) {
if (j[k]) sx[index] = key->x[k][XSTEP]; else sx[index]=0.0;
index++;
}
k=0; (j[k])++;
if (i!=endforM1)
while(j[k]>=2) { assert(j[k]==2); j[k]=0; k++; (j[k])++;}
}
for (index=indexx=i=0; i<endfor;
i++, index+=key->timespacedim, indexx+=Stimespacedim) {
for (k=0; k<Stimespacedim; k++) {
register Real dummy;
dummy = 0.0;
for (g=0; g<key->timespacedim; g++)
dummy += param[start_aniso[k]+g] * sx[index+g];
sxx[indexx + k] = dummy;
}
}
}
void GetRangeCornerDistances(key_type *key, Real *sxx, int Stimespacedim,
int Ssimuspatialdim, Real *min, Real *max) {
long endfor,indexx,index;
int d,i,k;
endfor = 1 << key->timespacedim;
*min = RF_INF;
*max = RF_NEGINF;
for (index=i=0; i<endfor; i++, index+=Stimespacedim)
for (indexx=k=0; k<i; k++, indexx+=Stimespacedim) {
register double sum = 0.0, dummy;
for (d=0; d<Ssimuspatialdim; d++) {
dummy = sxx[index+d] - sxx[indexx+d];
sum += dummy * dummy;
}
if ((sum > 0) && (sum < *min)) *min = sum;
if ((sum > 0) && (sum > *max)) *max = sum;
}
*min = sqrt(*min);
*max = sqrt(*max);
}
void niceFFTnumber(int *n) {
*n = (int) NiceFFTNumber( (unsigned long) *n);
}
bool ok_n(int n, int *f, int nf) // taken from fourier.c of R
{
int i;
for (i = 0; i < nf; i++)
while(n % f[i] == 0) if ((n = n / f[i]) == 1) return TRUE;
return n == 1;
}
int nextn(int n, int *f, int nf) // taken from fourier.c of R
{
while(!ok_n(n, f, nf)) n++;
return n;
}
bool HOMEMADE_NICEFFT=true;
unsigned long NiceFFTNumber(unsigned long n) {
unsigned long i,ii,j,jj,k,kk,l,ll,min, m=1, f[4]={2,3,5,7};
if (HOMEMADE_NICEFFT) {
if (n<=1) return n;
for (i=0; i<4; i++)
while ( ((n % f[i])==0) && (n>10000)) { m*=f[i]; n/=f[i]; }
if (n>10000) {
while (n>10000) {m*=10; n/=10;}
n++;
}
min = 10000000;
for (i=0, ii=1; i<=14; i++, ii<<=1) { // 2^14>10.000
if (ii>=n) { if (ii<min) min=ii; break; }
for (j=0, jj=ii; j<=9; j++, jj*=3) {// 3^9>10.000
if (jj>=n) { if (jj<min) min=jj; break; }
//for (k=0, kk=jj; k<=6; k++, kk*=5) {// 5^6>10.000
//if (kk>=n) { if (kk<min) min=kk; break; }
//for (l=0, ll=kk; l<=5; l++, ll*=7) {// 7^5>10.000
// instead of (if 7 is included)
for (l=0, ll=jj; l<=6; l++, ll*=5) {// 5^5>10.000
if (ll>=n) {
if (ll<min) min=ll;
break;
}
//}
}
}
}
return m*min;
} else { // not HOMEMADE_NICEFFT
#define F_NUMBERS 3
int f[F_NUMBERS]={2,3,5};
return nextn(n, f, F_NUMBERS);
}
}
int eigenvalues(Real *C, int dim, double *ev)
{
// SVD decomposition
double *V,*e, *U, *G, *D;
long longrow = dim, job=11, Error;
int error, d;
error=0;
D = V = e = U = G = NULL;
if ((U =(double *) malloc(sizeof(double) * longrow * longrow))==NULL){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;
}
if ((G = (double *) malloc(sizeof(double) * longrow))==NULL){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;
}
if ((D = (double *) malloc(sizeof(double) * longrow * longrow))==NULL){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;
}
if ((V =(double *) malloc(sizeof(double) * longrow * longrow))==NULL){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;
}
if ((e = (double *) malloc(sizeof(double) * longrow))==NULL){
error=ERRORMEMORYALLOCATION;goto ErrorHandling;
}
for (d=dim * dim -1; d>=0; d--) D[d]=C[d]; // dsvdc destroys the
// input matrix !!!!!!!!!!!!!!!!!!!!
{
int row,jobint;
row=longrow; jobint=job;
F77_CALL(dsvdc)(D,&row,&row,&row,ev,e,U,&row,V,&row,G,&jobint,
&error);
}
ErrorHandling:
if (D!=NULL) free(D);
if (V!=NULL) free(V);
if (e!=NULL) free(e);
if (U!=NULL) free(U);
if (G!=NULL) free(G);
return error;
}
Computing file changes ...
