https://github.com/cran/RandomFields
Tip revision: 4877e49dad8ee6b04e79289f69ff7f2186f11506 authored by Martin Schlather on 20 January 2012, 00:00:00 UTC
version 2.0.54
version 2.0.54
Tip revision: 4877e49
getNset.cc
/*
Authors
Martin Schlather, martin.schlather@math.uni-goettingen.de
(library for simulation of random fields)
Copyright (C) 2001 -- 2011 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.
RO
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.
*/
extern "C" {
#include <R.h>
#include <Rdefines.h>
}
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "RF.h"
#include "primitive.h"
#include <R_ext/Linpack.h>
void KEY_DELETE(key_type *key) {
KEY_DELETE_NOTTREND(key);
TREND_DELETE(&(key->trend));
}
void KEY_NULL(key_type *key)
{
KEY_NULLNOTTREND(key);
TREND_NULL(&(key->trend));
}
void KEY_DELETE_NOTTREND(key_type *key)
{
LOC_DELETE(&(key->loc));
METHOD_DELETE(&(key->meth));
// PrintModelInfo(key->cov);
if (key->cov!=NULL) COV_DELETE(&(key->cov));
// KEY_NULLNOTTREND(key);
}
void KEY_NULLNOTTREND(key_type *key) {
memcpy(&(key->gp), &GLOBAL, sizeof(globalparam));
memcpy(&(key->gpdo), &GLOBAL, sizeof(globalparam));
simu_type *simu = &(key->simu);
simu->stop = simu->active = false;
simu->expected_number_simu = 0;
LOC_NULL(&(key->loc));
// trend left as is
key->cov = NULL;
key->meth= NULL;
}
void TREND_DELETE(trend_type *trend)
{
if (trend->TrendFunction!=NULL) free(trend->TrendFunction);
if (trend->LinearTrend!=NULL) free(trend->LinearTrend);
// TREND_NULL(trend);
}
void TREND_NULL(trend_type *trend)
{
trend->mean = 0.0;
trend->TrendModus = -1;
trend->lTrendFct = 0;
trend->lLinTrend = 0;
trend->TrendFunction = NULL;
trend->LinearTrend = NULL;
}
void LOC_NULL(location_type *loc) {
int d;
loc->spatialdim = loc->timespacedim = -1;
loc->totalpoints = loc->spatialtotalpoints = 0;
for (d=0; d<MAXSIMUDIM; d++) {
loc->xgr[d]=NULL;
loc->length[d] = -1;
}
loc->x = NULL;
loc->T[0] = loc->T[1] = loc->T[2] = 0.0;
}
void LOC_DELETE(location_type *loc) {
if (loc->x!=NULL) {
free(loc->x);
loc->x = NULL;
} else
if (loc->xgr[0] != NULL) {
free(loc->xgr[0]);
loc->xgr[0] = NULL;
}
}
void METHOD_DELETE(method_type **Meth) {
method_type *meth = *Meth;
int i;
if (meth==NULL) {
return;
}
for (i=0; i<MAXSUB; i++) {
if (meth->sub[i] != NULL) METHOD_DELETE(meth->sub + i);
}
if (meth->destruct!=NULL) meth->destruct(&(meth->S));
if (meth->caniso != NULL) free(meth->caniso);
if (meth->cproj != NULL) free(meth->cproj);
if (meth->loc != NULL) {
if (meth->loc->Time) {
if (meth->space != NULL) free(meth->space);
if (meth->sptime != NULL) free(meth->sptime);
} else {
assert(meth->space == meth->sptime);
if (meth->space != NULL) {
free(meth->space);
}
}
}
free(*Meth);
*Meth = NULL;
}
void METHOD_NULL( method_type *meth) {
int i;
meth->gp = meth->gpdo = NULL; // never leave unset -- used in METHOD_SET
meth->simu = NULL;
meth->loc = NULL;
// method
meth->nr = -Forbidden-1;
meth->compatible = false;
meth->nsub = 0;
for (i=0; i<MAXSUB; i++) meth->sub[i] = NULL;
meth->domethod = NULL;
meth->destruct = NULL;
meth->S = NULL;
meth->cov = NULL;
meth->caniso = NULL;
meth->cscale = 1.0;
meth->cvar = -1.0;
meth->cproj = NULL;
meth->xdimout = -1;
meth->type = TypeAny;
meth->hanging = NULL;
meth->space = meth->sptime = NULL;
// grani
}
// Achtung! removeonly muss zwingend mit dem originalen SUB + i aufgerufen
// werden. Der **Cov ueberschrieben wird. Ansonsten muessen die
// Kommentarzeilen geloescht werden !
void removeOnly(cov_model **Cov) {
cov_model *cov = *Cov,
*next = cov->sub[0];
if (cov->calling == NULL) next->calling = NULL;
else {
cov_model *calling = cov->calling;
// for (i=0; i<MAXSUB; i++) if (calling->sub[i] == cov) break;
// assert (i<MAXSUB);
// calling->sub[i] = cov->sub[0];
next->calling = calling;
}
*Cov = next;
COV_DELETE_WITHOUTSUB(&cov);
}
void COV_DELETE_WITHOUTSUB(cov_model **Cov) {
cov_model *cov = *Cov;
int i, j,
last = (cov->nr < 0) ? MAXPARAM : CovList[cov->nr].kappas;
if (cov->q != NULL) {
free(cov->q);
cov->q = NULL;
cov->qlen = 0;
}
if (cov->MLE != NULL) free(cov->MLE);
// print("last = %d %s\n", last, CovList[cov->nr].name);
for (i=0; i<last; i++) {
if (cov->p[i] != NULL) {
if (CovList[cov->nr].kappatype[i] >= LISTOF) {
listoftype *list = (listoftype*) (cov->p[i]);
for (j=0; j<cov->nrow[i]; j++) {
free(list->p[j]);
}
}
free(cov->p[i]);
cov->p[i] = NULL;
cov->ncol[i] = cov->nrow[i] = 0;
}
}
// important check in combination with above; can be easily removed or
// generalised !!!!
for (; i<MAXPARAM; i++) {
if (cov->p[i] != NULL) {
// print("last=%d current=%d max=%d\n", last, i, MAXPARAM);
// PrintModelInfo(*Cov); //
// assert(cov->p[i] == NULL);
}
}
free(*Cov);
*Cov = NULL;
}
void COV_DELETE(cov_model **Cov) {
cov_model *cov = *Cov;
int i;
if (cov == NULL) {
warning("*cov is NULL in COV_DELETE");
return;
}
for (i=0; i<MAXSUB; i++) { // cov->sub[i] // seit 1.10.07: luecken erlaubt
// bei PSgen !
if (cov->sub[i] != NULL) COV_DELETE(cov->sub + i);
}
// for (; i<MAXSUB; i++) assert(cov->sub[i] == NULL);
COV_DELETE_WITHOUTSUB(Cov);
}
void COV_NULL(cov_model *cov) {
int i;
for (i=0; i<MAXPARAM; i++) {
cov->p[i] = NULL;
cov->ncol[i] = cov->nrow[i] = 0;
}
cov->q = NULL;
cov->qlen = 0;
for (i=0; i<MAXSUB; i++) cov->sub[i] = NULL;
cov->nsub = 0;
cov->calling = NULL;
cov->tsdim = cov->xdim =cov->maxdim = UNSET;
cov->normalmix = cov->finiterange = cov->diag = cov->semiseparatelast
= cov->separatelast = false;
for (i=0; i<=Nothing; i++) cov->user[i] = cov->pref[i] = PREF_BEST;
// mpp und extrG muessen auf NONE sein !
for (; i<=Forbidden; i++) cov->user[i] = cov->pref[i] = PREF_NONE;
cov->derivatives = -1;
cov->tbm2num = false;
cov->spec.density = NULL;
cov->statIn = STAT_MISMATCH;
cov->isoIn = ISO_MISMATCH;
cov->vdim = M_MISMATCH;
cov->nr = -1;
cov->MLE = NULL;
cov->anyNAdown = cov->anyNAscaleup = TriBeyond;
cov->manipulating_x=false;
cov->spec.nmetro = -1;
cov->spec.sigma = -1.0;
cov->hess = false;
}
void addModel(cov_model **pcov, int covnr) {
cov_model *cov;
int i;
if (covnr >= GATTER && covnr <= LASTGATTER &&
*pcov != NULL && (*pcov)->nr == GATTER)
return; // keine 2 Gatter hintereinander
if ((cov = (cov_model*) malloc(sizeof(cov_model)))==0)
ERR("memory allocation error");
COV_NULL(cov);
cov->nr = covnr;
cov->nsub = 1;
if (*pcov!=NULL) {
cov->calling = (*pcov)->calling;
(*pcov)->calling = cov;
}
cov->sub[0] = *pcov;
*pcov = cov;
cov->tbm2num = NOT_IMPLEMENTED;
if (covnr == GATTER && *pcov != NULL)
for (i=0; i<=Forbidden; i++) {
cov->user[i] = cov->sub[0]->user[i];
cov->pref[i] = cov->sub[0]->pref[i];
}
}
int loc_set(double *x, double *T,
int spatialdim, /* spatial dim only ! */
int lx, bool Time, bool grid,
location_type *loc, globalparam *gp) {
int d, i, j;
unsigned long totalBytes;
// preference lists, distinguished by grid==true/false and dimension
// lists must end with Nothing!
loc->timespacedim = spatialdim + (int) Time;
loc->spatialdim = spatialdim;
loc->lx = lx;
if (spatialdim<1 || loc->timespacedim>MAXSIMUDIM) return ERRORDIM;
loc->grid= grid;
// coordinates
assert(loc->xgr[0]==NULL);
assert(loc->x == NULL);
totalBytes = sizeof(double) * lx * spatialdim;
if (loc->grid) {
if ((loc->xgr[0]=(double*) malloc(totalBytes))==NULL){
return ERRORMEMORYALLOCATION;
}
memcpy(loc->xgr[0], x, totalBytes);
for (d=1; d<spatialdim; d++) loc->xgr[d]= &(loc->xgr[0][d * lx]);
for (; d<MAXSIMUDIM; d++) loc->xgr[d]=NULL;
} else {
if ((loc->x=(double*) malloc(totalBytes))==NULL){
return ERRORMEMORYALLOCATION;
}
for (d=0; d<spatialdim; d++) loc->xgr[d] = x + d * lx;
for (j=i=0; i<lx; i++) {
for (d=0; d<spatialdim; d++) {
loc->x[j++] = loc->xgr[d][i];
}
}
for (d=0; d<MAXSIMUDIM; d++) loc->xgr[d]=NULL;
}
//Time
if ((loc->Time = Time)) {
memcpy(loc->T, T, sizeof(double) * 3);
//if (loc->grid)
loc->xgr[loc->spatialdim] = loc->T;
}
if (loc->grid) {
if ((lx!=3) ||
(InternalGetGridSize(loc->xgr, loc->timespacedim, loc->length))) {
PRINTF("%d\n", lx);
ERR("problem with the coordiates");
return ERRORCOORDINATES;
}
for (loc->spatialtotalpoints=1, d=0; d<loc->spatialdim; d++) {
loc->spatialtotalpoints *= loc->length[d];
}
loc->totalpoints = loc->spatialtotalpoints;
if (loc->Time) loc->totalpoints *= loc->length[loc->spatialdim];
} else { // not grid
loc->totalpoints = loc->spatialtotalpoints = loc->length[0]= lx;
if (loc->Time) {
double *Tx[MAXSIMUDIM];
Tx[0] = loc->T;
if (InternalGetGridSize(Tx, 1 , loc->length + loc->spatialdim)) {
return ERRORCOORDINATES;
}
loc->totalpoints *= loc->length[loc->spatialdim];
}
// just to say that considering these values does not make sense
for (d=1; d<loc->spatialdim; d++) loc->length[d]=0;
// correct value for higher dimensions (never used)
}
for (d=loc->timespacedim; d<MAXSIMUDIM; d++)
loc->length[d] = (int) RF_NAN; // 1
return NOERROR;
}
void Aniso(double *x, double *aniso, int origdim, int dim, double *newx) {
double dummy;
int i, j, k;
// hier kann die Geschwindigkeit verbessert werden,
// indem die Typen unterschieden werden
if (aniso == NULL) {
assert(dim == origdim);
memcpy(newx, x, sizeof(double) * origdim);
} else {
for (k=i=0; i<dim; i++) {
dummy = 0.0;
for (j=0; j<origdim; j++) {
dummy += x[j] * aniso[k++];
}
newx[i] = dummy;
}
}
}
void analyse_matrix(double *aniso, int row, int col,
bool *diag, bool *quasidiag, int *idx,
bool *semiseparatelast,
bool *separatelast) {
// diag, falls durch einfuegen von spalten diag-Matrix erhalten
// werden kann
// see also Type -- can it be unified ????
/*
-> i
| * 0 0
j 0 0 *
0 * 0
*/
bool notquasidiag=true, *taken;
int j, k, startidx, i, last;
if (aniso == NULL) {
*diag = *quasidiag = *separatelast = *semiseparatelast;
for (i=0; i<col; i++) idx[i] = i;
return;
}
taken = (bool *) malloc(row * sizeof(bool));
for (j=0; j<row; j++) {
taken[j]=false;
idx[j] = -1;
}
for (k=startidx=i=0; i<col; i++) {
for (j=0; j<row; j++, k++) if (aniso[k] != 0.0) break;
if (j < row) {
if ((notquasidiag = taken[j])) break;
taken[j] = true;
idx[j] = i;
for (j++, k++ ; j<row; j++) {
if ((notquasidiag = aniso[k++] != 0.0)) break;
}
}
if (notquasidiag) break;
}
if ((*diag = *quasidiag = !notquasidiag)) {
if (idx[0] == -1) idx[0] = 0;
for (j=1; j<row; j++) {
if (idx[j] <= idx[j-1]) {
if (idx[j] == -1) idx[j] = idx[j-1] + 1; else break;
}
}
*diag = j >= row;
}
if (!(*semiseparatelast = *diag)) {
last = col * row - 1;
for (k=last - row + 1; k<last; k++)
if (!(*separatelast = aniso[k++] == 0.0)) break;
}
if (!(*separatelast = *semiseparatelast)) {
// entered here only if last is set!
for (k=row - 1; k<last; k+=row)
if (!(*separatelast = aniso[k++] == 0.0)) break;
}
free(taken);
}
void addkappa(int i, const char *n, SEXPTYPE t) {
cov_fct *C = CovList + currentNrCov - 1;
assert(n[0] != '\0' &&
n[0] != 'k' && // reserved for standard parameter names
n[0] != 'u' && // reserved for standard submodel names
(n[0] != 'm' || n[1] != 'e') // reserved for me-thod
);
assert(i < MAXPARAM);
strcopyN(C->kappanames[i], n, PARAMMAXCHAR);
C->kappatype[i] = t;
}
void kappanames(const char* n1, SEXPTYPE t1) {
addkappa(0, n1, t1);
}
void kappanames(const char* n1, SEXPTYPE t1, const char* n2, SEXPTYPE t2) {
addkappa(0, n1, t1);
addkappa(1, n2, t2);
}
void kappanames(const char* n1, SEXPTYPE t1, const char* n2, SEXPTYPE t2,
const char* n3, SEXPTYPE t3) {
addkappa(0, n1, t1);
addkappa(1, n2, t2);
addkappa(2, n3, t3);
}
void kappanames(const char* n1, SEXPTYPE t1, const char* n2, SEXPTYPE t2,
const char* n3, SEXPTYPE t3, const char* n4, SEXPTYPE t4) {
addkappa(0, n1, t1);
addkappa(1, n2, t2);
addkappa(2, n3, t3);
addkappa(3, n4, t4);
}
void kappanames(const char* n1, SEXPTYPE t1, const char* n2, SEXPTYPE t2,
const char* n3, SEXPTYPE t3, const char* n4, SEXPTYPE t4,
const char* n5, SEXPTYPE t5) {
addkappa(0, n1, t1);
addkappa(1, n2, t2);
addkappa(2, n3, t3);
addkappa(3, n4, t4);
addkappa(4, n5, t5);
}
void kappanames(const char* n1, SEXPTYPE t1, const char* n2, SEXPTYPE t2,
const char* n3, SEXPTYPE t3, const char* n4, SEXPTYPE t4,
const char* n5, SEXPTYPE t5, const char* n6, SEXPTYPE t6) {
addkappa(0, n1, t1);
addkappa(1, n2, t2);
addkappa(2, n3, t3);
addkappa(3, n4, t4);
addkappa(4, n5, t5);
addkappa(5, n6, t6);
}
/*
real 0m47.752s
user 0m47.479s
sys 0m0.204s
make: *** [vrun] Error 1
real 0m47.996s
user 0m47.559s
sys 0m0.364s
real 4m18.616s
user 4m12.044s
sys 0m2.284s
make: *** [vrun] Error 1
real 4m25.319s
user 4m17.472s
sys 0m3.008s
schlather@DD787F2J:~/R/RF> time make v
echo -e "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"
*/
void addsub(int i, const char *n) {
cov_fct *C = CovList + currentNrCov - 1;
assert(n[0] != 'k' && n[0] != 'u');
assert(i < MAXSUB);
strcopyN(C->subnames[i], n, PARAMMAXCHAR);
}
void subnames(const char* n1) {
addsub(0, n1);
}
void subnames(const char* n1, const char* n2) {
addsub(0, n1);
addsub(1, n2);
}
void subnames(const char* n1, const char* n2, const char* n3) {
addsub(0, n1);
addsub(1, n2);
addsub(2, n3);
}
void subnames(const char* n1, const char* n2, const char* n3, const char* n4) {
addsub(0, n1);
addsub(1, n2);
addsub(2, n3);
addsub(3, n4);
}
void subnames(const char* n1, const char* n2, const char* n3, const char* n4,
const char* n5) {
addsub(0, n1);
addsub(1, n2);
addsub(2, n3);
addsub(3, n4);
addsub(4, n5);
}
void subnames(const char* n1, const char* n2, const char* n3, const char* n4,
const char* n5, const char* n6) {
addsub(0, n1);
addsub(1, n2);
addsub(2, n3);
addsub(3, n4);
addsub(4, n5);
addsub(5, n6);
}
void ErrCov(double *x, cov_model *cov, double *v) {
// int *f; PRINTF("%d\n", f[0]);
error("unallowed or undefined call of function");
}
void ErrCovNonstat(double *x, double *y, cov_model *cov, double *v) {
error("unallowed or undefined call of non-stationary function");
}
char InternalName[]="-";
void kappasize1(int i, cov_model *cov, int *nrow, int *ncol) {
*nrow = *ncol = 1;
}
void createmodel(const char *name, int kappas, size_fct kappasize,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, init_meth initmethod,
do_meth domethod, int vdim, pref_shorttype pref) {
int i;
if (currentNrCov==-1) InitModelList();
assert(CovList!=NULL);
assert(currentNrCov>=0);
if (currentNrCov>=MAXNRCOVFCTS) {
PRINTF("Error. Model list full.\n");
assert(false);
}
//printf("%d\n", PL);
if (PL > 8) PRINTF("%d %s vdim=%d\n", currentNrCov, name, vdim);
strcopyN(CovList[currentNrCov].name, name, MAXCHAR);
strcopyN(CovNames[currentNrCov], name, MAXCHAR);
assert(strcmp(InternalName, name));
CovList[currentNrCov].name[MAXCHAR-1]='\0';
if (strlen(name)>=MAXCHAR) {
PRINTF("Warning! Covariance name is truncated to `%s'",
CovList[currentNrCov].name);
}
cov_fct *C = CovList + currentNrCov;
assert(kappas >= 0);
C->kappas = kappas;
C->minsub = C->maxsub = 0;
C->primitive = true;
memcpy(C->pref, pref, sizeof(pref_shorttype));
C->stationary = stationary;
C->isotropy = isotropy;
C->vdim = vdim;
for (i=0; i<kappas; i++) {
sprintf(C->kappanames[i], "k%d", i); // default (repeated twice)
C->kappatype[i] = REALSXP;
}
C->kappasize = (kappasize == NULL) ? kappasize1 : kappasize;
C->range= range;
assert(range != NULL);
C->check = check;
assert(C->check != NULL);
for (i=0; i<SimulationTypeLength; i++)
C->implemented[i] = NOT_IMPLEMENTED;
C->internal = false;
C->naturalscale=NULL;
C->invSq = C->cov = C->D = C->D2 = C->tbm2 = C->D3 = C->D4 = ErrCov;
C->derivs = -1;
C->nonstat_cov = ErrCovNonstat;
C->nabla = C->hess = ErrCov;
C->coinit = C->ieinit = NULL;
C->alternative = NULL;
C->spectral=NULL;
C-> initspectral=NULL;
C->mppinit = NULL;
C->randomcoin = NULL;
C->mppget = NULL;
C->mppgetstat = NULL;
C->mpplocations = MPP_MISMATCH;
C->avef = NULL;
C->avelogg = NULL;
C->mineigenvalue = NULL;
C->drawlogmix = NULL;
C->logmixweight = NULL;
C->hyperplane=NULL;
// C->initspecial=NULL;
// C->special=NULL;
C->initmethod=initmethod;
C->domethod=domethod;
currentNrCov++;
}
int IncludePrim(const char *name, int kappas,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range) {
createmodel(name, kappas, NULL, stationary, isotropy, check, range,
NULL, NULL, UNIVARIATE, PREF_ALL);
return currentNrCov - 1;
}
int IncludePrim(const char *name, int kappas, size_fct kappasize,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range) {
createmodel(name, kappas, kappasize, stationary, isotropy, check, range,
NULL, NULL, UNIVARIATE, PREF_ALL);
return currentNrCov - 1;
}
int IncludePrim(const char *name, int kappas,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, int vdim) {
createmodel(name, kappas, NULL, stationary, isotropy, check, range,
NULL, NULL, vdim, PREF_ALL);
return currentNrCov - 1;
}
int IncludePrim(const char *name, int kappas, size_fct kappasize,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, int vdim) {
createmodel(name, kappas, kappasize, stationary, isotropy, check, range,
NULL, NULL, vdim, PREF_ALL);
return currentNrCov - 1;
}
int IncludePrim(const char *name, int kappas,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, pref_type pref,
int vdim) {
createmodel(name, kappas, NULL, stationary, isotropy, check, range,
NULL, NULL, vdim, pref);
return currentNrCov - 1;
}
int IncludePrim(const char *name, int kappas, size_fct kappasize,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range,pref_type pref, int vdim) {
createmodel(name, kappas, kappasize, stationary, isotropy, check, range,
NULL, NULL, vdim, pref);
return currentNrCov - 1;
}
void make_internal() {
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
C->internal = true;
}
// extern ?!
int IncludeModel(const char *name, char minsub, char maxsub, int kappas,
size_fct kappasize,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, pref_type pref,
init_meth initmethod, do_meth domethod, bool internal,
int vdim) {
createmodel(name, kappas, kappasize, stationary, isotropy, check, range,
initmethod, domethod, vdim, pref);
// assert(maxsub > 0); // check deleted 25. nov 2008 due to nugget
assert(maxsub >= minsub && maxsub <= MAXSUB);
int i,
nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
C->minsub = minsub;
C->maxsub = maxsub;
C->primitive = false;
C->internal = internal;
for (i=0; i<maxsub; i++) {
sprintf(C->subnames[i], "u%d", i+1); // default (repeated twice)
}
return nr;
}
int IncludeModel(const char *name, char minsub, char maxsub, int kappas,
size_fct kappasize,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, pref_type pref) {
return
IncludeModel(name, minsub, maxsub, kappas, kappasize, stationary, isotropy,
check, range, pref, NULL, NULL, false, UNIVARIATE);
}
int IncludeModel(const char *name, char minsub, char maxsub, int kappas,
stationary_type stationary, isotropy_type isotropy,
checkfct check, rangefct range, pref_type pref) {
return
IncludeModel(name, minsub, maxsub, kappas, NULL, stationary, isotropy,
check, range, pref, NULL, NULL, false, UNIVARIATE);
}
void addCov(covfct cf, covfct D, covfct D2, natscalefct naturalscale) {
int nr = currentNrCov - 1;
assert((nr>=0) && (nr<currentNrCov) && cf!=NULL);
cov_fct *C = CovList + nr;
assert(cf != NULL);
C->cov = cf;
if (C->derivs < 0) C->derivs = 0;
assert(C->nonstat_cov == ErrCovNonstat);
assert(C->vdim == UNIVARIATE || C->vdim == SUBMODELM || D==NULL);
if (D != NULL) {
assert(cf != NULL);
if (C->cov!=NULL && C->derivs < 1) C->derivs = 1;
C->D=D;
assert(C->isotropy == ISOTROPIC || C->isotropy == SPACEISOTROPIC ||
C->isotropy == PREVMODELI);
C->implemented[TBM2] = NUM_APPROX;
C->implemented[TBM3] = true;
}
if (D2 != NULL) {
assert(D != NULL);
C->D2 = D2;
if (C->cov!=NULL && C->D != NULL && C->derivs < 2) C->derivs = 2;
}
if (C->maxsub > 0 && currentNrCov>GATTER+10)
assert(naturalscale == NULL);
C->naturalscale=naturalscale;
C->implemented[Direct] = cf != NULL;
C->implemented[CircEmbed] = cf != NULL &&
(C->stationary == STATIONARY || C->stationary == PREVMODELS);
C->implemented[Sequential] = C->implemented[CircEmbed] && C->vdim <= 1;
}
void addCov(covfct cf, covfct D, natscalefct naturalscale) {
addCov(cf, D, NULL, naturalscale);
}
void addCov(covfct cf, covfct D, covfct D2, covfct D3, covfct D4,
natscalefct naturalscale) {
addCov(cf, D, D2, naturalscale);
cov_fct *C = CovList + currentNrCov - 1;
C->D3 = D3;
C->D4 = D4;
assert(C->derivs == 2 && D3!=NULL && D4!=NULL);
C->derivs = 4;
}
void addCov(nonstat_covfct cf) {
int nr = currentNrCov - 1;
assert((nr>=0) && (nr<currentNrCov) && cf!=NULL);
cov_fct *C = CovList + nr;
assert(cf != NULL);
C->nonstat_cov = cf;
C->implemented[Direct] = true;
C->implemented[CircEmbed] =
C->stationary == STATIONARY || C->stationary == PREVMODELS;
C->implemented[Sequential] = C->implemented[CircEmbed] && C->vdim <= 1;
if (C->derivs < 0) C->derivs = 0;
}
void nablahess(covfct nabla, covfct hess) {
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert((nr>=0) && (nr<currentNrCov));
assert(C->vdim == UNIVARIATE || C->vdim == SUBMODELM || nabla==NULL);
assert(C->cov != NULL && nabla!=NULL && hess != NULL);
C->nabla=nabla;
C->hess = hess;
// if (C->derivs>0) {
// printf("%s: %d\n", C->name, C->derivs);
// assert(C->derivs >= 2);
// }
}
void addCov(aux_covfct auxcf)
{
int nr = currentNrCov - 1;
assert((nr>=0) && (nr<currentNrCov) && auxcf!=NULL);
cov_fct *C = CovList + nr;
assert(C->cov == ErrCov && C->nonstat_cov==ErrCovNonstat &&
(C->stationary == AUXMATRIX || C->stationary == AUXVECTOR || true));
C->aux_cov = auxcf;
if (C->derivs < 0) C->derivs = 0;
}
int addFurtherCov(covfct cf, covfct D, covfct D2) {
assert(currentNrCov > 0);
cov_fct *C = CovList + currentNrCov;
memcpy(C, C - 1, sizeof(cov_fct));
assert(C->vdim == UNIVARIATE || C->vdim == SUBMODELM || D == NULL);
strcopyN(CovNames[currentNrCov], InternalName, MAXCHAR);
C->name[0] = InternalName[0];
strcopyN(C->name + 1, CovList[currentNrCov-1].name, MAXCHAR - 1);
if (cf != NULL) {
C->cov = cf;
C->derivs = 0;
}
if (D != NULL) {
assert(cf != NULL);
C->D = D;
C->derivs = 1;
// printf("%s %d\n", C->name, C->isotropy);
assert(C->isotropy == ISOTROPIC || C->isotropy == SPACEISOTROPIC ||
C->isotropy == ZEROSPACEISO ||
C->isotropy == PREVMODELI);
C->implemented[TBM2] = NUM_APPROX;
C->implemented[TBM3] = true;
}
if (D2 != NULL) {
assert(D != NULL);
C->D2 = D2;
C->derivs = 2;
}
C->internal = true; // addCov is used without previous call of IncludeModel
currentNrCov++;
return currentNrCov - 1;
}
int addFurtherCov(covfct cf, covfct D) {
return addFurtherCov(cf, D, NULL);
}
int addFurtherCov(nonstat_covfct cf) {
assert(currentNrCov > 0);
cov_fct *C = CovList + currentNrCov;
memcpy(C, C - 1, sizeof(cov_fct));
strcopyN(CovNames[currentNrCov], InternalName, MAXCHAR);
C->name[0] = InternalName[0];
strcopyN(C->name + 1, CovList[currentNrCov-1].name, MAXCHAR - 1);
C->derivs = -1;
if (cf != NULL) {
C->nonstat_cov = cf;
C->derivs = 0;
}
C->D = ErrCov;
C->internal = true; // addCov is used without previous call of IncludeModel
C->mppget = NULL;
C->mppgetstat = NULL;
currentNrCov++;
return currentNrCov - 1;
}
int addFurtherCov(nonstat_covfct cf, covfct D, covfct D2) {
assert(currentNrCov > 0);
cov_fct *C = CovList + currentNrCov;
memcpy(C, C - 1, sizeof(cov_fct));
strcopyN(CovNames[currentNrCov], InternalName, MAXCHAR);
C->name[0] = InternalName[0];
strcopyN(C->name + 1, CovList[currentNrCov-1].name, MAXCHAR - 1);
C->derivs = -1;
if (cf != NULL) {
C->nonstat_cov = cf;
C->derivs = 0;
}
C->D = ErrCov;
C->internal = true; // addCov is used without previous call of IncludeModel
C->mppget = NULL;
C->mppgetstat = NULL;
currentNrCov++;
return currentNrCov - 1;
}
void addLocal(getlocalparam coinit, getlocalparam ieinit) {
int nr = currentNrCov - 1;
assert(nr>=0 && nr < currentNrCov) ;
cov_fct *C = CovList + nr;
assert(C->D!=ErrCov);
if ((C->implemented[CircEmbedIntrinsic] = ieinit != NULL)) {
assert(C->D2 != NULL);
C->ieinit = ieinit;
}
if ((C->implemented[CircEmbedCutoff] = coinit != NULL)) {
C->coinit = coinit;
}
}
void addCallLocal(altlocalparam alt) {
int nr = currentNrCov - 1;
assert(nr>=0 && nr < currentNrCov) ;
cov_fct *C = CovList + nr;
C->alternative = alt;
}
int addTBM(covfct tbm2) {
// must be called always AFTER addCov !!!!
int nr = currentNrCov - 1;
assert((nr>=0) && (nr<currentNrCov));
cov_fct *C = CovList + nr;
assert(C->stationary == STATIONARY || C->stationary == PREVMODELS);
assert(C->vdim == UNIVARIATE || C->vdim == SUBMODELM);
C->tbm2=tbm2;
if (tbm2 != NULL) {
// addTBM is called from the other addTBM's -- so tbm2 might
// be NULL
assert(C->isotropy==ISOTROPIC || C->isotropy==SPACEISOTROPIC ||
C->isotropy==PREVMODELI);
assert(C->D != ErrCov);
C->implemented[TBM2] = IMPLEMENTED;
}
// IMPLEMENTED must imply the NUM_APPROX to simplify the choice
// between TBM2 and Tbm2Num
return nr;
}
void addTBM(covfct tbm2, spectral_init initspectral, spectral_do spectral) {
int nr = addTBM(tbm2);
cov_fct *C = CovList + nr;
assert(C->vdim == UNIVARIATE || C->vdim == SUBMODELM);
C->initspectral=initspectral;
C->spectral=spectral;
C->implemented[SpectralTBM] = true;
assert(C->stationary == STATIONARY || C->stationary == PREVMODELS);
}
void addTBM(spectral_init initspectral, spectral_do spectral) {
addTBM((covfct) NULL, initspectral, spectral);
}
void addInv(covfct invSq) {
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert((nr>=0) && (nr<currentNrCov));
// printf("%s %d %d\n", C->name, invSq, invstableSq);
C->invSq=invSq;
}
void addHyper(hyper_pp_fct hyper_pp) {
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert((nr>=0) && (nr<currentNrCov));
C->hyperplane=hyper_pp;
C->implemented[Hyperplane] = hyper_pp!=NULL;
}
//void addSpecialMeth(initstandard initspecial, dometh special) {
/// int nr = currentNrCov - 1;
// cov_fct *C = CovList + nr;
// C->initspecial=initspecial;
// C->special=special;
// if ((special!=NULL) || (initspecial!=NULL))
// assert((special!=NULL) && (initspecial!=NULL));
// C->implemented[Special] = true;
//}
void addMarkov(int *variable) {
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
C->implemented[Markov] = true;
*variable = nr;
}
void addMPP(mpp_init mppinit, mpp_model randomcoin, mpp_get mppget,
sd_fct sd, int loc, bool timesep){
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
// print("addmpp %s\n", C->name);
assert(mppinit!=NULL && randomcoin!=NULL && mppget!=NULL);
C->mppinit = mppinit;
C->randomcoin = randomcoin;
C->mppget = mppget;
assert(C->stationary >= COVARIANCE);
C->sd = sd;
C->mpplocations = loc;
C->implemented[RandomCoin] = true;
}
void addMPP(mpp_init mppinit, mpp_model randomcoin, mpp_get mppget,
sd_fct sd, int loc){
addMPP(mppinit, randomcoin, mppget, sd, loc, false);
}
void addMPP(mpp_init mppinit, mpp_model randomcoin, mpp_get mppget,
sd_fct sd, int loc,
ave_fct avef, ave_logfct avelogg) {
addMPP(mppinit, randomcoin, mppget, sd, loc, true);
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert(avef != NULL && avelogg != NULL);
assert(C->stationary == STATIONARY);
C->avef = avef;
C->avelogg = avelogg;
C->implemented[Average] = true;
}
void addMPP(mpp_init mppinit, mpp_model randomcoin, mpp_getstat mppgetstat,
sd_fct sd, int loc, bool timesep){
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert(mppinit!=NULL && randomcoin!=NULL && mppgetstat!=NULL);
C->mppinit = mppinit;
C->randomcoin = randomcoin;
C->mppgetstat = mppgetstat;
// printf("%d\n", C->stationary);
assert(C->stationary < COVARIANCE || C->stationary == PREVMODELS );
C->sd = sd;
C->mpplocations = loc;
C->implemented[RandomCoin] = true;
}
void addMPP(mpp_init mppinit, mpp_model randomcoin, mpp_getstat mppgetstat,
sd_fct sd, int loc){
addMPP(mppinit, randomcoin, mppgetstat, sd, loc, false);
}
void addMPP(mpp_init mppinit, mpp_model randomcoin, mpp_getstat mppgetstat,
sd_fct sd, int loc,
ave_fct avef, ave_logfct avelogg) {
addMPP(mppinit, randomcoin, mppgetstat, sd, loc, true);
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert(avef != NULL && avelogg != NULL);
assert(C->stationary == STATIONARY || C->stationary == PREVMODELS);
C->avef = avef;
C->avelogg = avelogg;
C->implemented[Average] = true;
}
void addSpecial(realfct mineigen){
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
C->mineigenvalue = mineigen;
}
void addGaussMixture(draw_logmix drawlogmix,
log_mixweight logmixweight) {
int nr = currentNrCov - 1;
cov_fct *C = CovList + nr;
assert(drawlogmix != NULL && logmixweight != NULL);
C->drawlogmix = drawlogmix;
C->logmixweight = logmixweight;
}
int getmodelnr(char *name) {
// == -1 if no matching function is found
// == -2 if multiple matching fctns are found, without one matching exactly
// == -3 if internal name is passed
// if more than one match exactly, the last one is taken (enables overwriting
// standard functions)
if (currentNrCov==-1) InitModelList();
if (!strcmp(name, InternalName)) return MATCHESINTERNAL;
return Match(name, CovNames, currentNrCov);
}
int Match(char *name, name_type List, int n) {
// == -1 if no matching name is found
// == -2 if multiple matching fctns are found, without one matching exactly
unsigned int ln;
int Nr;
Nr=0;
ln=strlen(name);
while ( Nr < n && strncmp(name, List[Nr], ln)) Nr++;
if (Nr < n) {
if (ln==strlen(List[Nr])) // exactmatching -- take first -- changed 1/7/07
return Nr;
// a matching function is found. Are there other functions that match?
int j;
bool multiplematching=false;
j=Nr+1; // if two or more covariance functions have the same name
// the last one is taken
while (j<n) {
while ( (j<n) && strncmp(name, List[j], ln)) {j++;}
if (j<n) {
if (ln==strlen(List[j])) { // exactmatching -- take first
return j;
}
else {multiplematching=true;}
}
j++;
}
if (multiplematching) {return MULTIPLEMATCHING;}
} else return NOMATCHING;
return Nr;
}
void UserGetNatScaling(double *natscale) {
GetNaturalScaling(STORED_MODEL[MODEL_USER], natscale);
}
void GetNaturalScaling(cov_model *cov, double *natscale)
{ // called also by R
// values of naturalscaling:
//#define NATSCALE_EXACT 1
//#define NATSCALE_APPROX 2
//#define NATSCALE_MLE 3 /* check fitvario when changing !! */
// +10 if numerical is allowed
if (cov->nr >= GATTER && cov->nr <= LASTGATTER) cov = cov->sub[0];
cov_fct *C = CovList + cov->nr;
*natscale = 0.0;
// printf("%s %d\n", C->name, C->maxsub);
if (C->maxsub!=0) XERR(ERRORFAILED);
if (C->isotropy != ISOTROPIC || C->stationary != STATIONARY ||
C->vdim != UNIVARIATE)
XERR(ERRORANISOTROPIC);
if (C->naturalscale!=NULL) {
*natscale = C->naturalscale(cov);
if (ISNAN(*natscale) || ISNA(*natscale) || *natscale != 0.0) {
return;
}
}
if (NS != NATSCALE_ORNUMERIC) XERR(ERRORRESCALING);
double x,newx,yold,y,newy;
covfct cf;
int wave, i;
if ((C->cov)==nugget) XERR(ERRORRESCALING);
if ((cf=C->cov)==NULL) XERR(ERRORNOTDEFINED);
// already calculated ?
// parami=KAPPA; // do not compare mean,variance, etc.
// if (oldcovnr==*covnr) {
// for (; parami<=LASTKAPPA; parami++) {
// if (oldp[parami]!=p[parami]) {
// break;
// }
// }
// if (parami > LASTKAPPA) {
// *natscale=OldNatScale;
// return;
// }
// }
// for (; parami<=LASTKAPPA; parami++) {oldp[parami]=p[parami];}
/* **************************************
Now, find a quick and good solution for NewInvScale --
*/
wave = 0;
x = 1.0;
cf(&x, cov, &yold);
if (ISNA(yold) || ISNAN(yold)) {*natscale = RF_NAN; return;}
if (yold > 0.05) {
double leftx;
x *= 2.0;
cf(&x, cov, &y);
while (y > 0.05) {
if (yold<y){
wave++;
if (wave>10) XERR(ERRORWAVING);
}
yold = y;
x *= 2.0;
if (x>1E30) XERR(ERRORRESCALING); // or use a separate ERROR
cf(&x, cov, &y);
}
leftx = x * 0.5;
for (i=0; i<3 /* good choice?? */ ;i++) {
if (y==yold) XERR(ERRORWAVING); // should never appear
newx = x + (x-leftx)/(y-yold)*(0.05-y);
cf(&newx, cov, &newy);
if (newy > 0.05) {
leftx = newx;
yold = newy;
} else {
x = newx;
y = newy;
}
}
if (y==yold) XERR(ERRORWAVING); // should never appear
*natscale = 1.0 / ( x + (x-leftx)/(y-yold)*(0.05-y) );
} else {
double rightx;
x *= 0.5;
cf(&x, cov, &y);
while (y < 0.05) {
if (yold>y){ wave++; if (wave>10) XERR(ERRORWAVING); }
yold = y;
x *= 0.5;
if (x<1E-30) XERR(ERRORRESCALING); //or use a separate ERROR
cf(&x, cov, &y);
}
rightx = x * 2.0;
for (i=0; i<3 /* good choice?? */ ;i++) {
if (y==yold) XERR(ERRORWAVING); // should never appear
newx = x + (x-rightx)/(y-yold)*(0.05-y);
cf(&x, cov, &newy);
if (newy < 0.05) {
rightx = newx;
yold = newy;
} else {
x = newx;
y = newy;
}
}
if (y==yold) XERR(ERRORWAVING); // should never appear
*natscale = 1.0 / ( x + (x-rightx)/(y-yold)*(0.05-y));
}
}
void Getxsimugr(coord_type x, double *aniso, int timespacedim, double **xsimugr) {
// bisher nur fuer Diagonalmatrizen
int n, i, w;
if (aniso == NULL) {
for(w=0; w<timespacedim; w++) {
for (i=0; i<3; i++) {
xsimugr[w][i] = x[w][i];
}
}
} else {
for(n=w=0; w<timespacedim; w++, n+=timespacedim+1) {
for (i=0; i<3; i++) {
xsimugr[w][i] = aniso[n] * x[w][i];
}
}
}
}
void GetCoordinates(int *keyNr, double *x, int *err)
{
key_type *key;
if (currentNrCov==-1) {*err=-1;goto ErrorHandling;}
if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {*err=-2; goto ErrorHandling;}
key = &(KEY[*keyNr]);
if (!key->simu.active) {*err=-3; goto ErrorHandling;}
if (key->trend.TrendModus==0) {*err=-4; goto ErrorHandling;}
*err = ERRORNOTPROGRAMMED; // TODO !!
return;
ErrorHandling:
return;
}
void covcpy(cov_model **localcov, cov_model *cov, bool insertgatter,
bool keepuser) {
assert(cov != NULL);
int i,
covnr = cov->nr,
n = -1;
cov_model *current;
cov_fct *C = CovList + cov->nr;
if (insertgatter && covnr >= GATTER && covnr <= LASTGATTER)
error("covcpy detects # at forbidden place -- please contact author");
if (*localcov != NULL) error("local cov not NULL");
if ((*localcov = (cov_model*) malloc(sizeof(cov_model)))==0)
ERR("memory allocation error");
if (insertgatter) { // covnr != DOLLAR &&
COV_NULL(*localcov);
(*localcov)->nr = GATTER;
(*localcov)->nsub = 1;
(*localcov)->isoIn = cov->isoIn;
(*localcov)->statIn = cov->statIn;
(*localcov)->vdim = cov->vdim;
if (((*localcov)->sub[0] = (cov_model*) malloc(sizeof(cov_model)))==0)
ERR("memory allocation error");
current = (*localcov)->sub[0];
} else current = *localcov;
memcpy(current, cov, sizeof(cov_model)); // replaces COV_NULL(*localcov);
if (!keepuser) {
for (i=0; i<Forbidden; i++)
current->user[i] = PREF_BEST;
}
current->calling = *localcov;
for (i=0; i<MAXPARAM; i++) {
if (cov->p[i] == NULL) {
continue;
}
if (C->kappatype[i] >= LISTOF) {
current->p[i] = (double*) malloc(sizeof(listoftype));
int j, len = cov->nrow[i];
listoftype *p = (listoftype*) (cov->p[i]),
*q = (listoftype*) (current->p[i]);
for (j=0; j<len; j++) {
if (C->kappatype[i]== LISTOF + REALSXP) {
n = sizeof(double);
} else assert(false);
n *= p->nrow[j] * p->ncol[j];
q->nrow[j] = p->nrow[j];
q->ncol[j] = p->ncol[j];
q->p[j] = (double*) malloc(n);
memcpy(q->p[i], p->p[i], n);
}
} else {
if (C->kappatype[i]==REALSXP) {
n = sizeof(double);
} else if (C->kappatype[i]==INTSXP) {
n = sizeof(int);
} else assert(false);
n *= cov->nrow[i] * cov->ncol[i];
current->p[i] = (double*) malloc(n);
memcpy(current->p[i], cov->p[i], n);
}
}
if (cov->q != NULL) {
n = sizeof(double) * current->qlen;
current->q = (double*) malloc(n);
memcpy(current->q, cov->q, n);
} else assert(current->qlen==0);
for (i=0; i<MAXSUB; i++) {
current -> sub[i] = NULL;
if (cov->sub[i] == NULL) continue;
if (insertgatter &&
cov->sub[i]->nr >= GATTER && cov->sub[i]->nr <= LASTGATTER)
covcpy(current->sub + i, cov->sub[i]->sub[0], insertgatter, keepuser);
else
covcpy(current->sub + i, cov->sub[i], insertgatter, keepuser);
current->sub[i]->calling = current;
}
}
//void covcpy(cov_model **localcov, cov_model *cov) {
// covcpy(localcov, cov, true, true);
//}
double *getAnisoMatrix(method_type *meth) {
double *ani;
int i, bytes,
origdim = meth->loc->timespacedim,
dimP1 = origdim + 1;
if (meth->caniso != NULL) {
bytes = sizeof(double) * origdim * meth->xdimout;
ani = (double *) malloc(bytes);
// PrintMethodInfo(meth);
// printf("%ld %ld %d %d %d\n", ani, meth->caniso, bytes,
// origdim, meth->xdimout); assert(false);
memcpy(ani, meth->caniso, bytes);
} else {
double a = 1.0 / meth->cscale;
if (meth->cproj == NULL) {
bytes = origdim * origdim;
ani = (double *) calloc(bytes, sizeof(double));
for (i=0; i<bytes; i+=dimP1) ani[i] = a;
} else {
bytes = origdim * meth->xdimout;
ani = (double *) calloc(bytes, sizeof(double));
for (i=0; i<meth->xdimout; i++) ani[i + meth->cproj[i] * origdim] = a;
}
}
return ani;
}
double *getAnisoMatrix(cov_model *cov, int dim) {
double *ani;
int i, bytes, dimout = dim,
dimP1 = dim + 1,
*proj = (int *)cov->p[DPROJ];
if (cov->p[DANISO] != NULL) {
bytes = sizeof(double) * dim * dimout;
ani = (double *) malloc(bytes);
memcpy(ani, cov->p[DANISO], bytes);
} else {
double a;
a = (cov->p[DSCALE] == NULL) ? 1.0 : 1.0 / cov->p[DSCALE][0];
if (proj == NULL) {
bytes = dim * dimout;
ani = (double *) calloc(bytes, sizeof(double));
for (i=0; i<bytes; i+=dimP1) ani[i] = a;
} else {
dimout = cov->nrow[DPROJ];
bytes = dim * dimout;
ani = (double *) calloc(bytes, sizeof(double));
for (i=0; i<dimout; i++) ani[i + proj[i] * dim] = a;
}
}
return ani;
}
double GetDiameter(double *origcenter, double *origmin, double *origmax,
double *aniso, int origdim, int spatialdim,
double *min, double *max, double *center) {
// calculates twice the distance between origcenter %*% aniso and
// all 2^spatialdim corners spanned by min and max; returns maximum distance.
// NOTE: origcenter is not alway (min + max)/2 since origcenter might be
// given by the user and not calculated automatically !
double *lx, *sx, distsq, dummy,
diameter=0.0;
bool *j;
int d;
j = (bool*) malloc( (origdim + 1) * sizeof(double));
lx = (double*) malloc( origdim * sizeof(double));
sx = (double*) malloc(spatialdim * sizeof(double));
Aniso(origcenter, aniso, origdim, spatialdim, center);
for (d=0; d<origdim; d++) {
j[d]=false;
lx[d]=origmin[d];
}
j[origdim] = false;
for (d=0; d<spatialdim; d++) {
min[d] = R_PosInf;
max[d] = R_NegInf;
}
while(true) {
d=0;
while(j[d]) {
lx[d]=origmin[d];
j[d++]=false;
}
if (d==origdim) break;
j[d]=true;
lx[d]=origmax[d];
Aniso(lx, aniso, origdim, spatialdim, sx);
// suche maximale Distanz von center zu transformierter Ecke
distsq = 0.0;
for (d=0; d<spatialdim; d++) {
if (min[d] > sx[d]) min[d] = sx[d];
if (max[d] < sx[d]) max[d] = sx[d];
dummy = center[d] - sx[d];
distsq += dummy * dummy;
// printf("%d sx=%f %f %f diam=%f\n", d, sx[d], min[d], max[d], distsq);
}
if (distsq > diameter) diameter = distsq;
}
free(j);
free(lx);
free(sx);
return 2.0 * sqrt(diameter);
}
double GetDiameter(double *origcenter, double *min, double *max,
double *aniso, int tsdim, int spatialdim) {
double diam,
*dummymin, *dummymax, *dummycenter;
dummymin = (double*) malloc(spatialdim * sizeof(double));
dummymax = (double*) malloc(spatialdim * sizeof(double));
dummycenter = (double*) malloc(spatialdim * sizeof(double));
diam = GetDiameter(origcenter, min, max,
aniso, tsdim, // meth->loc->timespacedim,
spatialdim,
dummymin, dummymax, dummycenter);
free(dummymin);
free(dummymax);
free(dummycenter);
return diam;
}
void GetMinMax(method_type *meth, double *min, double *max, double *center,
int MaxDim) {
// untersucht die Original-Werte fuer Koordinaten fuer center!!
// aber gibt diameter in den neuen Koordinaten wieder.
location_type *loc = meth->loc;
int d,
origdim = loc->timespacedim,
spatialdim = loc->spatialdim;
if (loc->grid) {
// diameter of the grid
// printf("origdim =%d\n", origdim);
for (d=0; d<origdim; d++) {
if (loc->xgr[d][XEND] > loc->xgr[d][XSTART]) {
min[d] = loc->xgr[d][XSTART];
max[d] = loc->xgr[d][XSTART] + loc->xgr[d][XSTEP] *
(double) (loc->length[d] - 1);
} else {
min[d] = loc->xgr[d][XSTART] + loc->xgr[d][XSTEP] *
(double) (loc->length[d] - 1);
max[d] = loc->xgr[d][XSTART];
}
}
// assert(!false);
} else { // not simugrid
double *xx=loc->x;
int i,
endfor = loc->length[0] * loc->timespacedim;
for (d=0; d<MaxDim; d++) {min[d]=RF_INF; max[d]=RF_NEGINF;}
// to determine the diameter of the grid determine as approxmation
// componentwise min and max corner
for (i=0; i<endfor; ) {
for (d=0; d<spatialdim; d++, i++) {
//temporal part need not be considered, but for ease included
if (xx[i]<min[d]) min[d] = xx[i];
if (xx[i]>max[d]) max[d] = xx[i];
}
}
if (loc->Time) {
assert(d == origdim - 1);
if (loc->T[XEND] > loc->T[XSTART]) {
min[d] = loc->T[XSTART];
max[d] = loc->T[XSTART] + loc->T[XSTEP] *
(double) (loc->length[d] - 1);
} else {
min[d] = loc->T[XSTART] + loc->T[XSTEP] * (double) (loc->length[d] - 1);
max[d] = loc->T[XSTART];
}
}
} // !simugrid
for (d=0; d<origdim; d++) { // note here no time component!
center[d] = 0.5 * (max[d] + min[d]);
// printf("%d center %f\n", d, center[d]);
}
for ( ; d<MaxDim; d++) {
min[d] = 0.0;
max[d] = 0.0;
center[d] = 0.0;
}
// assert(false);
}
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 /= 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;
}
#define F_NUMBERS1 3
#define F_NUMBERS2 3
bool HOMEMADE_NICEFFT=false;
unsigned long NiceFFTNumber(unsigned long n) {
int f[F_NUMBERS1]={2,3,5};
unsigned long i,ii,j,jj,l,ll,min, m=1;
if (HOMEMADE_NICEFFT) {
if (n<=1) return n;
for (i=0; i<F_NUMBERS1; 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
int f[F_NUMBERS2]={2,3,5};
return nextn(n, f, F_NUMBERS2);
}
}
void cpyMethod(method_type *meth, method_type *ziel, bool cp_aniso) {
int xdimout = meth->xdimout,
tsdim = meth->loc->timespacedim;
ziel->gp = meth->gp;
ziel->gpdo = meth->gpdo;
ziel->simu = meth->simu;
ziel->loc = meth->loc;
// nr (by METHOD_NULL)
// compatible (by METHOD_NULL)
// nsub (by METHOD_NULL)
// domethod (by METHOD_NULL)
// destruct (by METHOD_NULL)
// S (by METHOD_NULL)
// cov (by METHOD_NULL)
if (cp_aniso) {
if (meth->cproj != NULL) {
ziel->cproj = (int*) malloc(xdimout * sizeof(int));
assert(ziel->cproj != NULL);
memcpy(ziel->cproj, meth->cproj, xdimout * sizeof(int));
}
if (meth->caniso != NULL) {
ziel->caniso = (double*) malloc(xdimout * tsdim * sizeof(double));
assert(ziel->caniso != NULL);
memcpy(ziel->caniso, meth->caniso, xdimout * tsdim * sizeof(double));
}
}
ziel->cscale = meth->cscale;
ziel->cvar = meth->cvar;
// cproj (by METHOD_NULL)
ziel->xdimout = meth->xdimout;
ziel->type = meth->type;
// hanging (by METHOD_NULL)
// space (by METHOD_NULL)
// sptime (by METHOD_NULL)
// grani (not set !)
}
void expandgrid(coord_type xgr, int *len, double **xx, int nrow){
double *x,* y; /* current point within grid, but without
anisotropy transformation */
int pts, w, k, total, i, d, dimM1 = nrow - 1,
*yi, /* counter for the current position in the grid */
ncol = nrow;
for (pts=1, i=0; i<nrow; i++) {
pts *= len[i];
}
y = (double*) malloc(ncol * sizeof(double));
yi = (int*) malloc(ncol * sizeof(int));
total = ncol * pts;
x = *xx = (double*) malloc(sizeof(double) * total);
for (w=0; w<ncol; w++) {y[w]=xgr[w][XSTART]; yi[w]=0;}
for (k=0; k<total; ){
for (d=0; d<ncol; d++, k++) {
x[k] = y[d];
}
i = 0;
(yi[i])++;
y[i] += xgr[i][XSTEP];
while(yi[i]>=len[i]) {
yi[i]=0;
y[i] = xgr[i][XSTART];
if (i<dimM1) {
i++;
(yi[i])++;
y[i] += xgr[i][XSTEP];
} else {
assert(k==total);
}
}
}
free(y);
free(yi);
}
void expandgrid(coord_type xgr, int *len, double **xx, double* aniso,
int nrow, int ncol){
double *x, * y; /* current point within grid, but without
anisotropy transformation */
int pts, w, k, total, n, i, d,
*yi, /* counter for the current position in the grid */
dimM1 = nrow - 1;
for (pts=1, i=0; i<nrow; i++) {
pts *= len[i];
}
total = ncol * pts;
x = *xx = (double*) malloc(sizeof(double) * total);
y = (double*) malloc(ncol * sizeof(double));
yi = (int*) malloc(ncol * sizeof(int));
for (w=0; w<nrow; w++) {y[w]=xgr[w][XSTART]; yi[w]=0;}
for (k=0; k<total; ){
if (aniso==NULL) {
for (d=0; d<ncol; d++) x[d] = y[d];
} else {
for (n=d=0; d<ncol; d++, k++) {
x[k] = 0.0;
for(w=0; w<nrow; w++) {
x[k] += aniso[n++] * y[w];
}
}
}
i = 0;
(yi[i])++;
y[i] += xgr[i][XSTEP];
while(yi[i]>=len[i]) {
yi[i]=0;
y[i] = xgr[i][XSTART];
if (i<dimM1) {
i++;
(yi[i])++;
y[i] += xgr[i][XSTEP];
} else {
assert(k==total);
}
}
}
free(y);
free(yi);
}
void grid2grid(coord_type xgr, static_grid grani, double *aniso, int dim) {
// diag is assumed
double factor;
int n, d, i;
if (aniso == NULL) {
for(d=0; d<dim; d++) {
for (i=0; i<3; i++) {
grani[d][i] = xgr[d][i];
}
}
} else {
for(n=d=0; d<dim; d++, n+=dim+1) {
factor = aniso[n];
// factor = d < cov->nrow ? aniso[n] : 0.0; // hier zwingend diag !!
for (i=0; i<3; i++) {
grani[d][i] = xgr[d][i] * factor;
}
}
}
}
void xtime2x(double *x, int nx, double *T, int timelen,
double **newx, int timespacedim) {
double *y, // umhaengen zwingend notwendig, da u.U. **xx und *newx
// der gleiche aufrufende Pointer ist, d.h. es wird "ueberschrieben"
*z, t;
int j, k, i, d,
spatialdim = timespacedim - 1;
z = *newx = (double*) malloc(sizeof(double) * timespacedim * nx * timelen);
for (k=j=0, t=T[XSTART]; j<timelen; j++, t += T[XSTEP]){
y = x;
for (i=0; i<nx; i++) {
for (d=0; d<spatialdim; d++, y++) {
z[k++] = *y;
}
z[k++] = t;
}
}
}
void xtime2x(double *x, int nx, double *T, int len, double **newx,
double *aniso, int nrow, int ncol) {
double *y, // umhaengen zwingend notwendig, da u.U. **xx und *newx
// der gleiche aufrufende Pointer ist, d.h. es wird "ueberschrieben"
*z, dummy, t;
int j, k, i, d, n, endfor, w,
spatialdim = nrow - 1,
nxspdim = nx * spatialdim,
ncolM1 = ncol - 1;
z = *newx = (double*) malloc(sizeof(double) * ncol * nx);
if (aniso == NULL) {
assert(ncol == nrow);
for (k=j=0, t=T[XSTART]; j<len; j++, t += T[XSTEP]) {
y = x;
for (i=0; i<nx; i++) {
for (d=0; d<ncolM1; d++, y++) {
z[k++] = *y; //auf keinen Fall nach vorne holen
}
z[k++] = t;
}
}
} else {
for (k=j=0, t=T[XSTART]; j<len; j++, t += T[XSTEP]){
y = x;
for (i=0; i<nxspdim; i+=spatialdim) {
endfor = i + spatialdim;
for (n=d=0; d<ncol; d++) {
dummy = 0.0;
for(w=i; w<endfor; w++) {
dummy += aniso[n++] * y[w];
}
// printf("%d %d %f\n", k, n, t);
z[k++] = dummy + aniso[n++] * t; //auf keinen Fall nach vorne holen
}
}
}
}
}
void x2x(double *x, int nx, double **newx,
double *aniso, int nrow, int ncol) {
double *y = x, // umhaengen zwingend notwendig, da u.U. **xx und *newx
// der gleiche aufrufende Pointer ist, d.h. es wird "ueberschrieben"
*z, dummy;
int k, i, d, n, endfor, w,
nxnrow = nx * nrow;
z = *newx = (double*) malloc(sizeof(double) * ncol * nx);
if (aniso == NULL) {
assert(nrow == ncol);
memcpy(z, y, sizeof(double) * nrow * ncol);
} else {
for (k=i=0; i<nxnrow; i+=nrow) {
endfor = i + nrow;
for (n=d=0; d<ncol; d++) {
dummy = 0.0;
for(w=i; w<endfor; w++) {
dummy += aniso[n++] * y[w];
}
z[k++] = dummy;
}
}
}
}
double *matrixmult(double *m1, double *m2, int dim1, int dim2, int dim3) {
double dummy,
*m0 = (double*) malloc(sizeof(double) * dim1 * dim3);
int i,j,k;
for (i=0; i<dim1; i++) {
for (k=0; k<dim3; k++) {
dummy = 0.0;
for (j=0; j<dim2; j++) {
dummy += m1[i + j * dim1] * m2[j + k * dim2];
}
m0[i + dim1 * k] = dummy;
}
}
return m0;
}
matrix_type Type(double *m, int nrow, int ncol) {
// see also analyse_matrix: can it be unified ???
matrix_type type = TypeIso; // default
double elmt;
int i, endfor, endfor2;
if (m==NULL) {
assert(ncol == nrow);
return TypeIso;
}
if (ncol > nrow) {
endfor = nrow * ncol;
for (i=nrow * ncol; i < endfor; i++)
if (m[i]!= 0.0) return TypeAny;
ncol = nrow;
}
i = nrow * ncol -1;
if (i==0) return TypeIso;
if (ncol==1) {
endfor = nrow - 1;
for (i=1; i<endfor; i++) if (m[i] != 0.0) return TypeAny;
if (m[endfor] == 0.0) return TypeIso;
if (m[0] == 0.0) return TypeIso;
return TypeAny;
}
elmt = m[i];
endfor = i - nrow;
// printf("endfor %d, i=%d\n", endfor, i);
for (i--; i>endfor; i--) { // last row but very last element: not all zero ?
if (m[i] != 0.0) return TypeAny;
}
endfor--;
for (i = 0; i < endfor; ) {
if (m[i++] != elmt) type=TypeDiag; // not iso
endfor2 = i + nrow;
for ( ; i<endfor2; )
if (m[i++] != 0.0) return TypeTimesep; // time compent of x is not
// mixed with spatial components of x; (other way round may be possible
}
if (m[i++] != elmt) type=TypeDiag; // vorletzte Spalte, vorletzte Zeile,
// d.h. vorletztes Diagonalelement
if(m[i++] != 0.0) return TypeTimesep; // vorletzte Spalte, letzte Zeile
// printf("%d %d %d %d\n", i , nrow, (ncol-1), type);
assert(i == nrow * (ncol-1));
return (type);
}
void Transform2NoGridExt(method_type *meth, bool timesep, bool gridexpand,
double **Space, double **Sptime) {
// this function transforms the coordinates according to the anisotropy
// matrix
location_type *loc = meth->loc;
int origdim = loc->timespacedim,
dim = meth->cov->tsdim,
*length = loc->length;
double *aniso=NULL,
*T = loc->T,
*x = loc->x;
aniso = getAnisoMatrix(meth);
if (loc->grid && meth->type <= TypeDiag) {
if (gridexpand) {
if (*Sptime == NULL) {
expandgrid(loc->xgr, length, Sptime, aniso, origdim, dim);
}
if (!loc->Time) *Space = *Sptime;
} else {
grid2grid(loc->xgr, meth->grani, aniso, origdim);
}
} else if (loc->grid) { // aniso not diag
if (!loc->Time) { // grid no time
assert(*Sptime == *Space);
if (*Space == NULL) {
expandgrid(loc->xgr, length, Space, aniso, origdim, dim);
*Sptime = *Space;
}
} else { // grid and time
if (timesep) {
if (*Space == NULL)
expandgrid(loc->xgr, length, Space, aniso, origdim, dim - 1);
} else {
if (*Sptime == NULL)
expandgrid(loc->xgr, length, Sptime, aniso, origdim, dim);
}
}
} else { // nogrid
if (!loc->Time) { // no grid no time
assert(*Sptime == *Space);
if (*Space == NULL) {
x2x(x, length[0], Space, aniso, origdim, dim);
*Sptime = *Space;
}
} else { // no grid but time
if (timesep) {
x2x(x, length[0], Space, aniso, origdim, dim-1);
} else {
xtime2x(x, length[0], T, length[dim-1],Sptime, aniso, origdim, dim);
}
}
}
if (aniso != NULL) free(aniso);
}
void Transform2NoGrid(method_type *meth, bool timesep) {
Transform2NoGridExt(meth, timesep, false, &(meth->space), &(meth->sptime));
}
int InternalGetGridSize(double *x[MAXSIMUDIM], int dim, int *lx)
{
int d;
for (d=0; d<dim; d++) {
if ((x[d][XEND]<=x[d][XSTART]) || (x[d][XSTEP]<=0.0)) {
return ERRORCOORDINATES;
}
lx[d] = 1 + (int)((x[d][XEND]-x[d][XSTART])/x[d][XSTEP]);
}
return NOERROR;
}
double *selectlines(double *m, int *sel, int nsel, int nrow, int ncol) {
// selects lines in matrix m according to sel
int j;
double *red_matrix,
*Red_Matrix = (double *) malloc(sizeof(double) * nsel * ncol),
*endfor = Red_Matrix + nsel * ncol;
for (red_matrix=Red_Matrix; red_matrix<endfor; m+=nrow) {
for (j=0; j<nsel; j++, red_matrix++) {
*red_matrix = m[sel[j]];
}
}
return Red_Matrix;
}
int *selectlines(int *m, int *sel, int nsel, int nrow, int ncol) {
int j;
int *red_matrix,
*Red_Matrix = (int *) malloc(sizeof(int) * nsel * ncol),
*endfor = Red_Matrix + nsel * ncol;
for (red_matrix=Red_Matrix; red_matrix<endfor; m+=nrow) {
for (j=0; j<nsel; j++, red_matrix++) {
*red_matrix = m[sel[j]];
}
}
return Red_Matrix;
}
void cum_dollar(method_type *meth, int origdim, cov_model *cov,
method_type *ziel) {
/*
cpy meth to ziel with dollar...
schiebe scale/proj/aniso in meth->cscale/cproj/caniso
Sowie: setze meth->xdimout
*/
int i, total,
nproj = cov->nrow[DPROJ],
*proj = (int *) cov->p[DPROJ];
assert(meth != ziel);
if (nproj > 0) { // DSCALE is given at the same time !
if (meth->caniso == NULL) {
ziel->cscale = meth->cscale * cov->p[DSCALE][0];
if (meth->cproj == NULL) {
ziel->cproj = (int*) malloc(sizeof(int) * nproj);
memcpy(ziel->cproj, proj, nproj * sizeof(int));
} else {
ziel->cproj = selectlines(meth->cproj, proj, nproj, -9999, 1);
}
ziel->type = TypeDiag;
} else { //cov->p[DPROJ] != NULL, mcaniso != NULL
ziel->caniso = selectlines(meth->caniso, proj,
nproj, origdim, meth->xdimout);
}
ziel->xdimout = nproj;
} else {
if (cov->p[DANISO] == NULL) {
total = origdim * meth->xdimout;
ziel->xdimout = meth->xdimout;
if (meth->caniso != NULL) {
double invscale = 1.0 / cov->p[DSCALE][0];
ziel->caniso = (double *) malloc(total * sizeof(double));
for (i=0; i<total; i++) ziel->caniso[i] = meth->caniso[i] * invscale;
} else {
ziel->cscale = meth->cscale * cov->p[DSCALE][0];
}
} else { // cov->p[DANISO] != NULL
int ncol = cov->ncol[DANISO];
total = origdim * cov->ncol[DANISO];
if (meth->caniso == NULL) {
if (meth->cproj == NULL) {
ziel->caniso = (double *) malloc(total * sizeof(double));
// printf("%d %d %d\n", total, ncol, cov->nrow[DANISO]);
// assert(false);
memcpy(ziel->caniso, cov->p[DANISO], total * sizeof(double));
} else {
double *ca, *endfor,
*aniso = cov->p[DANISO];
ziel->caniso =
(double *) calloc(origdim * ncol, sizeof(double));
// timespacedim x ncol
endfor = ziel->caniso + ncol * origdim;
for (ca = ziel->caniso; ca<endfor; ca += origdim) {
for (i=0; i<nproj; i++, aniso++) {
ca[meth->cproj[i]] = *aniso;
}
}
free(meth->cproj);
meth->cproj = NULL;
}
} else {
assert(meth->xdimout == cov->nrow[DANISO]);
ziel->caniso = matrixmult(meth->caniso, cov->p[DANISO], origdim,
meth->xdimout, ncol);
// caniso = caniso %*%cov->p
}
if (origdim != ncol) ziel->type = TypeAny;
else ziel->type = Type(ziel->caniso, origdim, ncol);
ziel->xdimout = ncol;
}
}
}
