https://github.com/cran/RandomFields
Tip revision: e5a7a2f272b7834f96c925ced7acfa0c6456a87f authored by Martin Schlather on 17 April 2017, 22:09:51 UTC
version 3.1.50
version 3.1.50
Tip revision: e5a7a2f
getNset.cc
/*
Authors
Martin Schlather, schlather@math.uni-mannheim.de
(library for simulation of random fields)
Copyright (C) 2001 -- 2015 Martin Schlather,
This program is free software; you can redist ribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 3
of the License, or (at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/
#include <R.h>
#include <Rdefines.h>
#include <Rmath.h>
#include <stdio.h>
//#include <stdlib.h>
#include <string.h>
#include <R_ext/Linpack.h>
#include "RF.h"
#include "primitive.h"
#include "Coordinate_systems.h"
void LOC_SINGLE_NULL(location_type *loc, int len) {
int d;
loc->spatialdim = loc->timespacedim = loc->lx = loc->ly = loc->xdimOZ = -1;
for (d=0; d<MAXSIMUDIM; d++) {
loc->xgr[d] = loc->ygr[d] = NULL;
}
loc->totalpoints = loc->spatialtotalpoints = 0;
loc->grid = loc->distances = loc->Time = false;
loc->delete_x = loc->delete_y = true;
loc->x = loc->y = loc->caniso = NULL;
loc->T[0] = loc->T[1] = loc->T[2] = 0.0;
loc->i_row = loc->i_col = I_COL_NA;
loc->cani_ncol = loc->cani_nrow = NA_INTEGER;
loc->len = len;
}
void LOC_NULL(location_type **Loc, int len) {
int i;
for (i=0; i<len; i++) LOC_SINGLE_NULL(Loc[i], len);
}
location_type **LOCLIST_CREATE(int n) {
int i;
location_type **loc = (location_type**) CALLOC(n, sizeof(location_type*));
//printf("sizeof = %d %d %d\n", sizeof(location_type), sizeof(coord_type), n); BUG;
for (i=0; i<n; i++) loc[i] = (location_type*) MALLOC(sizeof(location_type));
LOC_NULL(loc, n);
assert(loc[0]->len >= 1 && loc[0]->len <= 1e6);
return loc;
}
void LOCLIST_CREATE(cov_model *cov, int n) {
cov->ownloc = LOCLIST_CREATE(n);
}
void LOC_SINGLE_DELETE(location_type **Loc) {
location_type *loc = *Loc;
if (loc != NULL) {
if (loc->x != NULL) {
//printf("single %d\n", loc->delete_x);
if (loc->delete_y) FREE(loc->y);
if (loc->delete_x) {
//printf("single del x %ld\n", loc->x);
assert(loc->x != NULL);
UNCONDFREE(loc->x);
}
//printf("don LOC_SINGLE_DEL\n");
}
FREE(loc->caniso);
// it may happen that both are set ! Especially after calling
// partial_loc_set in variogramAndCo.cc
if (loc->spatialdim>0) {
if (loc->delete_y) FREE(loc->ygr[0]); ///
if (loc->delete_x) FREE(loc->xgr[0]); ///
}
UNCONDFREE(*Loc);
}
}
void LOC_DELETE(location_type ***Loc) {
if (*Loc == NULL) return;
int i,
len = (*Loc)[0]->len;
//printf("len=%d %ld %ld\n", len, (*Loc) + 0 , (*Loc) + 1); assert(len <= 2);
for (i=0; i<len; i++) {
//printf("loc_del i=%d\n", i);
LOC_SINGLE_DELETE( (*Loc) + i );
}
UNCONDFREE(**Loc);
UNCONDFREE(*Loc);
}
listoftype * LIST_CREATE(int len, int type) {
if (len <= 0) BUG;
listoftype* q = (listoftype *) MALLOC(sizeof(listoftype));
q->lpx = (double **) CALLOC(len, sizeof(double *));
q->ncol = (int*) CALLOC(len, sizeof(int));
q->nrow = (int*) CALLOC(len, sizeof(int));
q->deletelist = true;
q->len = len;
q->type = type;
return q;
}
void LIST_DELETE(listoftype **x) {
if (x == NULL) return;
listoftype *q = *x;
if (q != NULL) {
// printf("list-delete %ld %d\n", (long int) q, q->len);
assert(q->lpx != NULL);
if (q->deletelist) {
for (int i=0; i<q->len; i++) {
//printf("list_del %d q-len=%d\n", i, q->len);
FREE(q->lpx[i]);
}
FREE(q->lpx);
FREE(q->ncol);
FREE(q->nrow);
}
FREE(*x);
}
}
void listpt(listoftype **To, listoftype *p, int len, SEXPTYPE type,
bool force_allocating) {
if (*To == NULL || force_allocating) {
*To = (listoftype *) MALLOC(sizeof(listoftype));
}
listoftype *q = *To;
q->lpx = p->lpx;
q->ncol = p->ncol;
q->nrow = p->nrow;
q->deletelist = false;
q->len = len;
q->type = type;
}
void listcpy(listoftype **To, listoftype *p, bool force_allocating) {
// force_allocating in case of "c ovcpy"
int size,
len = p->len,
sizeint = len * sizeof(int);
if (p->type == LISTOF + REALSXP) {
size = sizeof(double);
} else BUG;
if (*To == NULL || force_allocating) *To = LIST_CREATE(len, p->type);
listoftype *q = *To;
for (int j=0; j<len; j++) {
int n = size * p->nrow[j] * p->ncol[j];
//printf("list cpy j=%d %d\n", j, n);
if (q->lpx[j] == NULL) q->lpx[j] = (double*) MALLOC(n);
MEMCOPY(q->lpx[j], p->lpx[j], n);
}
MEMCOPY(q->nrow, p->nrow, sizeint);
MEMCOPY(q->ncol, p->ncol, sizeint);
}
//void GLOBAL_NULL() {
// ReSetGlobal();
//}
//
//void SetTemporarilyGlobal(globalparam *gp) {
// if (!GLOBALORIG.set) {
// MEMCOPY(&(GLOBALORIG.param), &GLOBAL, sizeof(globalparam));
// GLOBALORIG.set = true;
// }
// MEMCOPY(&GLOBAL, gp, sizeof(globalparam));
//}
//
//void ReSetGlobal() {
// if (GLOBALORIG.set) {
// MEMCOPY(&GLOBAL, &(GLOBALORIG.param), sizeof(globalparam));
// GLOBALORIG.set = false;
// }
//}
// 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 MPPPROPERTIES_NULL(mpp_properties *Mpp) {
// Mpp->refradius=
int i;
for (i=0; i<MAXMPPVDIM; i++) Mpp->maxheights[i] = RF_INF;
Mpp->unnormedmass = RF_NA;
Mpp->mM = Mpp->mMplus = NULL;
// Mpp->refsd =
//Mpp->totalmass = RF_NA;
//Mpp->methnr = -1;
//Mpp->loc_done = false;
}
void MPPPROPERTIES_DELETE(mpp_properties *Mpp) {
FREE(Mpp->mM);
FREE(Mpp->mMplus);
}
void COV_DELETE_WITHOUTSUB(cov_model **Cov) {
cov_model *cov = *Cov;
//printf("deleting %s\n", NAME(cov));
assert(cov != NULL);
int i, j,
last = (cov->nr < 0) ? MAXPARAM : CovList[cov->nr].kappas;
for (i=0; i<last; i++) {
int type = CovList[cov->nr].kappatype[i];
if (!PisNULL(i)) {
if (isRObject(type)) {
sexp_type *S = PSEXP(i);
if (S->Delete) R_ReleaseObject(S->sexp);
} else if (type >= LISTOF) {
// printf("cov_delete_without %s %d; %d %d %d\n", NAME(*Cov), cov->zaehler, NROW(i), NCOL(i), i);
// printf("cov_delete_without %s %d %d\n", NAME(*Cov), ((listoftype *) (cov->px[i]))->deletelist, i);
//if (PLIST(i)->deletelist) {
LIST_DELETE((listoftype **) (cov->px + i));
//}
}
PFREE(i);
cov->ncol[i] = cov->nrow[i] = SIZE_NOT_DETERMINED; // ==0
}
}
MPPPROPERTIES_DELETE(&(cov->mpp));
if (cov->ownkappanames != NULL) {
int kappas = CovList[cov->nr].kappas;
for (j=0; j<kappas; j++) FREE(cov->ownkappanames[j]);
UNCONDFREE(cov->ownkappanames);
}
QFREE;
// important check in combination with above; can be easily removed or
// generalised !!!!
FREE(cov->MLE);
//MPPPROPERTIES_DELETE(&(cov->mpp));
cov->prevloc = NULL;
LOC_DELETE(&(cov->ownloc));
if (cov->key != NULL) {
// printf("deleting key %s\n", CovList[cov->key->nr].name);
COV_DELETE(&(cov->key));
}
if (cov->rf != NULL && cov->origrf) UNCONDFREE(cov->rf);
ce_DELETE(&(cov->Sce));
localCE_DELETE(&(cov->SlocalCE));
approxCE_DELETE(&(cov->SapproxCE));
direct_DELETE(&(cov->Sdirect));
hyper_DELETE(&(cov->Shyper));
mixed_DELETE(&(cov->Smixed));
nugget_DELETE(&(cov->Snugget));
plus_DELETE(&(cov->Splus));
sequ_DELETE(&(cov->Ssequ));
// SPECTRAL_DELETE(&(cov->Sspectral));
trend_DELETE(&(cov->Strend));
tbm_DELETE(&(cov->Stbm));
br_DELETE(&(cov->Sbr));
get_DELETE(&(cov->Sget));
pgs_DELETE(&(cov->Spgs));
set_DELETE(&(cov->Sset));
polygon_DELETE(&(cov->Spolygon));
rect_DELETE(&(cov->Srect));
dollar_DELETE(&(cov->Sdollar));
gatter_DELETE(&(cov->Sgatter));
earth_DELETE(&(cov->Searth));
extra_DELETE(&(cov->Sextra));
RU_solve_DELETE(&(cov->Ssolve));
biwm_DELETE(&(cov->Sbiwm));
inv_DELETE(&(cov->Sinv));
scatter_DELETE(&(cov->Sscatter));
mcmc_DELETE(&(cov->Smcmc));
// SELECT_DELETE(&(cov->Sselect));
gen_DELETE(&(cov->Sgen));
likelihood_DELETE(&(cov->Slikelihood));
covariate_DELETE(&(cov->Scovariate));
simu_type *simu = &(cov->simu);
simu->active = simu->pair = false;
simu->expected_number_simu = 0;
UNCONDFREE(*Cov);
}
void COV_DELETE_WITHOUT_LOC(cov_model **Cov) {
cov_model *cov = *Cov;
int i,
nsub = CovList[cov->nr].maxsub;
if (cov == NULL) {
warning("*cov is NULL in COV_DELETE");
return;
}
//printf("COV_DEL %s\n", NAME(*Cov));
// PMI(cov, "delete");
assert(cov != NULL);
for (i=0; i<MAXPARAM; i++) { // cov->sub[i] // seit 1.10.07: luecken erlaubt
// bei PSgen !
if (cov->kappasub[i] != NULL) {
// printf("%d %d %ld\n",i, MAXPARAM, cov->kappasub[i]);
//PMI(cov);
//assert(false);
COV_DELETE_WITHOUT_LOC(cov->kappasub + i);
}
}
for (i=0; i<nsub; i++) { // cov->sub[i] // seit 1.10.07: luecken erlaubt
// bei PSgen !
// Achtung nach nsub koennen "blinde" Untermodelle kommen, die
// nicht geloescht werden duerfen!!
if (cov->sub[i] != NULL) COV_DELETE_WITHOUT_LOC(cov->sub + i);
}
COV_DELETE_WITHOUTSUB(Cov);
}
void COV_DELETE_(cov_model **Cov) {
cov_model *cov = *Cov;
// if (*Cov == NULL) crash(*Cov);
assert(*Cov != NULL);
// BUG;PMI(cov);
// printf("del %s\n", NAME(cov));
if (cov->calling == NULL) LOC_DELETE(&(cov->prevloc));
COV_DELETE_WITHOUT_LOC(Cov);
*Cov = NULL;
}
void SIMU_NULL(simu_type *simu) {
simu->pair = simu->active = false;
simu->expected_number_simu =0;
}
static int CovZaehler = 0;
void COV_ALWAYS_NULL(cov_model *cov) {
cov->zaehler = CovZaehler++;
cov->calling = NULL;
cov->prevloc = cov->ownloc = NULL;
cov->checked = false;
cov->MLE = NULL;
cov->key = NULL;
cov->rf = NULL;
cov->mpp.mM = cov->mpp.mMplus = NULL;
cov->mpp.moments = MISMATCH;
cov->Sce = NULL;
cov->SlocalCE = NULL;
cov->SapproxCE = NULL;
cov->Sdirect = NULL;
cov->Shyper = NULL;
cov->Smixed = NULL;
cov->Snugget = NULL;
cov->Splus = NULL;
cov->Ssequ = NULL;
cov->Stbm = NULL;
cov->Strend = NULL;
cov->Sbr = NULL;
cov->Sget = NULL;
cov->Spgs = NULL;
cov->Sset = NULL;
cov->Spolygon = NULL;
cov->Srect = NULL;
cov->Sdollar = NULL;
cov->Sgatter = NULL;
cov->Searth = NULL;
cov->Sextra = NULL;
cov->Ssolve = NULL;
cov->Sbiwm = NULL;
cov->Sinv = NULL;
cov->Sscatter = NULL;
cov->Smcmc = NULL;
//cov->Sselect = NULL;
cov->Sgen = NULL;
cov->Slikelihood = NULL;
cov->Scovariate = NULL;
cov->Sbistable = NULL;
cov->fieldreturn = cov->origrf = false;
cov->initialised = false;
}
void COV_NULL(cov_model *cov) {
int i, j;
COV_ALWAYS_NULL(cov);
cov->nr = cov->gatternr = cov->secondarygatternr = MISMATCH;
for (i=0; i<MAXPARAM; i++) {
PtoNULL(i);
cov->ncol[i] = cov->nrow[i] = 0;
cov->kappasub[i] = NULL;
}
for (i=0; i<MAXTAYLOR; i++) {
for (j=(int) TaylorConst; j<=(int) TaylorPow; j++)
cov->taylor[i][j] = cov->tail[i][j] = 0.0;
for ( ; j<=(int) TaylorExpPow; j++)
cov->tail[i][j] = 0.0;
}
cov->q = NULL;
cov->qlen = 0;
for (i=0; i<MAXSUB; i++) cov->sub[i] = NULL;
cov->nsub = 0;
cov->user_given = ug_internal; // to know which models are given by the
// user and which ones have been added internally
cov->role = ROLE_UNDEFINED;
cov->typus = UndefinedType;
cov->method = Forbidden;
cov->tsdim = cov->xdimprev = cov->xdimown =cov->maxdim = cov->xdimgatter =
UNSET;
cov->vdim[0] = cov->vdim[1] = MISMATCH;
cov->ownkappanames = NULL;
cov->domown = cov->domprev = DOMAIN_MISMATCH;
cov->isoown = cov->isoprev = ISO_MISMATCH;
cov->logspeed = RF_NA;
cov->delflag = 0;
cov->full_derivs = cov->rese_derivs = MISMATCH;
cov->ptwise_definite = pt_undefined;
cov->deterministic = true;
cov->monotone = MISMATCH;
// cov->total_n = MISMATCH;
// cov->total_sum = 0.0;
cov->finiterange = cov->loggiven =
// cov->diag = cov->semiseparatelast = cov->separatelast =
cov->matrix_indep_of_x = false;
cov->hess = cov->tbm2num = NOT_IMPLEMENTED;
for (i=0; i<=Nothing; i++) cov->pref[i] = PREF_BEST;
// mpp und extrG muessen auf NONE sein !
for (; i<=Forbidden; i++) cov->pref[i] = PREF_NONE;
cov->taylorN = cov->tailN = 0;
MPPPROPERTIES_NULL(&(cov->mpp));
SIMU_NULL(&(cov->simu));
// cov->localcov=NULL;
}
void FFT_NULL(FFT_storage *FFT)
{
if (FFT == NULL) return;
FFT->work = NULL;
FFT->iwork = NULL;
}
void FFT_destruct(FFT_storage *FFT)
{
FREE(FFT->iwork);
FREE(FFT->work);
FFT_NULL(FFT);
}
void ce_DELETE(ce_storage **S) {
ce_storage *x = *S;
if (x != NULL) {
int l,
vdim = x->vdim,
vdimSQ = vdim * vdim;
if (x->c!=NULL) {
for(l=0; l<vdimSQ; l++) FREE(x->c[l]);
UNCONDFREE(x->c);
}
if (x->d!=NULL) {
for(l=0; l<vdim; l++) FREE(x->d[l]);
UNCONDFREE(x->d);
}
FFT_destruct(&(x->FFT));
FREE(x->aniso);
FREE(x->gauss1);
FREE(x->gauss2);
UNCONDFREE(*S);
}
}
void ce_NULL(ce_storage* x){
if (x == NULL) return;
FFT_NULL(&(x->FFT));
x->positivedefinite = FALSE;
x->trials = -1;
x->c = x->d = NULL;
x->smallestRe = x->largestAbsIm = RF_NA;
x->aniso = NULL;
x->stop = false;
// cur_call_odd = false;
// new_simulation_next = true;
x->gauss1 = x->gauss2 = NULL;
// int i;
//for (i=0; i<MAXCEDIM; i++) x->[i] = NULL;
}
void localCE_DELETE(localCE_storage**S)
{
localCE_storage* x = *S;
if (x!=NULL) {
FREE(x->correction);
UNCONDFREE(*S);
}
}
void localCE_NULL(localCE_storage* x){
// int i; for (i=0; i<MAXCEDIM; i++)
if (x == NULL) return;
x->correction = NULL;
}
void approxCE_DELETE(approxCE_storage **S) {
approxCE_storage* x = * S;
if (x != NULL) {
FREE(x->idx);
UNCONDFREE(*S);
}
}
void approxCE_NULL(approxCE_storage* x){
if (x == NULL) return;
x->idx = NULL;
}
void direct_DELETE(direct_storage ** S) {
direct_storage *x = *S;
if (x!=NULL) {
FREE(x->G);
UNCONDFREE(*S);
}
}
void direct_NULL(direct_storage *x) {
if (x == NULL) return;
x->G = NULL;
}
void hyper_DELETE(hyper_storage **S) {
hyper_storage *x = *S;
if (x != NULL) {
UNCONDFREE(*S);
}
}
void hyper_NULL(hyper_storage* x) {
if (x == NULL) return;
}
void mixed_DELETE(mixed_storage ** S) {
mixed_storage *x = *S;
if (x!=NULL) {
FREE(x->Xb);
FREE(x->mixedcov);
UNCONDFREE(*S);
}
}
void mixed_NULL(mixed_storage* x) {
x->Xb = NULL;
x->mixedcov = NULL;
x->initialized = false;
}
void nugget_NULL(nugget_storage *x){
if (x == NULL) return;
x->pos = NULL;
x->red_field = NULL;
}
void nugget_DELETE(nugget_storage ** S)
{
nugget_storage *x = *S;
if (x != NULL) {
FREE(x->pos);
FREE(x->red_field);
UNCONDFREE(*S);
}
}
void plus_NULL(plus_storage *x){
if (x != NULL) {
int i;
for (i=0; i<MAXSUB; i++) x->keys[i] = NULL;
}
}
void plus_DELETE(plus_storage ** S){
plus_storage *x = *S;
if (x != NULL) {
int i;
for (i=0; i<MAXSUB; i++)
if (x->keys[i] != NULL) COV_DELETE(x->keys + i);
UNCONDFREE(*S);
}
}
void sequ_NULL(sequ_storage *x){
if (x == NULL) return;
x->U11 = NULL;
x->U22 = NULL;
x->MuT = NULL;
x->G = NULL;
x->Cov21 = NULL;
x->Inv22 = NULL;
x->res0 = NULL;
}
void sequ_DELETE(sequ_storage ** S){
sequ_storage *x = *S;
if (x!=NULL) {
FREE(x->U11);
FREE(x->U22);
FREE(x->MuT);
FREE(x->G);
FREE(x->Cov21);
FREE(x->Inv22);
FREE(x->res0);
UNCONDFREE(*S);
}
}
void spectral_DELETE(spectral_storage **S)
{
spectral_storage *x = *S;
if (x!=NULL) {
// do NOT delete cov --- only pointer
// spectral_storage *x; x = *((spectral_storage**)S);
UNCONDFREE(*S);
}
}
void spectral_NULL(spectral_storage *x)
{
if (x!=NULL) {
}
}
void trend_DELETE(trend_storage ** S) {
trend_storage *x = *S;
if (x!=NULL) {
FREE(x->x);
FREE(x->xi);
FREE(x->evalplane);
FREE(x->powmatrix);
UNCONDFREE(*S);
}
}
void trend_NULL(trend_storage* x) {
x->x = NULL;
x->xi = NULL;
x->evalplane = NULL;
x->powmatrix = NULL;
}
void tbm_DELETE(tbm_storage **S)
{
tbm_storage *x = *S;
if (x!=NULL) {
UNCONDFREE(*S);
}
}
void tbm_NULL(tbm_storage* x) {
if (x == NULL) return;
x->ce_dim = 0;
x->err = ERRORFAILED;
}
void BRTREND_DELETE(double **BRtrend, int trendlen) {
int j;
if (BRtrend == NULL) return;
for (j=0; j<trendlen; j++) FREE(BRtrend[j]);
}
void br_DELETE(br_storage **S) {
br_storage *brS = *S;
if (brS != NULL) {
int i;
if (brS->trend != NULL) {
BRTREND_DELETE(brS->trend, brS->trendlen);
UNCONDFREE(brS->trend);
}
FREE(brS->shiftedloc);
FREE(brS->loc2mem);
if (brS->countvector != NULL) {
for (i=0; i<brS->vertnumber; i++) FREE(brS->countvector[i]);
UNCONDFREE(brS->countvector);
}
FREE(brS->areamatrix);
FREE(brS->logvertnumber);
FREE(brS->locindex);
FREE(brS->suppmin);
FREE(brS->suppmax);
FREE(brS->locmin);
FREE(brS->locmax);
FREE(brS->loccentre);
FREE(brS->mem2loc);
FREE(brS->newx);
if (brS->vario != NULL) COV_DELETE(&(brS->vario));
FREE(brS->lowerbounds);
if (brS->submodel != NULL) COV_DELETE(&(brS->submodel));
UNCONDFREE(*S);
}
}
void br_NULL(br_storage *brS) {
brS->trend = NULL;
brS->trendlen = 0;
brS->next_am_check = NA_INTEGER;
brS->shiftedloc = NULL;
brS->mem2loc = brS->locindex = brS->loc2mem = NULL;
brS->memcounter = 0;
brS->newx = brS->suppmin = brS->suppmax =
brS->locmin = brS->locmax = brS->loccentre = NULL;
int i;
for (i=0; i<MAXSUB; i++) {
brS->sub[i] = NULL;
}
brS->lowerbounds = NULL;
brS->vario = brS->submodel = NULL;
brS->countvector = NULL;
brS->areamatrix = NULL;
brS->logvertnumber = NULL;
}
void pgs_DELETE(pgs_storage **S)
{
pgs_storage *x = *S;
if (x!=NULL) {
// Huetchen
FREE(x->v);
FREE(x->y);
FREE(x->xgr[0]);
FREE(x->pos);
FREE(x->min);
FREE(x->max);
FREE(x->single);
FREE(x->total);
FREE(x->halfstepvector);
FREE(x->localmin);
FREE(x->localmax);
FREE(x->minmean);
FREE(x->maxmean);
// rf_interface.cc
FREE(x->supportmin);
FREE(x->supportmax);
FREE(x->supportcentre);
FREE(x->own_grid_start);
FREE(x->own_grid_step);
FREE(x->own_grid_len);
FREE(x->gridlen);
FREE(x->end);
FREE(x->start);
FREE(x->delta);
FREE(x->nx);
FREE(x->xstart);
FREE(x->x);
FREE(x->inc);
// variogramAndCo.cc:
FREE(x->endy);
FREE(x->startny);
FREE(x->ptrcol);
FREE(x->ptrrow);
FREE(x->C0x);
FREE(x->C0y);
FREE(x->cross);
FREE(x->z);
FREE(x->Val);
if (x->cov != NULL) {
cov_model *dummy = x->cov;
if (x->cov->Spgs != NULL && x->cov->Spgs->cov != NULL &&
x->cov->Spgs->cov->Spgs == x) {
x->cov->Spgs->cov = NULL;// um "Unendlich"-Schleife zu vermeiden!
} else {
assert(x->cov->Spgs == NULL || x->cov->Spgs->cov == NULL ||
x->cov->Spgs->cov->Spgs == x);
}
x->cov = NULL; // Sicherheitshalber
COV_DELETE(&(dummy));
}
UNCONDFREE(*S);
}
}
void pgs_NULL(pgs_storage* x) {
if (x == NULL) return;
x->log_density = x->intensity = x->globalmin = x->currentthreshold =
RF_NEGINF;
x->flat = x->estimated_zhou_c = x->logmean = false;
x->totalmass = RF_NA;
x->alpha = 1.0; // !! standard
x->size = x->rowscols = -1;
x->zhou_c = RF_NA;
x->sq_zhou_c = x->sum_zhou_c = 0;
x->n_zhou_c = 0;
// Huetchen
x->v = x->y = x->xgr[0] = NULL;
x->pos = x->min = x->max = NULL;
x->single = x->total = x->halfstepvector = NULL;
x->localmin = x->localmax =x->minmean = x->maxmean = NULL;
// rf_interface.cc
x->supportmin = x->supportmax = x->supportcentre = //
x->own_grid_start = x->own_grid_step = x->own_grid_len = NULL; //
x->gridlen = x->end = x->start = x->delta = x->nx = NULL; //
x->xstart = x->x = x->inc = NULL;
// variogramAndCo.cc
x->endy = x->startny = x->ptrcol = x->ptrrow = NULL;
x->C0x = x->C0y = x->cross = x->z = NULL;
x->Val = NULL;
//x->param_set = NULL;
x->cov = NULL;
}
void set_DELETE(set_storage **S) {
set_storage *x = *S;
if (x!=NULL) {
UNCONDFREE(*S);
}
}
void set_NULL(set_storage* x) {
if (x == NULL) return;
x->remote = NULL;
x->set = NULL;
x->variant = 0;
}
void polygon_DELETE(polygon_storage **S)
{
polygon_storage *x = *S;
if (x != NULL) {
if (x->vdual != NULL) {
int i;
for (i=0; i<x->n_vdual; i++) FREE(x->vdual[i]);
UNCONDFREE(x->vdual);
}
FREE(x->vprim);
if (x->P != NULL) {
freePolygon(x->P);
UNCONDFREE(x->P);
}
}
UNCONDFREE(*S);
}
void polygon_NULL(polygon_storage* x) {
if (x == NULL) return;
x->vprim = NULL;
x->vdual = NULL;
x->n_vdual = x->n_vertex = x->n_v = 0;
polygon *P = x->P;
if (P == NULL) BUG;
P->e = NULL;
P->v = NULL;
P->n = 0;
}
void rect_DELETE(rect_storage **S){
rect_storage *x = *S;
if (x!=NULL) {
FREE(x->value);
FREE(x->weight);
FREE(x->tmp_weight);
FREE(x->right_endpoint);
FREE(x->ysort);
FREE(x->z);
FREE(x->squeezed_dim);
FREE(x->asSign);
FREE(x->i);
UNCONDFREE(*S);
}
}
void rect_NULL(rect_storage* x) {
if (x == NULL) return;
x->tmp_n = x->nstep = -999;
x->value = x->weight = x->tmp_weight =
x->right_endpoint = x->ysort = x->z = NULL;
// for PMI
// x->inner = x->inner_const = x->inner_pow = x->outer =
// x->outer_const = x->outer_pow = x->outer_pow_const = x->step = RF_NA;
x->squeezed_dim = x->asSign = x->i = NULL;
}
void dollar_DELETE(dollar_storage **S)
{
dollar_storage *x = *S;
if (x!=NULL) {
FREE(x->z);
FREE(x->z2);
FREE(x->y);
FREE(x->y2);
FREE(x->save_aniso);
FREE(x->inv_aniso);
FREE(x->var);
FREE(x->nx);
FREE(x->len);
FREE(x->total);
FREE(x->cumsum);
UNCONDFREE(*S);
}
}
void dollar_NULL(dollar_storage* x) {
if (x == NULL) return;
x->z = x->z2 = x->y = x->y2 = x->save_aniso = x->inv_aniso = x->var = NULL;
x->cumsum = x->nx = x->total = x->len = NULL;
x->proj = 0;
}
void gatter_DELETE(gatter_storage **S)
{
gatter_storage *x = *S;
if (x!=NULL) {
FREE(x->z);
UNCONDFREE(*S);
}
}
void gatter_NULL(gatter_storage* x) {
if (x == NULL) return;
x->z = NULL;
// x->zsys = NULL;
}
void earth_DELETE(earth_storage **S)
{
earth_storage *x = *S;
if (x!=NULL) {
FREE(x->X);
FREE(x->Y);
FREE(x->U);
FREE(x->V);
UNCONDFREE(*S);
}
}
void earth_NULL(earth_storage* x) {
if (x == NULL) return;
x->X = x->Y = x->U = x->V = NULL;
}
void extra_DELETE(extra_storage **S)
{
extra_storage *x = *S;
if (x!=NULL) {
FREE(x->a);
FREE(x->b);
FREE(x->c);
LOC_DELETE(&(x->loc));
UNCONDFREE(*S);
}
}
void extra_NULL(extra_storage* x) {
if (x == NULL) return;
x->a = x->b = x->c = NULL;
x->loc =NULL;
}
void biwm_DELETE(biwm_storage **S)
{
biwm_storage *x = *S;
if (x!=NULL) {
UNCONDFREE(*S);
}
}
void biwm_NULL(biwm_storage* x) {
if (x == NULL) return;
}
void bistable_DELETE(bistable_storage **S)
{
bistable_storage *x = *S;
if (x!=NULL) {
UNCONDFREE(*S);
}
}
void bistable_NULL(bistable_storage* x) {
if (x == NULL) return;
}
void inv_DELETE(inv_storage **S) {
inv_storage *x = *S;
if (x!=NULL) {
FREE(x->v);
FREE(x->wert);
UNCONDFREE(*S);
}
}
void inv_NULL(inv_storage* x) {
if (x == NULL) return;
x->v = NULL;
x->wert = NULL;
}
void scatter_DELETE(scatter_storage **S) {
scatter_storage *x = *S;
if (x!=NULL) {
FREE(x->nx);
FREE(x->min);
FREE(x->max);
FREE(x->step);
FREE(x->x);
FREE(x->xmin);
FREE(x->value);
UNCONDFREE(*S);
}
}
void scatter_NULL(scatter_storage* x) {
if (x == NULL) return;
x->dim = x->vdim = -1;
x->min = x->max = x->nx = NULL;
x->step = x->value = x->xmin = x->x = NULL;
}
void mcmc_DELETE(mcmc_storage **S) {
mcmc_storage *x = *S;
if (x!=NULL) {
FREE(x->pos);
FREE(x->deltapos);
FREE(x->propdelta);
FREE(x->proposed);
UNCONDFREE(*S);
}
}
void mcmc_NULL(mcmc_storage* x) {
if (x == NULL) return;
x->pos = x->deltapos = x->propdelta = x->proposed = NULL;
//x->done = false;
x->integral = 0.0;
}
void get_NULL(get_storage *s){
s->orig = s->get_cov = NULL;
s->idx = NULL;
s->param_nr = -1;
}
void get_DELETE(get_storage **S){
get_storage *x = *S;
if (x != NULL) {
FREE(x->idx);
UNCONDFREE(*S);
}
}
void gen_NULL(gen_storage *x) {
int d;
if (x == NULL) return;
// x->mpp.newx = NULL;
// for (d=0; d<MAXMPPDIM; d++)
// x->window.min[d] = x->window.max[d] = x->window.centre[d] = RF_NA;
x->check = x->prodproc_random = x->dosimulate = true;
x->Sspectral.phistep2d = x->Sspectral.phi2d = x->Sspectral.prop_factor
= RF_NA;
x->Sspectral.grid = false;
x->spec.nmetro = -1;
x->spec.sigma = -1.0;
x->spec.density = NULL;
for (d=0; d<MAXTBMSPDIM; d++) {
x->spec.E[d] = x->spec.sub_sd_cum[d] = RF_NA;
}
}
void gen_DELETE(gen_storage **S) {
gen_storage *x = *S;
if (x != NULL) {
UNCONDFREE(*S);
}
}
void likelihood_info_DELETE(likelihood_info *x) {
// nicht: Var, ptvariance !
FREE(x->Matrix);
}
void likelihood_info_NULL(likelihood_info *x) {
if (x == NULL) return;
x->varmodel = model_undefined;
x->Var = NULL; // DO NOT FREE
x->Matrix = x->pt_variance = NULL; // DO NOT FREE
x->globalvariance = x->trans_inv = x->isotropic = false;
x->newxdim = x->neffect = 0;
for (int i=0; i<MAX_LIN_COMP; i++) {
x->effect[i] = 0;
x->nas[i] = 0;
}
}
void likelihood_NULL(likelihood_storage *x) {
if (x == NULL) return;
x->datasets = NULL;
x->X = x->YhatWithoutNA = NULL;
x->XtX = x->XCY = x->XitXi = x->C = x->CinvXY = x->matrix = x->betavec =
x->sumY = x->work = x->Cwork = x->Xwork = NULL;
x->sets = x->fixedtrends = x->dettrends = x->random = x->max_total_data =
x->maxbeta = 0;
for (int i=0; i<MAX_LIN_COMP; i++) {
x->betas[i] = x->nas[i] = x->nas_det[i] = x->nas_fixed[i] =
x->nas_random[i] = x->nas_boxcox = x->nas_boxcox_mu = 0;
x->cov_fixed[i] = x->cov_det[i] = x->cov_random[i] = NULL; // DO NOT FREE
}
x->where_fixed = NULL;
x->data_nas = NULL;
x->dettrend_has_nas = x->fixedtrend_has_nas = x->random_has_nas =
x->betas_separate = x->ignore_trend = x->data_has_nas = false;
for (int i=0; i<MAX_LIN_COMP;i++) {
x->betanames[i]=NULL;
}
likelihood_info_NULL(&(x->info));
}
void likelihood_DELETE(likelihood_storage **S) {
likelihood_storage *x = *S;
if (x != NULL) {
//printf("likelihood del\n");
LIST_DELETE(&(x->datasets));
if (x->X != NULL) for (int i=0; i<x->sets; i++) FREE(x->X[i]);
FREE(x->X);
if (x->YhatWithoutNA != NULL) for (int i=0; i<x->sets; i++)
FREE(x->YhatWithoutNA [i]);
FREE(x->YhatWithoutNA);
FREE(x->XCY);
FREE(x->XtX);
FREE(x->XitXi);
FREE(x->C);
FREE(x->CinvXY);
FREE(x->Cwork);
FREE(x->Xwork);
FREE(x->matrix);
FREE(x->betavec);
FREE(x->work);
//FREE(x->Yhat);
FREE(x->where_fixed);
FREE(x->sumY);
FREE(x->data_nas);
int end = x->betas[x->fixedtrends];
for (int i=0; i<end; i++) FREE(x->betanames[i]);
likelihood_info_DELETE(&(x->info));
UNCONDFREE(*S);
}
}
/*
likelihood_storage * likelihood_CREATEXXX(int n) { /// obsolete
likelihood_storage * x =
(likelihood_storage*) MALLOC(sizeof(likelihood_storage));
likelihood_NULL(x);
x->datasets = LIST_CREATE(n, LISTOF + REALSXP);
// x->withoutdet = LIST_CREATE(n, LISTOF + REALSXP);
if ((x->X = (double**) CALLOC(n, sizeof(double*))) == NULL ||
(x->XtX = (double*) CALLOC(n, sizeof(double))) == NULL ||
(x->XitXi = (double*) CALLOC(n, sizeof(double))) == NULL ||
(x->CinvXY = (double*) CALLOC(n, sizeof(double))) == NULL ||
(x->data_nas = (bool*) CALLOC(n, sizeof(bool))) == NULL)
err or("Memory allocation error for 'likelihood'");
x->sets = n;
}
*/
void covariate_NULL(covariate_storage *x) {
if (x == NULL) return;
x->loc = NULL;
x->matrix_err = MATRIX_NOT_CHECK_YET;
}
void covariate_DELETE(covariate_storage **S) {
covariate_storage *x = *S;
if (x != NULL) {
if (x->loc != NULL) LOC_DELETE(&(x->loc));
UNCONDFREE(*S);
}
}
void addModel(cov_model **pcov, int covnr, cov_model *calling, bool nullOK) {
cov_model *cov;
int i;
// printf("addmodel: %s %ld\n", CovList[covnr].name, (long int) *pcov);
cov = (cov_model*) MALLOC(sizeof(cov_model));
COV_NULL(cov);
//assert(cov->calling == NULL);
cov->nr = covnr;
if (*pcov != NULL) {
cov->nsub = 1;
cov->calling = (*pcov)->calling;
(*pcov)->calling = cov;
cov->sub[0] = *pcov;
for (i=0; i<=Forbidden; i++) {
// cov->user[i] = cov->sub[0]->user[i];
cov->pref[i] = cov->sub[0]->pref[i];
// printf(">> %s %s %d %d\n", NAME(cov), NAME(cov->sub[0]), i, cov->pref[i]);
}
}
// if (calling != NULL)
//printf("addmodel calling %s cov %s\n", NAME(calling), NAME(cov));
if (calling != NULL) cov->calling = calling;
else if (!nullOK && *pcov == NULL) {
PRINTF("Missing link for model '%s'.\n", NICK(cov));
BUG;
}
*pcov = cov;
}
void addModel(cov_model **pcov, int covnr) {
addModel(pcov, covnr, NULL, false);
}
void addModel(cov_model *pcov, int subnr, int covnr) {
assert(pcov != NULL);
bool newsub = pcov->sub[subnr] == NULL;
addModel(pcov->sub + subnr, covnr, pcov, false);
pcov->nsub += (int) newsub;
}
void addModelKappa(cov_model *pcov, int subnr, int covnr) {
assert(pcov != NULL);
addModel(pcov->kappasub + subnr, covnr, pcov, false);
}
void addModel(cov_model **pcov, int covnr, cov_model *calling) {
addModel(pcov, covnr, calling, false);
}
int setgrid(coord_type xgr, double *x, long lx, int spatialdim) {
if (lx!=3) {
// PRINTF("%d\n", lx);
// crash();
SERR("Problem with the coordinates (non-integer number of locations or non-positive step)")
}
int d;
unsigned long
totalBytes = sizeof(double) * lx * spatialdim; // nothing necessary for time
if (xgr[0] == NULL && (xgr[0] =(double*) MALLOC(totalBytes))==NULL)
return ERRORMEMORYALLOCATION;
//printf(" %ld %d, %d %d, \n", x, lx, spatialdim, totalBytes);
//printf("setgrid %d\n", totalBytes);
MEMCOPY(xgr[0], x, totalBytes);
// folgende Zeile nur beim ersten Mal zu setzen, aber egal
for (d=1; d<spatialdim; d++) {
xgr[d] = &(xgr[0][d * lx]);
if (xgr[d][XLENGTH] != (int) xgr[d][XLENGTH])
SERR2("grid length must be integer valued. Got %e in dimension %d.",
xgr[d][XLENGTH], d);
if (xgr[d][XLENGTH] < 1.0)
SERR2("grid length must be positive. Got %e in dimension %d.",
xgr[d][XLENGTH], d);
}
for (; d<MAXSIMUDIM; d++) xgr[d]=NULL;
/*
if (GLOBAL.internal.examples_reduced) {
for (d=0; d<spatialdim; d++) {
if (xgr[d][XLENGTH] > GLOBAL.internal.examples_reduced) {
warning("the size of the example has been reduced");
xgr[d][XLENGTH] = GLOBAL.internal.examples_reduced;
}
}
}
*/
return NOERROR;
}
int partial_loc_set(location_type *loc, double *x, double *y,
long lx, long ly, bool dist, int xdimOZ, double *T,
bool grid, bool cpy){
int d, err;
unsigned long totalBytes;
if (lx >= MAXINT || ly >= MAXINT) SERR("too many points");
STRCPY(ERROR_LOC, "setting the locations:");
if (((loc->x != NULL) && ((loc->y == NULL) xor (ly==0))) ||
((loc->xgr[0] != NULL) && ((loc->ygr[0] == NULL) xor (ly==0)))) {
// 28783440 51590192 0; 0 0 0 dist=0
//printf("%ld %ld %d; %ld %ld %d dist=%d\n", (long int) loc->x, (long int)loc->y, ly, (long int) loc->xgr[0], (long int) loc->ygr[0], ly, loc->distances);
SERR("domain structure of the first and second call do not match");
}
assert(x != NULL);
/*
if (GLOBAL.internal.examples_reduced) {
if (lx > GLOBAL.internal.examples_reduced) {
// lx = GLOBAL.internal.examples_reduced;
warning("The size of the example has been reduced.");
}
if (ly > GLOBAL.internal.examples_reduced) {
ly = GLOBAL.internal.examples_reduced;
warning("The size of y-coordinates in the example has been reduced.");
}
}
*/
loc->xdimOZ = xdimOZ; // ohne Zeit !!
loc->lx = lx;
loc->ly = ly;
if (ly > 0) {
//crash();
if (dist) { SERR("distances are not allowed if y is given"); }
}
// printf("partial_loc_set grid=%d dim=%d %f %f %f\n %f %f %f\n", grid, xdimOZ + (T!=NULL), x[0],x[1],x[2],x[3],x[4],x[5]);
loc->grid = grid;
loc->distances = dist;
if (loc->delete_y && loc->y != loc->x) FREE(loc->y);
if (loc->delete_x) FREE(loc->x);
loc->delete_x = loc->delete_y = cpy;
if (lx == 0) return NOERROR;
if (grid) {
loc->delete_x = true;
if ((err = setgrid(loc->xgr, x, lx, loc->spatialdim)) !=NOERROR) return err;
if (ly>0) {
if (x == y) {
for(d=0; d<loc->spatialdim; d++) loc->ygr[d] = loc->xgr[d];
loc->delete_y = false;
} else {
if ((err = setgrid(loc->ygr, y, ly, loc->spatialdim)) !=NOERROR)
return err;
}
}
double len = 1.0;
for (d=0; d<loc->spatialdim; d++) len *= loc->xgr[d][XLENGTH];
assert(len == (int) len);
if (len < MAXINT) loc->totalpoints = loc->spatialtotalpoints = (int) len;
else SERR("too many locations");
}
else if (dist) {
// lx : anzahl der Puntke involviert, d.h. x muss die Laenge lx ( lx -1) / 2
// haben.
if (lx > 0) {
if (cpy) {
totalBytes = sizeof(double) * lx * (lx - 1) / 2 * xdimOZ;
if ((loc->x=(double*) MALLOC(totalBytes))==NULL){
return ERRORMEMORYALLOCATION;
}
MEMCOPY(loc->x, x, totalBytes);
} else {
loc->x = x;
}
}
loc->totalpoints = loc->spatialtotalpoints = lx;
// (int) (1e-9 + 0.5 * (1.0 + SQRT(1.0 + 8.0 * lx)));
// if (0.5 * (loc->totalpoints * (loc->totalpoints - 1.0)) != lx) {
//printf("tot=%d %d %d\n", loc->totalpoints, (int) ( 0.5 * (loc->totalpoints * (loc->totalpoints - 1.0))), (int) lx); assert(false);
// SERR("distances do not have expected size");
//}
}
else { // not grid, not distances
if (cpy) {
totalBytes = sizeof(double) * lx * loc->xdimOZ;
//
assert(x != NULL);
// printf("PARTIAL : ");
assert(loc != NULL && loc->x == NULL);
//printf("\ntotalbyets %d %ld %d\n", totalBytes, lx, loc->xdimOZ);
if ((loc->x=(double*) MALLOC(totalBytes)) == NULL){
return ERRORMEMORYALLOCATION;
}
assert(loc->x != NULL);
//
MEMCOPY(loc->x, x, totalBytes);
// for (int k=0; k<lx *loc->xdimOZ; k++) loc->x[k] = x[k];
// for (int k=0; k<10; printf("%f ", x[k++]));
if (loc->ly>0) {
if (x == y) {
loc->y = loc->x;
loc->delete_y = false;
} else {
totalBytes = sizeof(double) * ly * loc->xdimOZ;
//printf("totalbytes y %d %ld %d\n", totalBytes, ly, loc->xdimOZ);
if ((loc->y=(double*) MALLOC(totalBytes))==NULL){
return ERRORMEMORYALLOCATION;
}
MEMCOPY(loc->y, y, totalBytes);
}
}
} else {
loc->x = x;
loc->y = y;
}
//printf("getnset %d\n", lx);
loc->totalpoints = loc->spatialtotalpoints = lx;
// just to say that considering these values does not make sense
}
if ((loc->Time) xor (T!=NULL)) {
//cov_model *cov; crash(cov);
// AERR(1);
SERR("partial_loc: time mismatch");
}
if (loc->Time) {
MEMCOPY(loc->T, T, sizeof(double) * 3);
if (grid) {
loc->xgr[loc->spatialdim] = loc->T;
if (ly>0) loc->ygr[loc->spatialdim] = loc->T;
}
if (loc->T[XLENGTH] <= 0) {
//crash();
SERR1("The number of temporal points is not positive. Check the triple definition of 'T' in the man pages of '%s'.", CovList[SIMULATE].nick);
}
/* else if (GLOBAL.internal.examples_reduced &&
loc->T[XLENGTH] > GLOBAL.internal.examples_reduced) {
loc->T[XLENGTH] = GLOBAL.internal.examples_reduced;
warning("The length of the tiem component in the example has been reduced.");
}
*/
if ((double) loc->totalpoints * loc->T[XLENGTH] >= MAXINT)
SERR("too many space-time locations");
loc->totalpoints *= (int) loc->T[XLENGTH];
}
return NOERROR;
}
int loc_set(double *x, double *y, double *T,
int spatialdim, /* spatial dim only ! */
int xdimOZ,
long lx, long ly, bool Time, bool grid,
bool distances,
location_type **Loc) {
int err;
//unsigned long totalBytes;
// preference lists, distinguished by grid==true/false and dimension
// lists must end with Nothing!
//printf("%d %d\n", spatialdim, Time);
if (xdimOZ < spatialdim) {
if (distances) {
if (xdimOZ != 1) SERR("reduced dimension is not one");
} else {
SERR3("dim (%d) of 'x' does not fit the spatial dim (%d); Time=%d",
xdimOZ,spatialdim, Time);
}
} else if (xdimOZ > spatialdim) {
SERR3("mismatch of dimensions (xdim=%d > space=%d; Time=%d)",
xdimOZ, spatialdim, Time);
}
int len = *Loc == NULL ? 1 : (*Loc)->len;
if (*Loc != NULL && (*Loc)->lx > 0) {
BUG;
LOC_SINGLE_DELETE(Loc);
*Loc = (location_type*) MALLOC(sizeof(location_type));
}
location_type *loc = *Loc;
LOC_SINGLE_NULL(loc, len);
loc->timespacedim = spatialdim + (int) Time;
loc->spatialdim = spatialdim;
loc->Time = Time;
if (spatialdim<1 || loc->timespacedim>MAXSIMUDIM) return ERRORDIM;
assert(x != NULL);
if ((err = partial_loc_set(*Loc, x, y, lx, ly, distances, xdimOZ,
Time ? T : NULL,
grid, true)) != NOERROR) XERR(err);
return NOERROR;
}
int loc_set(double *x, double *y, double *T,
int spatialdim, /* spatial dim only ! */
int xdimOZ,
long lx, long ly, bool Time, bool grid,
bool distances,
cov_model *cov) {
int err,
store = GLOBAL.general.set;
GLOBAL.general.set = 0;
// print("XXXXXXXXXXx\n");
// pmi(cov, 1);
LOC_DELETE(&(cov->ownloc));
assert(cov->ownloc == NULL);
LOCLIST_CREATE(cov, 1); // locown
assert(cov->ownloc != NULL);
assert(PLoc(cov) != cov->prevloc);
// PMI(cov);
// print("%ld %ld; %ld %f,%f,%f lx=%d %d time=%d grid=%d, %d %ld\n", cov->prevloc, cov->ownloc, x, x[0], x[1], x[2], lx, ly, Time , grid, distances, PLoc(cov));
err = loc_set(x, y, T, spatialdim, xdimOZ, lx, ly, Time, grid,
distances, PLoc(cov));
// Errorhandling:
GLOBAL.general.set = store;
return err;
}
int loc_set(double *x, double *T,
int spatialdim, // spatial dim only !
int xdimOZ, // original !
long lx, bool Time, bool grid,
bool distances,
location_type **Loc) {
return loc_set(x, NULL, T, spatialdim, xdimOZ, lx, 0, Time, grid,
distances, Loc);
}
int loc_set(double *x, double *T,
int spatialdim, /* spatial dim only ! */
int xdimOZ, /* original ! */
long lx, bool Time, bool grid,
bool distances,
cov_model *cov) {
return loc_set(x, NULL, T, spatialdim, xdimOZ, lx, 0, Time, grid,
distances, cov);
}
location_type ** loc_set(SEXP xlist, bool distances_ok){
bool
listoflists = (TYPEOF(xlist) == VECSXP &&
TYPEOF(VECTOR_ELT(xlist, 0)) == VECSXP);
int lx, err,
spatialdim,
xdimOZ = -1,
sets = !listoflists ? 1 : length(xlist);
bool
Time = false,
distances = false;
location_type **loc;
// printf("sets %d \n", sets);
loc = LOCLIST_CREATE(sets);
for (int i=0; i<sets; i++) {
SEXP
set = listoflists ? VECTOR_ELT(xlist, i) : xlist,
xx = VECTOR_ELT(set, XLIST_X),
yy = VECTOR_ELT(set, XLIST_Y),
TT = VECTOR_ELT(set, XLIST_T);
//print("set_addr=%ld : len.xx=%d %d %d\n", set, length(xx), length(yy), length(TT));
bool
dist = LOGICAL(VECTOR_ELT(set, XLIST_DIST))[0],
ggrid = LOGICAL(VECTOR_ELT(set, XLIST_GRID))[0];
int
xxdimOZ = ggrid ? ncols(xx) : nrows(xx),
llx = ggrid ? 3 : ncols(xx),
lly = length(yy) == 0 ? 0 : ggrid ? 3 : ncols(yy);
//printf("%d %d\n", ncols(xx), nrows(xx));
if (i==0) {
xdimOZ = xxdimOZ;
spatialdim = INTEGER(VECTOR_ELT(set, XLIST_SPATIALDIM))[0];
Time = LOGICAL(VECTOR_ELT(set, XLIST_TIME))[0];
distances = dist;
} else {
if (xdimOZ != xxdimOZ ||
spatialdim != INTEGER(VECTOR_ELT(set, XLIST_SPATIALDIM))[0] ||
Time != LOGICAL(VECTOR_ELT(set, XLIST_TIME))[0] ||
distances != dist
) BUG;
}
if (distances) {
if (distances_ok){
lx = (int) (1e-9 + 0.5 * (1 + SQRT(1. + 8 * llx)));
if (llx != lx * (lx - 1) / 2)
ERR("distance length not of form 'n * (n - 1) / 2'");
} else {
NotProgrammedYet("currently: distances in simulations");
distances = false;
}
} else {
lx = llx;
if (xxdimOZ != xdimOZ)
ERR("dimensions of the coordinates are not the same for the different sets");
}
if ((err = loc_set(REAL(xx), REAL(yy), REAL(TT), spatialdim, xdimOZ,
lx, lly, Time, ggrid,
distances, loc + i)) != NOERROR) {
LOC_DELETE(&loc);
XERR(err);
}
}
return loc;
}
/*
int loc_setX(cov_model *cov, long totalpoints) {
location_type *loc = ?;
if (cov->ownloc == NULL) {
LOCLIST_CREATE(cov, 1);
Loc(cov)->delete_ x = false;
} else {
loc = Loc(cov);
if (loc->xgr[0] != NULL || loc->x != NULL) BUG;
}
loc->totalpoints = totalpoints;
return NOERROR;
}
*/
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=NULL;
int j, k, startidx, i;
if (aniso == NULL) {
*diag = *quasidiag = *separatelast = *semiseparatelast = true;
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)) {
/*
* * 0
* * 0
* * *
*/
int last = col * row - 1;
for (k=last - row + 1; k<last; k++)
if (!(*separatelast = aniso[k] == 0.0)) break;
}
if (!(*separatelast = *semiseparatelast)) {
/*
* * 0
* * 0
0 0 *
*/
int last = col * row - 1;
for (k=row - 1; k<last; k+=row)
if (!(*separatelast = aniso[k] == 0.0)) break;
}
UNCONDFREE(taken);
}
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)
int match;
if (currentNrCov==-1) { InitModelList(); }
if (!STRCMP(name, InternalName)) return MATCHESINTERNAL;
if ((match = Match(name, CovNickNames, currentNrCov)) >= 0) return match;
return Match(name, CovNames, currentNrCov);
}
void MultiDimRange(int set, cov_model *cov, double *natscale) {
int wave, i, redxdim, d, idx,
xdimprev = cov->xdimprev,
vdim = cov->vdim[0],
store = GLOBAL.general.set,
err = NOERROR;
double y, yold, x[MAXGETNATSCALE], threshold, natsc, factor, Sign,
newx, xsave, newy,
*dummy = NULL,
rel_threshold = 0.05;
bool domain = cov->domown ==XONLY;
// covfct cf=NULL;
// nonstat_covfct ncf=NULL;
assert(MAXGETNATSCALE <= MAXCOVDIM); // ZERO
redxdim = cov->xdimown;
GLOBAL.general.set = set;
if (redxdim > MAXGETNATSCALE) {
err = -1; goto ErrorHandling;
}
if ((dummy = (double*) MALLOC(sizeof(double) * vdim * vdim)) == NULL) {
err = -2; goto ErrorHandling;
}
if (cov->full_derivs < 0) { err=ERRORNOTDEFINED; goto ErrorHandling; }
if (domain) {
COV(ZERO, cov, dummy);
} else {
NONSTATCOV(ZERO, ZERO, cov, dummy);
}
threshold = rel_threshold * dummy[0];
for (d=0; d<redxdim; d++) {
wave = 0;
for (i=0; i<xdimprev; i++) x[i] = 0.0;
idx = (redxdim == xdimprev || d==0) ? d : xdimprev-1;
x[idx] = 1.0;
if (domain) COV(x, cov, dummy) else NONSTATCOV(ZERO, x, cov, dummy);
yold = dummy[0];
if (ISNAN(yold)) { err = -3; goto ErrorHandling;}
if (yold > threshold) {
factor = 2.0;
Sign = 1.0;
} else {
factor = 0.5;
Sign = -1.0;
}
double otherx;
x[idx] *= factor;
if (domain) COV(x, cov, dummy) else NONSTATCOV(ZERO, x, cov, dummy);
y = dummy[0];
while (Sign * (y - threshold) > 0) {
if (yold<y){
if (wave++>10) { err=ERRORWAVING; goto ErrorHandling; }
}
yold = y;
x[idx] *= factor;
if (x[idx]>1E30) { err=ERRORRESCALING; goto ErrorHandling; }
if (domain) COV(x, cov, dummy) else NONSTATCOV(ZERO, x, cov, dummy);
y = dummy[0];
}
otherx = x[idx] / factor;
for (i=0; i<3 /* good choice?? */ ;i++) {
if (y==yold) { err=ERRORWAVING; goto ErrorHandling; }
newx = x[idx] + (x[idx]-otherx)/(y-yold)*(threshold-y);
xsave = x[idx];
x[idx] = newx;
if (domain) COV(x, cov, dummy) else NONSTATCOV(ZERO, x, cov, dummy);
newy = dummy[0];
x[idx] = xsave;
if (Sign * (newy - threshold) > 0) {
otherx = newx;
yold = newy;
} else {
x[idx] = newx;
y = newy;
}
}
if (y==yold) { err=ERRORWAVING; goto ErrorHandling; }
natsc = 1.0 / ( x[idx] +
(x[idx]-otherx)/(y-yold)*(threshold-y) );
if (redxdim == xdimprev || d==0) {
natscale[d] = natsc;
} else {
int beg, end;
if (redxdim == 2) {
if (d==0) {
beg = 0;
end = xdimprev-1;
} else {
beg = end = xdimprev-1;
}
} else {
beg = 0;
end = xdimprev;
}
for (i=beg; i<end; natscale[i++] = natsc);
}
}
ErrorHandling:
FREE(dummy);
GLOBAL.general.set = store;
switch(err) {
case NOERROR : break;
case -1 :
ERR("dimension of x-coordinates too high to detect natural scaling.");
case -2 :
ERR("not enough memory when determining natural scaling.");
case -3 :
ERR("NA in model evaluation detected");
default: XERR(err);
}
}
void UserGetNatScaling(double *natscale) {
GetNaturalScaling(KEY[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 !! */
cov_fct *C = CovList + cov->nr; // nicht gatternr
*natscale = 0.0;
if (C->maxsub!=0) XERR(ERRORFAILED);
if (!is_any(ISOTROPIC, C) || cov->isoown != ISOTROPIC
|| C->domain != XONLY ||
!isPosDef(cov->typus)
|| C->vdim != SCALAR)
ERR("anisotropic function not allowed");
if (C->finiterange == true) {
*natscale = 1.0;
return;
}
if (C->inverse!=NULL) {
C->inverse(&GLOBAL.gauss.approx_zero, cov, natscale);
*natscale = 1.0 / *natscale;
if (ISNAN(*natscale) || *natscale != 0.0) {
return;
}
}
if (NS != NATSCALE_ORNUMERIC) XERR(ERRORRESCALING);
if ((C->cov)==nugget) XERR(ERRORRESCALING);
// 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 --
*/
assert(cov->xdimown == 1);
MultiDimRange(0, cov, natscale);
}
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 TaylorCopy(cov_model *to, cov_model *from) {
int i, j;
to->taylorN = from->taylorN;
to->tailN = from->tailN;
for(i=0; i<to->taylorN; i++) {
for (j=0; j<=TaylorPow; j++) to->taylor[i][j] = from->taylor[i][j];
}
for(i=0; i<to->tailN; i++) {
for (j=0; j<=TaylorExpPow; j++) to->tail[i][j] = from->tail[i][j];
}
}
void paramcpy(cov_model *to, cov_model *from,
bool freeing, // of all the parameters
bool force_allocating, // in "c ovcpy" notwendig
bool copy_lists, // die Unterlisten der LISTOF-Parameter
bool recursive, bool copy_mpp) {
cov_fct *C = CovList + from->nr; // nicht gatternr
double **pto = to->px,
**pfrom = from->px;
int i, v,
n = -1,
*to_col = to->ncol,
*to_row = to->nrow,
*from_col = from->ncol,
*from_row = from->nrow;
bool same_model = std::abs(to->nr - from->nr) <= 1 ||
(isDollar(to) && isDollar(from));
if (!same_model) {
BUG;
}
for (i=0; i<MAXPARAM; i++) {
if (pfrom[i] == NULL) {
continue;
}
if (freeing) {
PARAMFREE(to, i);
to_col[i] = from_col[i];
to_row[i] = from_row[i];
}
if (C->kappatype[i] >= LISTOF) {
int len = from_row[i];
listoftype *p = PARAMLIST(from, i);
//printf("paramcpy %d %d %d %ld %d\n", copy_lists, from->zaehler, to->zaehler, (long int) (pto + i), pto[i] == NULL);
if (copy_lists) {
listcpy((listoftype **) (pto + i), p, force_allocating);
} else {
listpt((listoftype **) (pto + i), p, len, C->kappatype[i],
force_allocating);
}
} else if (isRObject(C->kappatype[i])) {
n = sizeof(sexp_type);
if (pto[i] == NULL || force_allocating) pto[i] = (double*) MALLOC(n);
MEMCOPY(pto[i], pfrom[i], n);
((sexp_type *) pto[i])->Delete = false;
} else {
if (C->kappatype[i]==REALSXP) {
n = sizeof(double);
} else if (C->kappatype[i]==INTSXP) {
n = sizeof(int);
} else BUG;
n *= from_row[i] * from_col[i];
// printf("i=%d to=%s %s from=%s %s %d\n", i,
//NICK(to), CovList[to->nr].kappanames[i],
// NICK(from), CovList[from->nr].kappanames[i], n);
if (pto[i] == NULL || force_allocating) pto[i] = (double*) MALLOC(n);
MEMCOPY(pto[i], pfrom[i], n);
}
}
if (copy_mpp) {
// PMI(to); PMI(from);
if (to->mpp.moments < 0 && alloc_mpp_M(to, from->mpp.moments)!=NOERROR)
ERR("error in allocating memory for Poisson point process data");
if (to->mpp.moments != from->mpp.moments) BUG;
//printf("size = %d\n", sizeof(mpp_properties));
assert(sizeof(mpp_properties) == 96);
mpp_properties *To = &(to->mpp), *From=&(from->mpp);
assert(To != NULL && From != NULL);
// To->sum_zhou_c = From->sum_zhou_c;
// To->sq_zhou_c = From->sq_zhou_c;
int vdim = from->vdim[0];
for (v=0; v<vdim; v++) To->maxheights[v] = From->maxheights[v];
To->unnormedmass = From->unnormedmass;
// To->zhou_c = From->zhou_c;
assert(To->mM != NULL && To->mMplus != NULL);
int nmP1 = To->moments + 1;
MEMCOPY(To->mM, From->mM, nmP1 * sizeof(double));
MEMCOPY(To->mMplus, From->mMplus, nmP1 * sizeof(double));
if (to->qlen != from->qlen) BUG;
if (from->qlen > 0) {
assert(to->q != NULL);
MEMCOPY(to->q, from->q, (to->qlen)* sizeof(double));
}
}
if (recursive) {
for (i=0; i<MAXSUB; i++) if (from->sub[i] != NULL) {
assert(to->sub[i] != NULL);
paramcpy(to->sub[i], from->sub[i], freeing,
force_allocating, copy_lists, recursive, copy_mpp);
}
}
//PMI(to, "to");
}
int covCpy(cov_model **localcov, bool sub, cov_model *cov, // err
location_type **prevloc, location_type **ownloc,
bool copy_lists, bool copy_randomparam,
bool allowCopyingInterface) {
assert(cov != NULL);
int i,
n = -1;
//cov_fct *C = CovList + cov->nr; // nicht gatternr
//PMI(cov);
if ((*localcov = (cov_model*) MALLOC(sizeof(cov_model)))==0)
return ERRORMEMORYALLOCATION;
cov_model *current = *localcov;
MEMCOPY(current, cov, sizeof(cov_model)); // replaces COV_NULL(*localcov);
COV_ALWAYS_NULL(current);
current->calling = NULL;
//printf("covCpy!\n");
paramcpy(current, cov, false, true, copy_lists, false, false);
if (cov->ownkappanames != NULL) {
int nkappas = CovList[cov->nr].kappas;
current->ownkappanames = (char**) CALLOC(nkappas, sizeof(char*));
for (i=0; i<nkappas; i++) {
if (cov->ownkappanames[i] != NULL) {
current->ownkappanames[i] =
(char*) MALLOC(sizeof(char) * (1 + STRLEN(cov->ownkappanames[i])));
STRCPY(current->ownkappanames[i], cov->ownkappanames[i]);
}
}
}
if (cov->q != NULL) {
n = sizeof(double) * current->qlen;
current->q = (double*) MALLOC(n); // QALLOC NOT APPROPRIATE
MEMCOPY(current->q, cov->q, n);
} else assert(current->qlen==0);
current->prevloc = ownloc != NULL ? ownloc :
cov->prevloc == prevloc ? prevloc : NULL;
if (current->prevloc == cov->prevloc && cov->calling==NULL) {
if (!isInterface(cov)) {
BUG;
}
if (!allowCopyingInterface) {PRINTF("\n\n***** unallowed copying ******\n"); BUG;}
}
for (i=0; i<MAXPARAM; i++) {
int err;
current->kappasub[i] = NULL;
if (cov->kappasub[i] == NULL || !copy_randomparam) continue;
err = covCpy(current->kappasub + i, true, cov->kappasub[i],
prevloc, ownloc, copy_lists, copy_randomparam, false);
if (err != NOERROR) return err;
current->kappasub[i]->calling = current;
}
if (sub) {
for (i=0; i<MAXSUB; i++) {
int err;
current->sub[i] = NULL;
if (cov->sub[i] == NULL) continue;
err = covCpy(current->sub + i, sub, cov->sub[i], prevloc, ownloc,
copy_lists, copy_randomparam, false);
if (err != NOERROR) return err;
current->sub[i]->calling = current;
}
} else {
for (i=0; i<MAXSUB; i++) current->sub[i] = NULL;
}
// PMI(*localcov, "end err=covCpy");
return NOERROR;
}
int covCpy(cov_model **localcov, cov_model *cov) { //err
bool cov2key = &(cov->key)==localcov;
// PMI(cov);
int
err = covCpy(localcov, true, cov,
cov->prevloc, NULL,
false, true, false);//err
// PMI(*localcov);
if (err == NOERROR)
(*localcov)->calling = cov2key || cov->calling==NULL ? cov : cov->calling;
// falls !cov2key && cov->calling == NULL; dann gibt es nur einen Verweis
// rueckwaerts! Muss aber sein, da sonst LOC_DELETE bei cov->calling==NULL
// versucht cov->prevloc zu loeschen. oder es muesste prevloc auf NULL
// gesetzt zu werden. Dies scheint eher contraproduktiv zu sein.
return err;
}
int covCpyWithoutRandomParam(cov_model **localcov, cov_model *cov) {//err
bool cov2key = &(cov->key)==localcov;
int
err = covCpy(localcov, true, cov, cov->prevloc, NULL, false, false, false);
if (err == NOERROR)
(*localcov)->calling = cov2key || cov->calling==NULL ? cov : cov->calling;
// falls !cov2key && cov->calling == NULL; dann gibt es nur einen Verweis
// rueckwaerts! Muss aber sein, da sonst LOC_DELETE bei cov->calling==NULL
// versucht cov->prevloc zu loeschen. oder es muesste prevloc auf NULL
// gesetzt zu werden. Dies scheint eher contraproduktiv zu sein.
return err;
}
int covCpy(cov_model **localcov, cov_model *cov, bool copy_lists) {//err
bool cov2key = &(cov->key)==localcov;
int
err = covCpy(localcov, true, cov, cov->prevloc, NULL, copy_lists, true,
false);
if (err == NOERROR)
(*localcov)->calling = cov2key || cov->calling==NULL ? cov : cov->calling;
// falls !cov2key && cov->calling == NULL; dann gibt es nur einen Verweis
// rueckwaerts! Muss aber sein, da sonst LOC_DELETE bei cov->calling==NULL
// versucht cov->prevloc zu loeschen. oder es muesste prevloc auf NULL
// gesetzt zu werden. Dies scheint eher contraproduktiv zu sein.
return err;
}
int covCpy(cov_model **localcov, cov_model *cov, //err
double *x, double *T, int spatialdim, int xdimOZ, long lx,
bool Time, bool grid, bool distances) {
bool cov2key = &(cov->key)==localcov;
int err;
location_type **loc = LOCLIST_CREATE(1);
err = loc_set(x, T, spatialdim, xdimOZ, lx, Time, grid, distances, loc);
if (err != NOERROR) return err;
if ((err = covCpy(localcov, true, cov, loc, NULL, false, true, false))
!= NOERROR) return err;
(*localcov)->prevloc = cov->prevloc;
(*localcov)->ownloc = loc;
(*localcov)->calling = cov2key || cov->calling==NULL ? cov : cov->calling;
return NOERROR;
}
cov_model *getRemote(cov_model *remotecov, cov_model *rmt, cov_model *target) {
cov_model *found;
int i;
if (rmt == target) return remotecov;
for (i=0; i<MAXPARAM; i++) {
if (rmt->kappasub[i] != NULL) {
if (remotecov->kappasub[i] == NULL) BUG;
found = getRemote(remotecov->kappasub[i], rmt->kappasub[i], target);
if (found != NULL) return found;
}
}
for (i=0; i<MAXSUB; i++) {
if (rmt->sub[i] != NULL) {
if (remotecov->sub[i] == NULL) BUG;
found = getRemote(remotecov->sub[i], rmt->sub[i], target);
if (found != NULL) return found;
}
}
return NULL;
}
void Ssetcpy(cov_model *localcov, cov_model *remotecov, cov_model *cov,
cov_model *rmt) {
int i;
if (cov->Sset != NULL) {
localcov->Sset = (set_storage*) MALLOC(sizeof(set_storage));
MEMCOPY(localcov->Sset, cov->Sset, sizeof(set_storage));
localcov->Sset->remote = getRemote(remotecov, rmt, cov->Sset->remote);
if (localcov->Sset->remote == NULL) BUG;
}
for (i=0; i<MAXPARAM; i++) {
if (cov->kappasub[i] != NULL) {
if (localcov->kappasub[i] == NULL) BUG;
Ssetcpy(localcov->kappasub[i], remotecov, cov->kappasub[i], rmt);
}
}
for (i=0; i<MAXSUB; i++) {
if (cov->sub[i] != NULL) {
if (localcov->sub[i] == NULL) BUG;
Ssetcpy(localcov->sub[i], remotecov, cov->sub[i], rmt);
}
}
}
int newmodel_covCpy(cov_model **localcov, int model, cov_model *cov, //err
double *x, double *y, double *T,
int spatialdim, /* spatial dim only ! */
int xdimOZ, long lx, long ly, bool Time, bool grid,
bool distances) {
Types type = CovList[model].Typi[0];
assert(type == InterfaceType && CovList[model].variants == 1); // anderes z.zt nicht benutzt
int i, err;
assert(*localcov == NULL);
addModel(localcov, model, NULL, true);
cov_model *neu = *localcov;
assert(neu->ownloc==NULL && neu->prevloc==NULL);
if (type == InterfaceType) neu->prevloc = LOCLIST_CREATE(1); // locown
else BUG;
assert(cov->prevloc != NULL);
assert(PLoc(cov) == cov->prevloc);
loc_set(x, y, T, spatialdim, xdimOZ, lx, ly, Time, grid, distances, neu);
if ((err = covCpy(neu->sub + 0, cov)) != NOERROR) return err;
// PMI(neu, "hier");
neu->sub[0]->calling = neu;
// PMI(neu);
//printf("newmodel %s %s\n", NICK(cov), TYPENAMES[CovList[cov->nr].Type]);
for (i=0; i<2; i++) {
if ((err = CHECK(neu, cov->tsdim, cov->xdimprev,
// CovList[cov->nr].Type,Martin:changed 27.11.13
//cov->typus,// auch nicht 20.5.14
type,
type == InterfaceType ? XONLY : cov->domprev,
type == InterfaceType ? CARTESIAN_COORD : cov->isoprev,
cov->vdim, ROLE_BASE)) != NOERROR) {
return err;
}
if (i==0 && (err = STRUCT(neu, NULL)) != NOERROR) return err;
}
return NOERROR;
}
int newmodel_covCpy(cov_model **localcov, int model, cov_model *cov) {//err
int err,
store = GLOBAL.general.set;
GLOBAL.general.set = 0;
location_type *loc = Loc(cov);
err = newmodel_covCpy(localcov, model, cov,
loc->grid ? loc->xgr[0] : loc->x,
loc->grid ? loc->ygr[0] : loc->y,
loc->grid ? loc->xgr[0] + 3 * loc->spatialdim : loc->T,
loc->spatialdim, loc->xdimOZ,
loc->grid ? 3 : loc->totalpoints,
loc->ly == 0 ? 0 : loc->grid ? 3 : loc->totalpoints,
loc->Time, loc->grid, loc->distances);
GLOBAL.general.set = store;
return err;
}
double *getAnisoMatrix(cov_model *cov, bool null_if_id, int *nrow, int *ncol) {
int
origdim = PrevLoc(cov)->timespacedim;
if (!isAnyDollar(cov) && null_if_id) { // probably not used anymore
*nrow = *ncol = origdim;
return NULL;
}
double *ani,
*aniso = P(DANISO),
a = PisNULL(DSCALE) ? 1.0 : 1.0 / P0(DSCALE);
int i, total,
*proj = PINT(DPROJ),
dimP1 = origdim + 1;
if (aniso != NULL) {
total = origdim * cov->ncol[DANISO];
int bytes = total * sizeof(double);
ani = (double *) MALLOC(bytes);
MEMCOPY(ani, aniso, bytes);
for (i=0; i<total; i++) {
ani[i] *= a;
}
*nrow = cov->nrow[DANISO];
*ncol = cov->ncol[DANISO];
} else if (proj != NULL) {
int projs = cov->nrow[DPROJ];
total = origdim * projs;
ani = (double *) CALLOC(total, sizeof(double));
for (i=0; i<projs; i++) {
ani[i * origdim + proj[i] - 1] = a;
}
*nrow = origdim;
*ncol = cov->nrow[DPROJ];
} else {
if (a == 1.0 && null_if_id) {
*nrow = *ncol = origdim;
return NULL;
}
total = origdim * origdim;
ani = (double *) CALLOC(total, sizeof(double));
for (i=0; i<total; i+=dimP1) ani[i] = a;
*nrow = *ncol = origdim;
}
return ani;
}
double *getAnisoMatrix(cov_model *cov, int *nrow, int *ncol) {
return getAnisoMatrix(cov, false, nrow, ncol);
}
double GetDiameter(location_type *loc, double *min, double *max,
double *center, bool docaniso) {
// 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 !
bool *j=NULL;
int d,
origdim = loc->timespacedim,
spatialdim = loc->spatialdim;
double distsq, dummy,
*lx=NULL, *sx=NULL,
diameter=0.0;
if (loc->grid) {
double
*origmin = (double*) MALLOC(origdim * sizeof(double)),
*origmax = (double*) MALLOC(origdim * sizeof(double)),
*origcenter = (double*) MALLOC(origdim * sizeof(double));
for (d=0; d<origdim; d++) {
if (loc->xgr[d][XSTEP] > 0) {
origmin[d] = loc->xgr[d][XSTART];
origmax[d] = loc->xgr[d][XSTART] + loc->xgr[d][XSTEP] *
(loc->xgr[d][XLENGTH] - 1.0);
} else {
origmin[d] = loc->xgr[d][XSTART] + loc->xgr[d][XSTEP] *
(double) (loc->xgr[d][XLENGTH] - 1.0);
origmax[d] = loc->xgr[d][XSTART];
}
origcenter[d] = 0.5 * (origmin[d] + origmax[d]);
}
if (!docaniso || loc->caniso == NULL) {
for (d=0; d<origdim; d++) {
center[d] = origcenter[d];
min[d] = origmin[d];
max[d] = origmax[d];
dummy = max[d] - min[d];
diameter += dummy * dummy;
}
} else { // caniso != NULL
j = (bool*) MALLOC( (origdim + 1) * sizeof(double));
lx = (double*) MALLOC( origdim * sizeof(double));
sx = (double*) MALLOC(spatialdim * sizeof(double));
xA(origcenter, loc->caniso, 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] = RF_INF;
max[d] = RF_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];
xA(lx, loc->caniso, 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;
}
if (distsq > diameter) diameter = distsq;
}
UNCONDFREE(j);
UNCONDFREE(lx);
UNCONDFREE(sx);
}
UNCONDFREE(origmin);
UNCONDFREE(origmax);
UNCONDFREE(origcenter);
} else { // not loc->grid
if (loc->caniso != NULL) BUG;
double *xx=loc->x;
int i,
endfor = loc->spatialtotalpoints * spatialdim;
for (d=0; d<spatialdim; 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[XSTEP] > 0) {
min[d] = loc->T[XSTART];
max[d] = loc->T[XSTART] + loc->T[XSTEP] * (loc->T[XLENGTH] - 1.0);
} else {
min[d] = loc->T[XSTART] + loc->T[XSTEP] * (loc->T[XLENGTH] - 1.0);
max[d] = loc->T[XSTART];
}
}
for (diameter=0.0, d=0; d<origdim; d++) {
center[d] = 0.5 * (max[d] + min[d]);
dummy = max[d] - min[d];
diameter += dummy * dummy;
}
} // !simugrid
return 2.0 * SQRT(diameter);
}
double GetDiameter(location_type *loc, double *min, double *max,
double *center) {
return GetDiameter(loc, min, max, center, true);
}
double GetDiameter(location_type *loc) {
double diam,
*dummymin=NULL, *dummymax=NULL, *dummycenter=NULL;
int origdim = loc->timespacedim;
dummymin = (double*) MALLOC(origdim * sizeof(double));
dummymax = (double*) MALLOC(origdim * sizeof(double));
dummycenter = (double*) MALLOC(origdim * sizeof(double));
diam = GetDiameter(loc, dummymin, dummymax, dummycenter, true);
UNCONDFREE(dummymin);
UNCONDFREE(dummymax);
UNCONDFREE(dummycenter);
return diam;
}
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) {
unsigned long i,ii,j,jj,l,ll,min, m=1;
if (HOMEMADE_NICEFFT) {
int f[F_NUMBERS1]={2,3,5};
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 expandgrid(coord_type xgr, double **xx, double* aniso,
int olddim,
int nrow, // matrix size of old dim
int ncol // new dim
){
double *x=NULL, * y=NULL; /* current point within grid, but without
anisotropy transformation */
int
*yi=NULL, /* counter for the current position in the grid */
dimM1 = olddim - 1;
long pts, w, k, total, n, i, d;
assert(olddim <= nrow);
if (aniso == NULL && olddim != ncol) BUG;
for (pts=1, i=0; i<olddim; i++) pts *= (long int) xgr[i][XLENGTH];
total = ncol * pts;
x = *xx = (double*) MALLOC(sizeof(double) * total);
y = (double*) MALLOC(olddim * sizeof(double));
yi = (int*) MALLOC(olddim * sizeof(int));
for (w=0; w<olddim; w++) {y[w]=xgr[w][XSTART]; yi[w]=0;}
for (k=0; k<total; ){
if (aniso==NULL) {
for (d=0; d<ncol; d++, k++) x[k] = y[d];
} else {
for (n=d=0; d<ncol; d++, k++, n+=nrow - olddim) {
x[k] = 0.0;
for(w=0; w<olddim; w++) {
x[k] += aniso[n++] * y[w];
}
}
}
i = 0;
(yi[i])++;
y[i] += xgr[i][XSTEP];
while(yi[i]>=xgr[i][XLENGTH]) {
yi[i]=0;
y[i] = xgr[i][XSTART];
if (i<dimM1) {
i++;
(yi[i])++;
y[i] += xgr[i][XSTEP];
} else {
assert(k==total);
}
}
}
UNCONDFREE(y);
UNCONDFREE(yi);
}
void grid2grid(coord_type xgr, double **grani, double *aniso, int origdim, int dim) {
// diag (+permutations) is assumed, so each row and each col has at most one non-zero
// component
double *pgr, *A;
int d, i,
origdimM1 = origdim - 1;
// print("\n\nALLOCATED %d\n", 3 * dim);
(*grani) = (double *) MALLOC(sizeof(double) * 3 * dim);
pgr = *grani;
if (aniso == NULL) {
for(d=0; d<dim; d++) {
for (i=0; i<3; i++, pgr++) {
*pgr = xgr[d][i];
}
}
} else {
for (d=0; d<dim; d++, pgr += 3) {
//printf("%d %d\n", d, origdim);
A = aniso + d * origdim;
// printf("%f %f %d\n", aniso[0], A[0], origdimM1);
for (i=0; i<origdimM1; i++, A++) if (*A != 0.0) break;
//printf("A=%f %d \n", A[0], i);
//printf("X=%f\n", xgr[i][XSTART]);
pgr[XSTART] = xgr[i][XSTART] * *A;
pgr[XSTEP] = xgr[i][XSTEP] * *A;
pgr[XLENGTH] = xgr[i][XLENGTH];
}
}
}
void xtime2x(double *x, int nx, double *T,
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,
timelen = (int) T[XLENGTH],
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,
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,
timelen = T[XLENGTH],
nxspdim = nx * spatialdim;
if (aniso == NULL) {
assert(nrow == ncol);
xtime2x(x, nx, T, newx, nrow);
return;
}
z = *newx = (double*) MALLOC(sizeof(double) * ncol * nx * timelen);
for (k=j=0, t=T[XSTART]; j<timelen; 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];
}
z[k++] = dummy + aniso[n++] * t; //auf keinen Fall nach vorne holen
}
}
}
}
void x2x(double *x, int nx, double **newx,
double *aniso, int physical_nrow, 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;
// printf("x2x %d x %d %d %d\n", nrow, ncol, nxnrow, nrow); BUG;
z = *newx = (double*) MALLOC(sizeof(double) * ncol * nx);
if (aniso == NULL) {
assert(nrow == ncol);
MEMCOPY(z, x, sizeof(double) * nx * ncol);
} else {
for (k=i=0; i<nxnrow; i+=nrow) {
endfor = i + nrow;
for (d=0; d<ncol; d++) {
n = physical_nrow * d;
dummy = 0.0;
for(w=i; w<endfor; w++) {
// printf(" n=%d %f\n", n, aniso[n]);
// printf(" w=%d %f %d\n", w, y[w], endfor);
dummy += aniso[n++] * y[w];
}
z[k++] = dummy;
}
}
}
}
bool isMiso(matrix_type type) {
return type == TypeMiso;
}
bool isMproj(matrix_type type) {
assert(TypeMtimesepproj == 2);
return type == TypeMproj || type <= TypeMtimesepproj;
}
bool isMdiag(matrix_type type) {
return type <= TypeMdiag;
}
bool isMtimesep(matrix_type type) {
assert(TypeMtimesep == 3);
return type <= TypeMtimesep;
}
matrix_type Type(double *M, int nrow, int ncol) {
// see also analyse_matrix: can it be unified ???
matrix_type type = TypeMiso; // default
// double elmt;
int i, j, k,
endfor = nrow * ncol;
double *m = M;
if (m==NULL || (ncol == 1 && nrow == 1)) {
assert(ncol == nrow);
return type;
}
if (ncol > nrow) {
for (i=ncol * ncol; i < endfor; i++)
if (m[i]!= 0.0) return TypeMany;
ncol = nrow; // !!
}
for (k=0; k<ncol; ) {
for (i=0; i<nrow; i++) if (!R_FINITE(m[i]) || m[i] != 0.0) break;
for (j=i+1; j<nrow; j++)
if (!R_FINITE(m[j]) || m[j] != 0.0) goto TimeInvestigation;
if (k == ncol - 1) {
/* // nur fuer strikte projektion auf kleinere Koordinaten
if (type == TypeMproj) {
int ii, kk;
for (ii=0; ii<nrow; ii++) {
double sum = 0.0;
for (kk=ii; kk<endfor; kk+= nrow) sum += M[kk];
printf("sum %d %f %d\n", ii, sum, type); //
if (sum != 1.0 && sum != 0.0) {
type = TypeMany;
break;
}
}
}
*/
if (i != k && (!R_FINITE(m[i]) || m[i] != 0.0)) {
return (type == TypeMproj || k==0 || type == TypeMiso)
// (type == TypeMiso && M[0] == 1.0))
&& R_FINITE(m[i]) //&& m[i] == 1.0
? TypeMproj
: TypeMany;
} else {
// korrekt ???
return type == TypeMproj
? (R_FINITE(m[i]) //&& (m[i] == 0.0 || m[i] == 1.0)
? TypeMtimesepproj : TypeMany)
: type == TypeMiso && (!R_FINITE(m[i]) || m[i] != M[0])
? TypeMdiag : type;
}
}
// k < ncol - 1; type != timesep and != timesepproj garantied
matrix_type newtype;
if ( i!=k && m[i] != 0.0) {
//if (m[i] != 1.0) goto TimeInvestigation;
newtype = TypeMproj;
} else {
newtype = R_FINITE(m[i]) && m[i]==M[0] ? TypeMiso : TypeMdiag;
}
if (newtype > type) type = newtype;
k++;
m += nrow;
continue;
TimeInvestigation:
if (k == ncol - 1) return TypeMany;
type = TypeMtimesep;
k = ncol - 1;
m = M + k * nrow;
}
return type;
}
void TransformLocExt(cov_model *cov, bool timesep,
int gridexpand, // false, true, GRIDEXPAND_AVOID
bool same_nr_of_points,
double **grani, double **SpaceTime,
double **caniso, int *Nrow, int *Ncol,//caniso obsolete
bool *Time, bool *grid, int *newdim, bool takeX,
bool involvedollar) {
// this fctn transforms the coordinates according to the anisotropy matrix
// tbm TransformLoc-Aufruf im prinzip kritisch
location_type *loc = Loc(cov);
bool isdollar = isAnyDollar(cov) && involvedollar;
matrix_type type;
int
nrow = -1,
ncol = -1,
origdim = loc->caniso == NULL ? loc->timespacedim : loc->cani_nrow,
dim = (isdollar ? (!PisNULL(DANISO) ? cov->ncol[DANISO] :
!PisNULL(DPROJ) ? cov->nrow[DPROJ] : origdim )
: origdim);
double *aniso,
*T = loc->T, *x = takeX ? loc->x : loc->y;
coord_type
*xgr= takeX ? &(loc->xgr) : &(loc->ygr);
if (x==NULL && (*xgr)[0] ==NULL) ERR("locations are all NULL");
*newdim = dim;
// printf("dim %d %d %d\n", dim , origdim, loc->spatialdim);
*caniso = NULL;
*Nrow = *Ncol = -1;
*grani = NULL;
*SpaceTime = NULL;
if (isdollar) {
aniso = getAnisoMatrix(cov, true, &nrow, &ncol);
} else {
aniso = NULL;
nrow = ncol = PrevLoc(cov)->timespacedim;
}
if (loc->caniso != NULL) {
if (aniso == NULL) {
unsigned long bytes = sizeof(double) * loc->cani_nrow * loc->cani_ncol;
aniso =(double*) MALLOC(bytes);
MEMCOPY(aniso, loc->caniso, bytes);
nrow = loc->cani_nrow;
ncol = loc->cani_ncol;
} else {
double *aniso_old = aniso;
assert(loc->cani_ncol == nrow);
aniso = matrixmult(loc->caniso, aniso_old, loc->cani_nrow, nrow, ncol);
nrow = loc->cani_nrow;
UNCONDFREE(aniso_old);
}
}
type = aniso == NULL ? TypeMiso : Type(aniso, origdim, dim);
// printf("type = %d\n", type);
*Time = loc->Time;
*grid = false;
assert(origdim == nrow);
assert(dim == ncol);
if (loc->grid) {
assert(*xgr != NULL);
if (gridexpand==true || (gridexpand==GRIDEXPAND_AVOID && !isMproj(type))) {
if (timesep && isMtimesep(type) && *Time) {
// space
//print("%ld %d %d;\n", aniso, nrow, ncol);
// printf("aniso=%f %f %f\n %f %f %f\n %f %f %f\n; %d %d\n", aniso[0], aniso[1],aniso[2], aniso[3], aniso[4], aniso[5], aniso[6], aniso[7], aniso[8], nrow, ncol);
expandgrid(*xgr, SpaceTime, aniso, nrow - 1, nrow, ncol - 1);
// time
grid2grid(*xgr + loc->spatialdim, grani,
aniso == NULL ? NULL : aniso + nrow * (ncol-1),
1, 1);
} else {
*Time = false;// time is also expanded, if given
expandgrid(*xgr, SpaceTime, aniso, nrow, nrow, ncol);
}
} else if (isMproj(type) && (!same_nr_of_points || ncol==nrow)) {
// grid wird multipliziert und/oder umsortiert
grid2grid(*xgr, grani, aniso, nrow, ncol);
*grid = true;
*Time = false; // no need to have time extra
} else { // !gridexpand, !isMproj, but still grid
// z.B. plusmalS.cc falls TBM_INTERN
// nur aniso auf grid multipliziert
int i,d,k;
(*grani) = (double *) MALLOC(sizeof(double) * 3 * origdim);
for (k=d=0; d<origdim; d++) {
for (i=0; i<3; i++) {
// printf(" !!!!!!! xgr %d %d %d %f\n", takeX, d, i,(*xgr)[d][i]);
(*grani)[k++] = (*xgr)[d][i];
}
}
*grid = true;
*caniso = aniso;
*Nrow = nrow;
*Ncol = ncol;
return; // hier rausgehen!
}
} else { // nogrid
//APMI(cov->calling->calling->calling->calling->calling->calling);
// APMI(cov);
// printf("Time = %d %d %d type=%d %d\n", *Time, timesep, isMtimesep(type), type,TypeMtimesepproj);
// printf("aniso %f %f %f\n", aniso[0], aniso[1], aniso[2]);
//if (dim==1) BUG;
if (! *Time) { // no grid no time
x2x(x, loc->spatialtotalpoints, SpaceTime, aniso, nrow, nrow, ncol);
} else if (timesep && isMtimesep(type)) { // no grid, but timesep
//PMI(cov);
//printf("%d %d %d\n", same_nr_of_points && ncol, nrow);
if (same_nr_of_points && ncol!=nrow) { // do not reduce
x2x(x, loc->spatialtotalpoints, SpaceTime, aniso, nrow, nrow-1, ncol-1);
grid2grid(&T, grani, aniso==NULL ? NULL : aniso + nrow*ncol - 1, 1, 1);
} else if (aniso != NULL && aniso[nrow * ncol - 1]==0.0) {//no time comp
x2x(x, loc->spatialtotalpoints, SpaceTime, aniso, nrow, nrow - 1, ncol);
*Time = false;
} else { // both time and space are left, and time is separated
x2x(x, loc->spatialtotalpoints, SpaceTime, aniso, nrow, nrow - 1,
ncol - 1);
grid2grid(&T, grani, aniso==NULL ? NULL : aniso + nrow*ncol - 1, 1, 1);
}
} else if (ncol == 1 && type == TypeMproj && aniso[nrow - 1] &&
!same_nr_of_points) { //only time component left
assert(nrow != ncol);
*Time = false;
*grid = true;
grid2grid(&T, grani, aniso==NULL ? NULL : aniso + nrow*ncol - 1, 1, 1);
} else { // no grid, but time, no timesep || not isMtimesep(type)
xtime2x(x, loc->spatialtotalpoints, T, SpaceTime, aniso, nrow, ncol);
*Time = false;
}
}
FREE(aniso);
// printf("NCOL %d\n", *ncol);
}
void TransformCovLoc(cov_model *cov, bool timesep, int gridexpand,
bool same_nr_of_points, bool involvedollar) {
location_type *loc = PrevLoc(cov); // transform2nogrid
assert(loc != NULL);
bool Time, grid;
int err,
spacedim=-1,
nrow=-1,
ncol=-1;
double *xgr=NULL,
*x = NULL,
*caniso = NULL;
//printf("%s %s %d\n", NAME(cov), NAME(cov->sub[0]), Loc(cov)->spatialdim);
if ((loc->y != NULL && loc->y != loc->x) ||
(loc->ygr[0] != NULL && loc->ygr[0] != loc->xgr[0])) {
ERR("unexpected y coordinates");
}
assert(cov->prevloc != NULL);
TransformLocExt(cov, timesep, gridexpand, same_nr_of_points, &xgr, &x,
&caniso, &nrow, &ncol, &Time, &grid, &spacedim, true,
involvedollar);
if (Time) spacedim--;
assert(cov->ownloc == NULL);
if (spacedim > 0) {
// pmi(cov,0); printf("%d %d time=%d %d\n", grid, spacedim, Time, loc->totalpoints);
err = loc_set(grid ? xgr : x,
grid ? xgr + 3 * spacedim : xgr,
spacedim,
spacedim,
grid ? 3 : loc->spatialtotalpoints,
Time, grid, false, cov);
} else {
assert(Time);
err = loc_set(xgr, NULL, 1, 1, 3, false, true, false, cov);
}
// falls not gridexpand und nicht diag bzw. proj
location_type *ownloc = Loc(cov);
assert(ownloc->caniso == NULL);
ownloc->caniso = caniso;
ownloc->cani_nrow = nrow;
ownloc->cani_ncol = ncol;
caniso = NULL;
FREE(x);
FREE(xgr);
if (err != NOERROR) ERR("error when transforming to no grid");
}
void TransformLoc(cov_model *cov, bool timesep, int gridexpand,
bool involvedollar) {
TransformCovLoc(cov, timesep, gridexpand, true, involvedollar);
}
void TransformLocReduce(cov_model *cov, bool timesep, int gridexpand,
bool involvedollar) {
TransformCovLoc(cov, timesep, gridexpand, false, involvedollar);
}
int TransformLoc(cov_model *cov, double **xx, bool involvedollar) {
bool Time, grid;
int newdim, nrow, ncol;
double *caniso = NULL,
*SpaceTime = NULL;
TransformLocExt(cov, false, true, true, &SpaceTime, xx,
&caniso, &nrow, &ncol, &Time, &grid, &newdim, true,
involvedollar);
assert(caniso == NULL && SpaceTime ==NULL);
return newdim;
}
int TransformLoc(cov_model *cov, double **xx, double **yy,
bool involvedollar) {
location_type *loc = Loc(cov);
bool Time, grid;
int newdim, nrow, ncol;
double *caniso = NULL,
*SpaceTime = NULL;
TransformLocExt(cov, false, true, true, &SpaceTime, xx, &caniso,&nrow, &ncol,
&Time, &grid, &newdim, true, involvedollar);
assert(caniso == NULL && SpaceTime ==NULL);
if (loc->y == NULL && loc->ygr[0] == NULL) *yy = NULL;
else TransformLocExt(cov, false, true, true, &SpaceTime, yy, &caniso, &nrow,
&ncol, &Time, &grid, &newdim, false, involvedollar);
assert(caniso == NULL && SpaceTime ==NULL);
return newdim;
}
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;
}
bool TypeConsistency(Types requiredtype, Types deliveredtype) {
if (deliveredtype == UndefinedType) BUG;
switch(requiredtype) {
case TcfType: return isTcf(deliveredtype);
case PosDefType: return isPosDef(deliveredtype);
case VariogramType: return isVariogram(deliveredtype);
case NegDefType: return isNegDef(deliveredtype);
case ProcessType : return (isProcess(deliveredtype) ||
isTrend(deliveredtype));
case GaussMethodType:return isGaussMethod(deliveredtype);
case BrMethodType: return isBRuserProcess(deliveredtype);
case PointShapeType: return isPointShape(deliveredtype);
case RandomType : return isRandom(deliveredtype);
case ShapeType : return isShape(deliveredtype);
case TrendType : return isTrend(deliveredtype);
case InterfaceType: return isInterface(deliveredtype);
case RandomOrShapeType : BUG; // sollte nur bei den Argumenten auftauchen
case OtherType : return isOther(deliveredtype);
default : BUG;
}
BUG;
return FALSE;
}
int TypeConsistency(Types requiredtype, cov_model *cov, int depth) {
// if (cov == NULL) crash();
assert(cov != NULL);
cov_fct *C = CovList + cov->nr; // nicht gatternr
// printf("entering TC %s requiredtype=%s undefined!\n", NAME(cov), TYPENAMES[requiredtype]);
if (isUndefined(C)) {
assert(C->TypeFct != NULL);
// printf(" undefined!\n");
int tc = C->TypeFct(requiredtype, cov, depth);
// printf("tc=%d %s\n", tc, NAME(cov));
return tc;
}
int i;
for (i=0; i<C->variants; i++) {
// printf("TC %d %s: %d %s %s %d %d \n", i, NAME(cov), TypeConsistency(requiredtype, C->Typi[i]), TYPENAMES[requiredtype], TYPENAMES[C->Typi[i]], depth, depth<=0 );
if (TypeConsistency(requiredtype, C->Typi[i])
&& (depth <= 0 || atleastSpecialised(cov->isoown, C->Isotropy[i]))
) {
// printf("ok ? %d \n", i +1);
return i + 1;
}
}
// printf("not ok %s\n", NAME(cov));
return false;
}
int get_internal_ranges(cov_model *cov, cov_model *min, cov_model *max,
cov_model *pmin, cov_model *pmax,
cov_model *openmin, cov_model *openmax) {
cov_fct *C = CovList + cov->nr; // nicht gatternr
int
err = NOERROR, k = 0, i = 0, // compiler dummy values
kappas = C->kappas ;
range_type range;
rangefct getrange = C->range;
SEXPTYPE *type = C->kappatype;
// print("entry %s \n", NAME(cov));
//PMI(cov, "gir");
if (kappas > 0) {
getrange(cov, &range);
// PMI(cov);
assert (cov->maxdim < 0 || cov->maxdim >= cov->tsdim);
err = NOERROR;
for (i=0; i<kappas; i++) {
// printf(" i=%d %s\n", i, KNAME(i));
int len=cov->ncol[i] * cov->nrow[i];
double dmin, dmax, dpmin, dpmax, dopenmin, dopenmax,
value = RF_NA;
dpmin = range.pmin[i];
dpmax = range.pmax[i];
dmin = range.min[i];
dmax = range.max[i];
dopenmin = (double) range.openmin[i];
dopenmax = (double) range.openmax[i];
// unsigned int xxx = POW(2, 32) - 1;
// signed int yyy = (signed int) (double) xxx;
// print("%u %d NA=%d %f NA_LOG=%d '%s'\n", xxx, yyy, NA_INTEGER, (double) NA_INTEGER, NA_LOGICAL, NA_STRING); BUG;
if (type[i] == INTSXP) {
if (dmin < -MAXINT) {
#ifdef RANDOMFIELDS_DEBUGGING
PRINTF("%s <min %s\n", NICK(cov), KNAME(i));
BUG;
#else
dmin = (double) -MAXINT;
#endif
}
if (dmax > MAXINT) {
#ifdef RANDOMFIELDS_DEBUGGING
PRINTF("%s >max %s\n", NICK(cov), KNAME(i));
BUG;
#else
dmax = (double) MAXINT;
#endif
}
}
// printf("%s %d %d %d %d\n", NAME(cov), i, type[i], REALSXP, INTSXP);
for (k=0; k<len; k++) {
if (type[i] == REALSXP) {
value = P(i)[k];
PARAM(min, i)[k] = dmin;
PARAM(max, i)[k] = dmax;
PARAM(pmin, i)[k] = dpmin;
PARAM(pmax, i)[k] = dpmax;
PARAM(openmin, i)[k] = dopenmin;
PARAM(openmax, i)[k] = dopenmax;
} else if (type[i] == INTSXP) {
value = PINT(i)[k] == NA_INTEGER ? RF_NA : (double) PINT(i)[k];
PARAMINT(min, i)[k] = dmin;
PARAMINT(max, i)[k] = dmax;
PARAMINT(pmin, i)[k] = dpmin;
PARAMINT(pmax, i)[k] = dpmax;
PARAMINT(openmin, i)[k] = dopenmin;
PARAMINT(openmax, i)[k] = dopenmax;
} else if (type[i] == LISTOF + REALSXP) {
listoftype
*p = PARAMLIST(min, i);
if (p->deletelist) {
int j,
lenj = p->nrow[k] * p->ncol[k];
double
*qmin = PARAMLIST(min, i)->lpx[k],
*qmax = PARAMLIST(max, i)->lpx[k],
*ppmin = PARAMLIST(pmin, i)->lpx[k],
*ppmax = PARAMLIST(pmax, i)->lpx[k],
*omin = PARAMLIST(openmin, i)->lpx[k],
*omax = PARAMLIST(openmax, i)->lpx[k];
for (j=0; j<lenj; j++) {
qmin[j] = dmin;
qmax[j] = dmax;
ppmin[j] = dpmin;
ppmax[j] = dpmax;
omin[j] = dopenmin;
omax[j] = dopenmax;
}
} // elses list had not been copied
value = RF_NA;
// error cannot appear as range is (-infty, +infty)
} else if (isRObject(type[i])) {
continue;
} else return ERRORUNKOWNSXPTYPE;
if (ISNAN(value)) {
continue;
}
dmin = range.min[i];
dmax = range.max[i];
if (value < dmin
|| value > dmax
|| (range.openmin[i] && value == dmin)
|| (range.openmax[i] && value == dmax)
) {
SPRINTF(ERRORSTRING,
"value (%f) of '%s' in '%s' not within interval %s%f,%f%s",
value, C->kappanames[i], NICK(cov),
range.openmin[i] ? "(" : "[", dmin, dmax,
range.openmax[i] ? ")" : "]"
);
return ERRORM;
}
} // for k
} // for i in kappas;
} // if (kappa > 0)
for (i=0; i<MAXPARAM; i++) {
if (cov->kappasub[i] != NULL) {
err = get_internal_ranges(cov->kappasub[i],
min->kappasub[i], max->kappasub[i],
pmin->kappasub[i], pmax->kappasub[i],
openmin->kappasub[i], openmax->kappasub[i]);
if (err != NOERROR) return err;
}
}
for (i=0; i<MAXSUB; i++) {
if (cov->sub[i] != NULL) {
err = get_internal_ranges(cov->sub[i], min->sub[i], max->sub[i],
pmin->sub[i], pmax->sub[i],
openmin->sub[i], openmax->sub[i]);
if (err != NOERROR) return err;
}
}
return NOERROR;
}
void strround(double x, char *s){
if (x==RF_INF) SPRINTF(s, "Inf"); else
if (x==RF_NEGINF) SPRINTF(s, "-Inf"); else
if (x == FLOOR(x + 0.5)) SPRINTF(s, "%d", (int) x);
else SPRINTF(s, "%f", x);
}
void addmsg(double value, const char *Sign, double y, char *msg) {
char str1[30], str2[30];
if ( FABS(value - y) <= 1e-8 * y) {
SPRINTF(msg, "%12.12e %s %12.12e", value, Sign, y);
} else {
strround(value, str1);
strround(y, str2);
SPRINTF(msg, "%s %s %s", str1, Sign, str2);
}
}
void addone(char *str, double x) {
char str2[30];
strround(x, str2);
SPRINTF(str, "%s, %s", str, str2);
}
void addpair(char *str, double x, double y) {
if (x == FLOOR(x + 0.5) && y==floor(y + 0.5))
SPRINTF(str, "%s, (%d,%d)", str, (int) x, int(y));
else
SPRINTF(str, "%s, (%f,%f)", str, x, y);
}
int check_within_range(cov_model *cov, bool NAOK) {
// sets also maxdim and finiterange in cov !
cov_fct *C = CovList + cov->nr; //nicht gatternr
int len,
err = NOERROR,
k = 0,
i = 0, // compiler dummy values
kappas = C->kappas;
range_type range;
SEXPTYPE *type = C->kappatype;
rangefct getrange = C->range;
char Msg[255] = "";
double min, max,
value= RF_NA;
if (GLOBAL_UTILS->basic.skipchecks) return NOERROR;
getrange(cov, &range);
if (cov->maxdim >=0 && cov->maxdim < cov->tsdim) {
GERR2("Max. dimension is %d. Got %d", cov->maxdim, cov->tsdim);
}
for (i=0; i<kappas; i++) {
//printf("%s %d %s\n", NAME(cov), i, KNAME(i));
if (cov->kappasub[i] != NULL) continue;
if (!STRCMP(C->kappanames[i], FREEVARIABLE)) {
// i.e. equal
if (PisNULL(i)) continue;
}
if (type[i] >= LISTOF) {
// to simplify things -- otherwise in simu.cc
// code must also be changed.
assert(range.min[i] == RF_NEGINF && range.max[i]==RF_INF);
continue;
}
// full range !!
len = cov->ncol[i] * cov->nrow[i];
min = range.min[i];
max = range.max[i];
// print("%s i=%d [%f, %f] size=%d %s\n",
// C->name, i, range.min[i], range.max[i], len, C->kappanames[i]);
for (k=0; k<len; k++) {
if (type[i] == REALSXP) value = P(i)[k];
else if (type[i] == INTSXP)
value = PINT(i)[k] == NA_INTEGER ? RF_NA : (double) PINT(i)[k];
else if (isRObject(type[i])) continue;
else GERR2("%s [%d] is not finite", C->kappanames[i], k+1);
if (ISNAN(value)) {
if (NAOK) {
continue;
} else GERR2("%s[%d] is not finite.", C->kappanames[i], k+1);
}
// print("k=%d %f\n", k, value);
// print("range %d %d %d op=%d %f %f %f\n",
// kappas, i, k, range.openmin[i], value, min, max);
err = ERRORUNSPECIFIED;
if (range.openmin[i] && value <= min) {
addmsg(value, ">", min, Msg);
goto ErrorHandling;
} else if (value < min) {
addmsg(value, ">=", min, Msg);
goto ErrorHandling;
}
if (range.openmax[i] && value >= max) {
addmsg(value, "<", max, Msg);
goto ErrorHandling;
} else if (value > max) {
addmsg(value, "<=", max, Msg);
goto ErrorHandling;
}
err = NOERROR;
}
} // kappas
ErrorHandling:
// PMI(cov, "range");
if (err != NOERROR) {
if (PL>PL_ERRORS)
PRINTF("error in check range: %s kappa%d err=%d ('%s' does not hold.)\n",
C->name, i+1, err, Msg);
if (err == ERRORUNSPECIFIED)
SERR4("%s[%d] = %s does not hold (dim=%d).",
C->kappanames[i], k+1, Msg, cov->tsdim); // + return
}
return err;
}
int get_ranges(cov_model *cov, cov_model **min, cov_model **max,
cov_model **pmin, cov_model **pmax,
cov_model **openmin, cov_model **openmax) {
int err;
// returns a reliable value only in the very first
// entry of each parameter vector/matrix int err;
if ((err = covCpy(min, cov, true)) != NOERROR) return err;
if ((err = covCpy(max, cov, true)) != NOERROR) return err;
if ((err = covCpy(pmin, cov, true)) != NOERROR) return err;
if ((err = covCpy(pmax, cov, true)) != NOERROR) return err;
if ((err = covCpy(openmin, cov, true)) != NOERROR) return err;
if (( err = covCpy(openmax, cov, true)) != NOERROR) return err;
(*min)->calling = (*max)->calling = (*pmin)->calling = (*pmax)->calling =
(*openmin)->calling = (*openmax)->calling = cov;
// (*min)->calling = NULL; ja nicht auf NULL setzen, da sonst angenommen
// wird, dass prevloc ein Original ist
return get_internal_ranges(cov, *min, *max, *pmin, *pmax, *openmin, *openmax);
}
int FieldReturn(cov_model *cov) {
location_type *loc = Loc(cov);
if (cov->rf != NULL) {
//PMI(cov, "fieldreturn");
assert(cov->fieldreturn && cov->origrf);
if (cov->origrf) {
UNCONDFREE(cov->rf);
}
}
if ((cov->rf =
(double*) MALLOC(sizeof(double) * loc->totalpoints * cov->vdim[0]))
== NULL) return ERRORMEMORYALLOCATION;
cov->fieldreturn = cov->origrf = true;
return NOERROR;
}
void SetLoc2NewLoc(cov_model *cov, location_type **loc) {
int i,
maxsub = CovList[cov->nr].maxsub;
if (cov->ownloc != NULL) {
return;
}
for (i=0; i<MAXPARAM; i++) {
if (cov->kappasub[i] != NULL) SetLoc2NewLoc(cov->kappasub[i], loc);
}
cov->prevloc = loc;
for (i=0; i<maxsub; i++) {
if (cov->sub[i] != NULL) SetLoc2NewLoc(cov->sub[i], loc);
}
if (cov->key != NULL) SetLoc2NewLoc(cov->key, loc);
}
bool verysimple(cov_model *cov) {
assert(cov != NULL);
cov_fct *C = CovList + cov->nr; // kein gatternr, da isotrop
int i, k, total,
kappas = C->kappas;
for (k=0; k<kappas; k++) {
if (cov->kappasub[k] != NULL) return false;
total = cov->ncol[k] * cov->nrow[k];
if (C->kappatype[k] == REALSXP) {
for (i=0; i<total; i++)
if (ISNAN(P(k)[i]) || ISNA(P(k)[i])) return false;
} else if (C->kappatype[k] == INTSXP) {
for (i=0; i<total; i++) if (P(k)[i] == NA_INTEGER) return false;
} else return false;
}
return true;
}
