/*
 Authors
 Martin Schlather, martin.schlather@cu.lu

 library for unconditional simulation of stationary and isotropic random fields

 Copyright (C) 2001 -- 2004 Martin Schlather, 

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/


/*

  PRINTING LEVELS
  ===============
  error messages: 1
  forcation : 2
  minor tracing information : 3--5
  large debugging information: >10

*/

/*
  calculate the transformation of the points only once and store the result in
  a register (indexed by ncov) of key, not of s->. Advantages: points have to
  calculated only once for several tries; if zonal anisotropy then specific 
  methods may be called;  
*/

#include <math.h>  
#include <stdio.h>  
#include <stdlib.h>
#include <sys/timeb.h>
#include <assert.h>
#include <string.h>
#include "RFsimu.h"
#include <unistd.h>

void DeleteKey(int *keyNr)
{  
  key_type *key;
 
  int d, m;
  if (GENERAL_PRINTLEVEL>=4) { 
    PRINTF("deleting stored parameters of register %d ...\n", *keyNr);
  }
  if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {return;}
  key = &(KEY[*keyNr]);
  if (key->x[0]!=NULL) free(key->x[0]);
  for (d=0; d<MAXDIM; d++) {key->x[d]=NULL;}
  for (m=0; m<MAXCOV; m++) {
    if (key->destruct[m]!=NULL) key->destruct[m](&(key->S[m]));
    key->destruct[m]=NULL;
  }
  key->TrendModus = -1;
  if (key->TrendFunction!=NULL) free(key->TrendFunction); 
  key->TrendFunction=NULL;
  if (key->LinearTrend!=NULL) free(key->LinearTrend); 
  key->LinearTrend=NULL;
  if (key->destructX!=NULL) key->destructX(&(key->SX));
  key->destructX=NULL;
  key->active=false;
}

void DeleteAllKeys() {
  int i;
  for (i=0;i<MAXKEYS;i++) {DeleteKey(&i);}
}


void printkey(key_type *key) {
  int i,actparam,j;  
  char op_sign[][2]={"+","*"};

  PRINTF("\ngrid=%d, active=%d, anisotropy=%d, tbm_method=%d lx=%d sp_dim=%d,\ntotalpoints=%d, distr=%s Time=%d [%f %f %f]\ntimespacedim=%d traditional=%d compatible=%d mean=%f ncov=%d\n",	
	 key->grid, key->active, key->anisotropy,
	 key->tbm_method, key->lx, key->spatialdim, key->totalpoints,
	 DISTRNAMES[key->distribution], 
	 key->Time, key->T[0], key->T[1], key->T[2], 
	 key->timespacedim, key->traditional, key->compatible,
	 key->mean, key->ncov);
  PRINTF("\n");
  for (i=0; i<key->ncov; i++) {
    PRINTF("%d: meth=%s, cov=%d (%s) ^S=%d ^destruct=%d\n left=%d unimeth=%d (%s) [%d], param=", i,
	   METHODNAMES[key->method[i]],
	   key->covnr[i], CovList[key->covnr[i]].name,
	   key->S[i],key->destruct[i], 
	   key->left[i], 
	   key->unimeth[i], METHODNAMES[key->unimeth[i]], key->unimeth_used[i]
	   );
    for (j=0; j<key->totalparam; j++) PRINTF("%f,",key->param[i][j]);
    if (i < key->ncov - 1) PRINTF("\n%s\n",op_sign[key->op[i]]); 
    else PRINTF("\n\n");
  }
  //
  for (i=0; i<MAXDIM; i++) 
  PRINTF("%d : len=%d \n", i, key->length[i]);

  PRINTF("^x=%d ^y=%d ^z=%d ^SX=%d ^destructX=%d \n", 
	 key->x[0], key->x[1], key->x[2], key->SX, key->destructX
	 );

}

void printKEY(int *keyNr) {
  key_type *key;
  if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
    PRINTF("\nkeynr out of range!\n");
    return;
  }
  key = &(KEY[*keyNr]);
  PRINTF("\nNr=%d", *keyNr);
  printkey(key);
}


void printAllKeys() {
  int i;
  for (i=0;i<MAXKEYS;i++) {printKEY(&i);PRINTF("\n");}
}

int Transform2NoGrid(key_type *key, Real* param, int TrueDim, 
		     int *start_param, Real **xx)
{ 
  /* 
     this function transforms the coordinates according to the anisotropy
     matrix given in param (usually param is the param the first anisotropy
     matrix of the covariance product -- all the other anisotropy matrices
     have to multiple of the first, and this parameter is put into the
     covariance function as scale parameter 
     
     
     input  : key, param, TrueDim (got from GetTrueDim)
              start_param (vector of places where the parameters start)
     output : 
         **xx
           isotropic   : grid : NULL
                         arbitrary : x-coordinates multiplied with INVSCALE
           anisotropic : dim-reduced, x-coordinates multiplied with aniso matrix.
                         The index for the  temporal component runs the
                         slowest, if there is any

  */
  int i,j,k,d,w;
  unsigned long endfor, startfor, total;
  int dimM1;
  Real t, *x;
  dimM1 = key->timespacedim - 1;
  assert(xx!=NULL);
  if (key->anisotropy) {
    // used: TrueDim, start_param
    total = TrueDim * key->totalpoints;
    // total number of coordinates to store

    if ((x = *xx = (Real*) malloc(sizeof(Real) * total))==0)
      return ERRORMEMORYALLOCATION;
    /* determine points explicitly */
    if (key->Time && !key->grid) { // time component, no grid
      endfor = key->length[key->spatialdim];
      k=0;
      for (j=0, t=key->T[XSTART]; j<endfor; j++, t+=key->T[XSTEP])
	for (i=0; i<key->length[0]; i++)
	  for (d=0; d<TrueDim; d++, k++) {
 	    x[k] = 0.0;
	    for(w=0; w<key->spatialdim; w++)
	      x[k] += param[start_param[d]+w] * key->x[w][i];
	    x[k] += param[start_param[d]+w] * t;
	  }
    } else { 
      if (key->grid) {/* grid; with or without time component */
	Real y[MAXDIM]; /* current point within grid, but without
			   anisotropy transformation */
	int yi[MAXDIM]; /* counter for the current position in the grid */
	for (w=0; w<key->timespacedim; w++) {y[w]=key->x[w][XSTART]; yi[w]=0;}
	for (k=0; k<total; ){
	  for (d=0; d<TrueDim; d++, k++) {
	    x[k] = 0.0;
	    for(w=0; w<key->timespacedim; w++)
	      x[k] += param[start_param[d]+w] * y[w];
	  }
	  i = 0;
	  (yi[i])++;
	  y[i] += key->x[i][XSTEP];
	  while(yi[i]>=key->length[i]) {
	    yi[i]=0;
	    y[i] = key->x[i][XSTART];
	    if (i<dimM1) {
	      i++;
	      (yi[i])++;
	      y[i] += key->x[i][XSTEP];
	    } else {
	      assert(k==total);
	    }
	  }
	}
      } else {/* no grid, no time component */
	k = 0;
	for (i=0; i<key->length[0]; i++)
	  for (d=0; d<TrueDim; d++, k++) {
	    x[k] = 0.0;
	    for(w=0; w<key->spatialdim; w++)
	      x[k] += param[start_param[d]+w] * key->x[w][i];
	  }
      }
    }
  } else { /* isotropic, no time component */
    assert(fabs(param[INVSCALE] * param[SCALE] - 1.0)<EPSILON);
    if (key->grid) {    
      *xx = NULL;
    } else{
      if ((x = *xx = (Real*) malloc(sizeof(Real) * TrueDim * key->totalpoints)
	   )==0)
	// *xx is returned !!
	return ERRORMEMORYALLOCATION;      
      for (k=i=0; i<key->length[0]; i++) {
	for (j=0; j<key->spatialdim; j++) {
	  x[k++] = key->x[j][i] * param[INVSCALE];
	}
      }
    }
  }
  return 0;
}


int InternalGetGridSize(Real *x[MAXDIM], int *dim, int *lx)
{  
  int d;
  for (d=0; d<*dim; d++) {
    if ((x[d][XEND]<=x[d][XSTART]) || (x[d][XSTEP]<=0))
      return ERRORCOORDINATES;
    lx[d] = 1 + (int)((x[d][XEND]-x[d][XSTART])/x[d][XSTEP]); 
  }
  return NOERROR;
}


// obsolete?!
void GetGridSize(Real *x,Real *y,Real *z,int *dim,int *lx,int *ly,int *lz)
{ 
  // should now be consistent with `seq' of R
  // if not, please report!!!!!!!!!!!! 
  Real *xx[MAXDIM];
  int lxx[MAXDIM], d;
  xx[0]=x; xx[1]=y; xx[2]=z;
  if (InternalGetGridSize(xx,dim,lxx)) {
    *lx = *ly = *lz = 0;
  } else {
    *lx = lxx[0];  *ly = lxx[1];  *lz=lxx[2];
  }
}

#define COV_VARIO(FORMULA,FCT)\
  int v, aniso, k, j, endfor;\
  Real  var,  zz, zw, result, d;\
  cov_fct *cov;\
  zw = result = 0.0;\
  if (anisotropy) {\
    Real z[MAXDIM];\
    int cov_isotropy, ncovM1;\
    ncovM1 = ncov - 1;\
    for (v=0; v<ncov; v++) {\
      zw = var = 1.0;\
      for(;v<ncov;v++) {\
        cov = &(CovList[covnr[v]]);\
        cov_isotropy=cov->isotropic;\
	var *= param[v][VARIANCE];\
	for (k=0; k<dim; k++) z[k]=0.0;\
	for(aniso=ANISO, k=0; k<dim; k++){\
	  for (j=0; j<dim; j++, aniso++) {\
	    z[k] += param[v][aniso] * x[j];\
            }\
        }\
 	if (cov_isotropy==ANISOTROPIC) zw *= cov->FCT(z,param[v],dim);\
          /* derzeit is dim in cov->FCT bis auf assert-checks unbenutzt */\
	else {\
          /* the following definiton allows for calculating the covariance */\
          /* funcion for spatialdim=0 and SPACEISOTROPIC model, namely as a */\
          /* purely temporal model; however, in CheckAndRescale such kind of */\
          /* calculations are prevend as being TRIVIAL(generally, independent */\
          /* of being SPACEISOTROPIC or not) -- unclear whether relaxation in */\
          /* CheckAndRescale possible */\
	  zz = 0.0;\
	  endfor = dim - 1;\
	  for (j=0; j<endfor; j++) zz += z[j] * z[j];\
	  if (cov_isotropy==FULLISOTROPIC) zz += z[endfor] * z[endfor];\
	  else z[1]=z[endfor];\
	  z[0] = sqrt(zz);\
	  zw *= cov->FCT(z, param[v], 1 + (int) (cov_isotropy==SPACEISOTROPIC));\
	}\
	if ( (v<ncovM1) && (op[v]==0)) break;\
      }\
      result += FORMULA;\
    }\
  } else {\
    if (dim==1) d=x[0];\
    else {\
      for (d=0.0, j=0; j<dim; j++) {d += x[j] * x[j];}\
      d = sqrt(d);\
    }\
    for (v=0; v<ncov; v++) {\
      zw = var = 1.0;\
      for(;v<ncov;v++) {\
	zz = d * param[v][INVSCALE];\
	var *= param[v][VARIANCE];\
        cov = &(CovList[covnr[v]]); /* cov needed in FORMULA for Vario */\
	zw *= cov->FCT(&zz,param[v],1);\
	if (op[v]==0) break; /* v=ncov does not matter since stopped anyway */\
      }\
      result += FORMULA;\
    }\
  }\
 return result;


Real CovFct(Real *x, int dim, int *covnr, int *op,
  param_type param, int ncov, bool anisotropy) {
  COV_VARIO(var * zw, cov);
}

Real LocalCovFct(Real *x, int dim, int *covnr, int *op,
  param_type param, int ncov, bool anisotropy) {
  COV_VARIO(var * zw, cov_loc);
}

Real Vario(Real *x, int dim, int *covnr, int *op,
	    param_type param, int ncov, bool anisotropy){
  // Vario defines some kind of multiplication (namely as the mulitplication
  // of the respective covariance functions)
  // 
  // make sure that no genuine variogram model are used in multiplication terms
  //
  // first part says how to add models
  // cov.variogram checks whether last model in the block is a variogram
  // this may allow relaxations towards multiplications of varioagrams in
  // the sense that the corresponding covariance functions are multiplied
  // and then 1-block is used as value for variogram (only covariance functions
  // may appear in this case
  COV_VARIO(cov->variogram ? var * zw : var * (1.0 - zw), cov); // 30.1.02
}


// checks a *single* basic covariance function
// p and param must be of the same kind (Real)!!! 
int CheckAndRescale( int covnr, int naturalscaling, Real *p,
		    int timespacedim, int anisotropy, Real* param) 
{ 
  // IMPORTANT!: p is the parameter vector of the R interface !!
  // output param : equals p, except naturalscaling
  int actparam, error, endfor, i, j;
  Real newscale;

  if ((covnr>=currentNrCov) || (covnr<0)) return ERRORNOTDEFINED;
  if (p[VARIANCE]<0.0) return ERRORNEGATIVEVAR;
  
  actparam = KAPPA + CovList[covnr].kappas;
  if (anisotropy) {
    endfor = ANISO + timespacedim * timespacedim;
    for (j=actparam, i=ANISO; i<endfor; i++, j++) param[i] = p[j];     
  } else {
    if (CovList[covnr].isotropic!=FULLISOTROPIC) return ERRORANISOTROPIC;
    param[SCALE] = p[actparam];	
  }
  
  { // is dimension OK ? 
    int j, TrueDim, start_param[MAXDIM], index_dim[MAXDIM];
    bool no_last_comp;
    switch (CovList[covnr].isotropic) {
    case FULLISOTROPIC : 
      GetTrueDim(anisotropy, timespacedim, param,  &TrueDim, &no_last_comp, 
		 start_param, index_dim);
      break;
    case SPACEISOTROPIC :
      if (anisotropy) {
	GetTrueDim(anisotropy, timespacedim, param,  &TrueDim, &no_last_comp, 
		 start_param, index_dim);
	if (!no_last_comp) TrueDim--;
      } else {
	return ERRORNOTANISO;
      }
      break;
    case ANISOTROPIC :
      if (!anisotropy) return ERRORNOTANISO;
      TrueDim = timespacedim;
      break;
    default : assert(false);
    } 
    if (TrueDim==0) return ERRORTRIVIAL;
    if (CovList[covnr].method(TrueDim, false)==Forbidden)
      return ERRORCOVNOTALLOWED;
  }

  // specific parameter check for the covariance function
  if ((CovList[covnr].check!=NULL) && 
      ((error=(CovList[covnr].check(p, timespacedim, Nothing))))!=NOERROR)
    return error;
 
  // getnaturalscaling
  GetNaturalScaling(&covnr,&(p[KAPPA]),&naturalscaling,&newscale,&error);
  if (error!=NOERROR) return error;

  if (anisotropy) {
    newscale = 1.0 / newscale;
    for (j=actparam, i=ANISO; i<endfor; i++, j++) param[i] = p[j] * newscale;
  } else {
    param[SCALE] *= newscale;	
    if (CovList[covnr].cov==nugget && param[SCALE]==0) param[SCALE] = 1;
    param[INVSCALE] = 1.0 / param[SCALE];
    if (param[SCALE] <= 0.0) return ERRORNEGATIVESCALE;
  }
  for (i=0; i<actparam; i++) param[i] = p[i]; /* variance, kappas */

  return NOERROR;
}


#define FIRST_CHECK_COV_VARIO(ERROR)\
  if (currentNrCov==-1) InitModelList(); assert(CovList!=NULL);\
  if ((*logicaldim!=*xdim) && (*anisotropy || (*xdim!=1))) \
    {ERROR=ERRORDIMMISMATCH; goto ErrorHandling;}\
  if ((*xdim<1) || (*xdim>MAXDIM)) {ERROR=ERRORDIM; goto ErrorHandling;}\
  if (*anisotropy) { \
    pAdd = *xdim * *xdim + 1; /* VARIANCE */\
  } else {\
    pAdd = 2;\
  }\
  for (i=pAdd * *ncov, v=0; v<*ncov; v++) {\
    if (CovList[covnr[v]].cov==NULL) {\
      ERROR=ERRORNOTPROGRAMMED; goto ErrorHandling;}\
    i += CovList[covnr[v]].kappas;\
    }\
  if (i!=*np) {ERROR=ERRORPARAMNUMBER; goto ErrorHandling;}

// note: number of parameters is not checked whether it is correct!!
// done by R function `PrepareModel'
void CovarianceNatSc(Real *x, int *lx, int *covnr, Real *p, int *np, 
		     int *logicaldim, /* timespacedim ! */
		     int *xdim,
		     int *ncov, int *anisotropy, int *op,
		     Real *result, int *naturalscaling)
{
  // note: column of x is point !!
  int error, i,pAdd, v;
  param_type param;

  error = NOERROR;
  FIRST_CHECK_COV_VARIO(error);  

  for (v=0; v<*ncov; v++) {
    if (((error=(CheckAndRescale(covnr[v], *naturalscaling, p, *logicaldim,
				*anisotropy, param[v])))!=NOERROR) ||
	CovList[covnr[v]].variogram) 
      goto ErrorHandling;    
    // increase of address !
    p += pAdd+CovList[covnr[v]].kappas; // pointer addition!
  }
  
  for (i=0; i<*lx; i++, x += *xdim) {
    result[i] = CovFct(x, *xdim, covnr, op, param, *ncov, *anisotropy);
  }
  return;
  
 ErrorHandling:
  if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,error);
  for (i=0;i<*lx;i++) { result[i] = RF_NAN; }
  return;	
}

void Covariance(Real *x,int *lx, int *covnr, Real *p, int *np, 
		int *logicaldim, /* timespacedim ! */
		int *xdim,
		int *ncov, int *anisotropy, int *op,
		Real *result)
{
  CovarianceNatSc(x,lx,covnr,p,np,logicaldim,xdim,
		  ncov, anisotropy, op,result, &GENERAL_NATURALSCALING);
}

static param_type Unchecked_param;
static int Unchecked_covnr[MAXCOV], Unchecked_DIM, Unchecked_op[MAXDIM],
  Unchecked_ncov;
static bool Unchecked_anisotropy;

void InitUncheckedCovFct(int *covnr,Real *p,int *np, 
			 int *logicaldim, /* timespacedim ! */
			 int *xdim,
			 int *ncov, int *anisotropy, int *op,
			 int *naturalscaling, int *error){
  int i, j,  pAdd, v;
  bool variogram;
 
  FIRST_CHECK_COV_VARIO(*error);
  for (v=0; v<*ncov; v++) {
    if ((*error=(CheckAndRescale(covnr[v], *naturalscaling, p, *logicaldim,
				 *anisotropy, Unchecked_param[v]))) || 
	CovList[covnr[v]].variogram) 
      goto ErrorHandling;
    p += pAdd+CovList[covnr[v]].kappas;
  }
  memcpy(Unchecked_covnr, covnr, sizeof(int) * *ncov);
  memcpy(Unchecked_op, op, sizeof(int) * *ncov);
  Unchecked_ncov = *ncov;
  Unchecked_anisotropy = (bool) *anisotropy;
  Unchecked_DIM = *xdim;
  return;
  
 ErrorHandling:
  if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,*error);
  return;	
}

void UncheckedCovFct(Real *x, int *lx, Real *result)
{
  int i;
  /*
     WARNING:  make sure that CovList is initialised,
                              covnr is correct,
                              p is fine according to covnr (length(p), value)
                              PracticalRange is included in p[SCALE] if 
			      necessary
  */
  for (i=0; i<*lx; i++, x += Unchecked_DIM) {
    result[i] = CovFct(x, Unchecked_DIM, Unchecked_covnr, 
		       Unchecked_op, Unchecked_param, 
		       Unchecked_ncov, Unchecked_anisotropy);
  }
}

void VariogramNatSc(Real *x,int *lx,int *covnr,Real *p,int *np,
		    int *logicaldim, /* timespacedim ! */
		    int *xdim,
		    int *ncov, int *anisotropy, int *op,
		    Real *result, int *naturalscaling)
{
  // note: column of x is point !!
  int error=NOERROR, i, j, pAdd, v;
  param_type param;

  FIRST_CHECK_COV_VARIO(error);

  for (v=0; v<*ncov; v++) {
    if (error=CheckAndRescale(covnr[v], *naturalscaling, p, *logicaldim, 
			      *anisotropy,
			      param[v])) goto ErrorHandling;
    if ( (v>0) && (op[v-1]) && 
	 (CovList[covnr[v]].variogram || CovList[covnr[v-1]].variogram)) {
      // only covariance functions may be multiplied!
      error=ERRORNOMULTIPLICATION; goto ErrorHandling;
    } 
    p += pAdd+CovList[covnr[v]].kappas;
  }
  for (i=0; i<*lx; i++, x += *xdim) {
    result[i] = Vario(x, *xdim, covnr, op, param, *ncov, *anisotropy);
  }
  return;
   
 ErrorHandling:
  if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,error);
  for (i=0;i<*lx;i++) { result[i] = RF_NAN; }
  return;	 
}

void Variogram(Real *x,int *lx,int *covnr,Real *p,int *np, 
	       int *logicaldim, /* timespacedim ! */
	       int *xdim,
	       int *ncov, int *anisotropy, int *op,
	       Real *result)
{
  VariogramNatSc(x,lx,covnr,p,np,logicaldim,xdim,ncov,anisotropy,op,
		 result,&GENERAL_NATURALSCALING);
}
  
void CovarianceMatrixNatSc(Real *dist,int *lx,int *covnr,Real *p,int *np, 
			   int *logicaldim, /* timespacedim ! */
			   int *xdim,
			   int *ncov, int *anisotropy, int *op,
			   Real *result, int *naturalscaling)
{
  // note: column of x is point !!
  int error=NOERROR, i, ii, j, pAdd, zaehler, v, endfor, lxP1;
  long lxq, ve, ho;
  bool variogram;
  param_type param;
  Real var;
 
 
  FIRST_CHECK_COV_VARIO(error);  
  for (v=0; v<*ncov; v++) {
    if ( ((error=(CheckAndRescale(covnr[v], *naturalscaling, p, *logicaldim, 
				 *anisotropy, param[v]))) != NOERROR) || 
	 CovList[covnr[v]].variogram) 
      goto ErrorHandling;
    p += pAdd + CovList[covnr[v]].kappas;  // pointer increment!
  }
   
  var = CovFct(ZERO, *xdim, covnr, op, param, *ncov, *anisotropy);
  lxP1 = *lx + 1;
  for (ii=*lx, i=0; ii>0; i+=lxP1, ii--) {
    result[i] = var;
    endfor = i + ii;
    for (ve = i + 1, ho = i + *lx; ve < endfor; ve++, ho += *lx, dist += *xdim){
      result[ve] = result[ho] = 
	CovFct(dist, *xdim, covnr, op, param, *ncov, *anisotropy);
    }
  }
  return;
  
 ErrorHandling:
  for(i=0; i<lxq; i++) result[i] = RF_NAN;
  if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,error); 
  return;	
}

void CovarianceMatrix(Real *dist, int *lx,int *covnr, Real *p, int *np,
		      int *logicaldim, /* timespacedim ! */
		      int *xdim,
		      int *ncov, int *anisotropy, int *op,
		      Real *result)
{
  CovarianceMatrixNatSc(dist, lx, covnr, p, np, logicaldim, xdim, ncov, 
			anisotropy, op,	result, &GENERAL_NATURALSCALING);
}


void VariogramMatrix(Real *dist,int *lx,int *covnr,Real *p,int *np, Real *result)
{ 
  // not programmed yet
  assert(false);
}

void insert_method(int* last_incompatible, key_type *key,
			 SimulationType m) {
//  printf("insert method\n");
  int idx;
  assert(key->S[key->n_unimeth]==NULL);
  assert(key->destruct[key->n_unimeth]==NULL);
  if (do_incompatiblemethod[m]!=NULL) {
    // first the non-compatible (i.e. which do not allow direct 
    // adding to the random field), then those which allow for 
    // direct adding
    // --if attention wouldn't be payed to the ordering of the methods,
    // the simulation algorithm had to create an intermediate field 
    // more frequently -- sometimes the fields are pretty large,
    // and then this procedure is quite helpful
    (*last_incompatible)++;

    key->unimeth_used[key->n_unimeth]=key->unimeth_used[*last_incompatible];
    key->unimeth[key->n_unimeth]=key->unimeth[*last_incompatible];
    key->S[key->n_unimeth]=key->S[*last_incompatible]; // necessary
    // only if traditional and further methods are tried due to
    // partial failure
    key->destruct[key->n_unimeth]=key->destruct[*last_incompatible];

    key->S[*last_incompatible] = NULL;
    key->destruct[*last_incompatible] = NULL;
    idx = *last_incompatible;
  } else idx=key->n_unimeth;
  key->unimeth[idx] = m;
  key->unimeth_used[idx] = false;
  key->n_unimeth++;
//  printf("end insert method\n");
}

void delete_method(int *M, int *last_incompatible, key_type *key) {
  // note that key->n_unimeth had been decreased before this function
  // is called
  int idx;

//  printf("delete method\n");
  if (key->destruct[*M]!=NULL) key->destruct[*M](&(key->S[*M]));
  key->destruct[*M]=NULL;
  if (do_incompatiblemethod[key->unimeth[*M]]!=NULL) {
    assert(*last_incompatible>=0);
    assert(*M<=*last_incompatible);
    key->unimeth[*M]=key->unimeth[*last_incompatible];
    key->unimeth_used[*M] = key->unimeth_used[*last_incompatible];
    key->S[*M] = key->S[*last_incompatible];
    key->S[*last_incompatible]=NULL;
    key->destruct[*M]=key->destruct[*last_incompatible];
    key->destruct[*last_incompatible]=NULL;
    idx = *last_incompatible;
    (*last_incompatible)--;
  } else idx = *M;
//  printf("idx=%d %d %d %d\n", idx, *M ,*last_incompatible, key->n_unimeth);
  key->n_unimeth--;
  if (idx==key->n_unimeth){ // deleting of the very last entry
    assert(*M==key->n_unimeth); 
  } else {
    assert(idx>=0 && idx<key->n_unimeth);
    key->unimeth[idx] = key->unimeth[key->n_unimeth];
    key->unimeth_used[idx] = key->unimeth_used[key->n_unimeth];
    key->S[idx] = key->S[key->n_unimeth];
    key->S[key->n_unimeth]=NULL;
    key->destruct[idx] = key->destruct[key->n_unimeth];
    key->destruct[key->n_unimeth]=NULL;
  }
  key->unimeth_used[key->n_unimeth] = false;
  (*M)--;
//  printf("end delete method\n");
}

void InitSimulateRF(Real *x, Real *T, 
		    int *spatialdim, /* spacial dim only ! */
		    int *lx, int *grid, 
		    int *Time,
	            int *covnr, Real *ParamList, int *nParam,
		    Real *mean,
		    int *ncov, int *anisotropy, int *op,
		    int *method, 
		    int *distr, /* still unused */
		    int *keyNr , int *error)
{ // grid  : boolean
  // keyNr : label where intermediate results are to be stored;
  //         can be chosen freely between 0 and MAXKEYS-1

  // NOTE: if grid and Time then the Time component is treated as if
  //       it was an additional space component
  //       On the other hand, SPACEISOTROPIC covariance function will
  //       always assume that the second calling component is a time component!

  bool user_defined,finished;
  int i,act_number, last_incompatible, extndd_nr, d, endfor, 
    pAdd, v, m, M;
  unsigned long totalBytes;
  cov_fct *cov;
  Real newscale, p[TOTAL_PARAM], *PL;
  bool method_used[(int) Nothing];
  // preference lists, distinguished by grid==true/false and dimension
  // lists must end with Nothing!
  SimulationType Merr;
  SimulationType pg1[]={CircEmbed, CircEmbedLocal, Direct, AdditiveMpp, Nothing}; 
  SimulationType pg2[]={CircEmbed, CircEmbedLocal, TBM2, SpectralTBM,  
			 AdditiveMpp, TBM3, Direct, Nothing};
  SimulationType pg3[]={CircEmbed, CircEmbedLocal, TBM3, AdditiveMpp, 
			Direct, Nothing};
  SimulationType png1[]={Direct, AdditiveMpp, Nothing};
  SimulationType png2[]={TBM2, SpectralTBM, AdditiveMpp, TBM3, Direct, 
			 Nothing};
  SimulationType png3[]={TBM3, AdditiveMpp, Direct, Nothing};
  SimulationType *preference_list,first_method[MAXCOV];


  key_type *key;
  if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
    *error=ERRORREGISTER; goto ErrorHandling;}
  key = &(KEY[*keyNr]);
  key->storing = GENERAL_STORING;
 
  // check parameters

  key->timespacedim = *spatialdim + (int) (bool) (*Time);
  if (currentNrCov==-1) InitModelList(); assert(CovList!=NULL); 
  if ((*spatialdim<1) || (key->timespacedim>MAXDIM) || (*Time>1) || (*Time<0)){
    *error=ERRORDIM; goto ErrorHandling;}    
  if (x==NULL) {*error=ERRORCOORDINATES;goto ErrorHandling;}
  for (i=0; i<*ncov; i++) 
    if ((covnr[i]<0) | (covnr[i]>=currentNrCov)) {
      *error=ERRORCOVNROUTOFRANGE; goto ErrorHandling; }
  if (*anisotropy) {
   pAdd = 1 + key->timespacedim * key->timespacedim;
   // not literally "total"param; there might be gaps in the sequence of the 
   // parameters (if not all seven KAPPA parameters are used)
   key->totalparam = ANISO + key->timespacedim * key->timespacedim; 
  }
  else {
    pAdd = 2;
    key->totalparam = SCALE + 2; // SCALE, INVSCALE 
  }
  
  // also checked by R function `PrepareModel'
  for (i=pAdd * *ncov, v=0; v<*ncov; v++) {
    i += CovList[covnr[v]].kappas;
  }
  if (i!=*nParam) {*error=ERRORPARAMNUMBER; goto ErrorHandling;}
  // bytes to be stored for each vector of coordinates:

  totalBytes =  sizeof(Real) * *lx * *spatialdim;
  if (key->active) {   
    assert(key->x[0]!=NULL);
    key->active = 
      (key->lx==*lx) && 
      (key->spatialdim==*spatialdim) && 
      (key->grid==(bool)*grid) && 
      (key->distribution==*distr) &&
      (key->ncov==*ncov) &&
      (!memcmp(key->covnr, covnr, sizeof(int) * *ncov)) &&
      (key->anisotropy==*anisotropy) &&
      (key->mean==*mean) &&
      ( (*Time == key->Time) &&
	(!(*Time) || ! memcmp(key->T, T, sizeof(Real) * 3))) &&
      (!memcmp(key->op, op, (*ncov>0) ? *ncov-1 : 0)) &&
      (!memcmp(key->x[0], x, totalBytes)) 
      ;       
  }

  if (!key->active) {
    // if not active, all the parameters have to be 
    //                        stored again
    // save all the parameters if key has not been active   
    key->lx=*lx;
    key->spatialdim = *spatialdim;
    key->grid=(bool) *grid;
    key->distribution=*distr; 
    key->ncov=*ncov;
    memcpy(key->covnr, covnr, sizeof(int) * *ncov);
    key->anisotropy = (bool) *anisotropy;    
    key->mean = *mean;
    // operators +, *
    memcpy(key->op, op, (*ncov>0) ? *ncov-1 : 0);
    // coordinates

    if (key->x[0]!=NULL) free(key->x[0]);
    if ((key->x[0]=(Real*) malloc(totalBytes))==NULL){
      *error=ERRORMEMORYALLOCATION;goto ErrorHandling;
    }
    memcpy(key->x[0], x, totalBytes);

    for (d=1; d<*spatialdim; d++) key->x[d]= &(key->x[0][d * *lx]);
    for (; d<MAXDIM; d++)  key->x[d]=NULL;
    
    //Time
    if (key->Time = *Time) {
      if (!key->anisotropy) { *error = ERRORTIMENOTANISO; goto ErrorHandling; }
      memcpy(key->T, T, sizeof(Real) * 3);
      if (key->grid) {
	key->x[key->spatialdim] = key->T;
      }
    }
  } else { // key->active
    // CHECK:
    for (d = key->timespacedim; d<MAXDIM; d++) {assert(key->x[d]==NULL);} 
  }
 
  // check the parameter lists -- note that naturalscaling changes the values !
  for (v=0, PL=ParamList; v<key->ncov; v++){
    int j;
    if ((*error=(CheckAndRescale(covnr[v], GENERAL_NATURALSCALING, PL, 
				 key->timespacedim, *anisotropy, p))) != NOERROR)
      goto ErrorHandling; 
    for (j=KAPPA+CovList[covnr[v]].kappas; j<=LASTKAPPA; j++) p[j]=0.0;
    if (!(key->active &=
	  (key->covnr[v]==covnr[v]) &&
	  (!memcmp(key->param[v], p, sizeof(Real) * key->totalparam))
	  )) {
      key->covnr[v]=covnr[v];
      memcpy(key->param[v], p, sizeof(Real) * key->totalparam);
    }
   PL+=pAdd+CovList[covnr[v]].kappas;
  }

  // method list must be treated separately since method is integer coded
  // when passed to InitSimulate !
  if (user_defined = (method[0]>=0)) {//either all or none shld be user defined
    for (v=0; v<key->ncov; v++) {
      if (method[v] >= (int) Nothing) {
	*error=ERRORNOTDEFINED; goto ErrorHandling;}
      if ((key->method[v]!=(SimulationType) method[v])) {	
	key->active = false;
	key->method[v]=(SimulationType) method[v];
      }
    }
  } else {
    // TO DO: not consistent -- unclear which method to choose
    if (!key->active) {
      for (v=0; v<key->ncov; v++) {
	key->method[v]=CovList[covnr[v]].method(key->spatialdim, key->grid);
      }	
    }
  }   
  
  if (key->active) { // detected that all relevant parameters agree, 
    //              no initialization is necessary
    if (GENERAL_PRINTLEVEL>=2) ErrorMessage(Nothing, USEOLDINIT);
    *error=NOERROR; 
   return; /* ***** END ***** */ 
  }

  // Delete stored, intermediate result if there is any
  for (v=0; v<MAXCOV; v++) {
    if (key->destruct[v]!=NULL) {
      key->destruct[v](&(key->S[v]));
      key->destruct[v] = NULL;
    }
  }
  if (key->destructX!=NULL) {
    key->destructX(&(key->SX));
    key->destructX = NULL;
  }
  
  // Check
  if (!key->anisotropy)
    for (v=0; v<*ncov; v++) {
      assert(fabs(key->param[v][SCALE] * key->param[v][INVSCALE]-1.0) < EPSILON);
   }

  // minor settings, usually overtaken from previous run, if all the 
  // other parameters are identical
  
  // no further choice of TBM_METHOD !!!!!!!!!!!!! TO BE CHANGED IN FUTURE
  // if (TBM_METHOD==Nothing) {key->tbm_method=CovList[*covnr].tbm_method;}
  //else 
  { key->tbm_method=TBM_METHOD;}
 
  if (key->grid) { 
    if ((*lx!=3) || 
	(InternalGetGridSize(key->x,&key->timespacedim,key->length))) {
      *error=ERRORCOORDINATES; goto ErrorHandling;
    }
    for (d=key->timespacedim; d<MAXDIM; d++) key->length[d] = 1;
    for (key->spatialtotalpoints=1, d=0; d<key->spatialdim; d++) {
      key->spatialtotalpoints *= key->length[d];
    }
    key->totalpoints = key->spatialtotalpoints;
    if (key->Time) key->totalpoints *= key->length[key->spatialdim];

    // preference list
    switch (key->timespacedim) {
    case 1 : preference_list = pg1; break;
    case 2 : preference_list = pg2; break;
    case 3: case 4 : preference_list = pg3; break;
    default : assert(false);
    }
  } else { // not grid
    key->totalpoints = key->spatialtotalpoints = key->length[0]= *lx;
   if (key->Time) {
      int Tdim; Real *Tx[MAXDIM];
      Tdim = 1;
      Tx[0] = key->T;
      if (InternalGetGridSize(Tx,&Tdim,&(key->length[key->spatialdim]))) {
	*error=ERRORCOORDINATES; goto ErrorHandling;
      } 
      key->totalpoints *= key->length[key->spatialdim];
    } 
    // just to say that considering these values does not make sense
    for (d=1; d<key->spatialdim; d++) key->length[d]=0; 
    // correct value for higher dimensions (never used)
    for (d=key->timespacedim; d<MAXDIM; d++) key->length[d] = 1;
    switch (key->timespacedim) {
    case 1 : preference_list = png1; break;
    case 2 : preference_list = png2; break;
    case 3 : case 4 : preference_list = png3; break;
    default : printkey(key); assert(false);
    }
  }

  // simple (traditional in the sence of version 1.0 of RandomFields) 
  // definition of the covariance function?
  key->traditional = (!key->anisotropy) && ((key->ncov==1) || 
    ((key->ncov==2) && ( (key->method[0]!=Nugget) &&
		     key->method[1]==Nugget) ));

  // check (and set) method in case of multiplicative covariance function
  // here, the internal first prefence might be overwritten
  endfor = key->ncov - 1;
  if (user_defined) {
    for (v=0; v<endfor; v++) {  
      if (op[v]) { // *
	if (method[v]!=method[v+1]) {*error=ERRORMETHODMIX; goto ErrorHandling;}
      }
    } 
  } else {
    for (v=0; v<endfor; v++) {  
      if (op[v]) { // *
	key->method[v] = key->method[v+1] = (key->grid) ? CircEmbed : TBM3; 
      }
    }
  }

  for (v=0; v<*ncov; v++) first_method[v] = key->method[v]; 
  // a first prefered method is tried, before the preference_list
  // is consulted; first_method ist storred to avoid that 
  // it is tried a second time when consulting the preference_list

  finished = false; 
  assert(*error==NOERROR); 
  extndd_nr = 0; // used if !traditional; has same role as act_number 
  //               in the traditional case

  for (v=0; v<*ncov; v++) key->left[v]=true;
  // indicates which covariance function are left to be initialized

  key->n_unimeth = 0;
  act_number = -1;
  last_incompatible = -1;

  while (!finished) {
    int err_occurred;

// printf("\nYY %d %d %d %d %d\n", key->n_unimeth, M, last_incompatible, user_defined, key->traditional); for (v=0; v<MAXCOV; v++) printf(" %d ", key->destruct[v]);

    *error = err_occurred = NOERROR;// necessary in case ncov=0, since then 
    //                          the following loop is not entered


    // create list of methods used
    for (i=0; i<(int) Nothing; i++) { method_used[i] = false; }
    for (v=0; v<key->ncov; v++) {
      if (key->left[v]) { 
      method_used[key->method[v]] = true;  
      // if user_defined take over the user_defined one
      // if !user_defined take specifically preferred one (for 
      //                     each covariance function)
      }
    }

    for (m=CircEmbed; m<Nothing; m++) 
      if (method_used[m])
	insert_method(&last_incompatible, key, (SimulationType) m); 
        // increments key->n_unimeth

    assert((key->destructX==NULL) && (key->SX==NULL));
    for (M=0; M<key->n_unimeth; M++) if (!key->unimeth_used[M]) {
      bool left[MAXCOV];
      key->unimeth_used[M] = true; // optimistic, corrected by delete_method,
      //                              if method fails
      for (v=0; v<key->ncov; v++) left[v] = key->left[v]; // sicherung,
      // da init_method ueberschreibt und *danach* entscheidet, ob
      // die methode funktioniert -- ueber key->method ist es zu schwierig
      // da eine methode mehrmals verwendet werden kann
      if (GENERAL_PRINTLEVEL>=4)
	PRINTF("Init %d %s %d %s\n",M,
	       METHODNAMES[key->unimeth[M]], key->n_unimeth, "--");

// printf("\nXX %d %d %d %d %d\n", key->n_unimeth, M, last_incompatible, user_defined, key->traditional); for (v=0; v<MAXCOV; v++) printf(" %d ", key->destruct[v])

      assert((key->destruct[M]==NULL));
      assert((key->S[M]==NULL));
      if (key->unimeth[M]>Special) {
	err_occurred=ERRORMETHODNOTALLOWED; goto ErrorHandling;}
      err_occurred = init_method[key->unimeth[M]](key, M);
      if (GENERAL_PRINTLEVEL>=5) {
	PRINTF("return code %d :", err_occurred);
	ErrorMessage(Nothing, err_occurred);
      }
      if (err_occurred) {
	if (err_occurred==NOERROR_REPEAT) {
	  insert_method(&last_incompatible, key, key->unimeth[M]);
	  err_occurred = NOERROR;
	  // incompatible here means that partial simulation results
	  // are not directly added to the final result by the algorithm
	  // E.g. TBM is incompatible since at the final a division is performed
	} else if (err_occurred==NOERROR_ENDOFLIST) {
	  delete_method(&M, &last_incompatible, key);
	  err_occurred = NOERROR;
	} else { 
	  Merr = key->unimeth[M];
	  *error = err_occurred; // if not traditional, error may occur anywhere
	  //                     and final *error might be 0

// printf("\n%d %d %d %d %d\n",key->n_unimeth,  M, last_incompatible, user_defined, key->traditional); for (v=0; v<MAXCOV; v++) printf(" %d ", key->destruct[v]);

	  if (user_defined || key->traditional) break;
	  delete_method(&M, &last_incompatible, key);

// printf("\nZZ %d %d %d %d %d\n",key->n_unimeth,  M, last_incompatible, user_defined, key->traditional); for (v=0; v<MAXCOV; v++) printf(" %d ",key->destruct[v]);

	  for (v=0; v<key->ncov; v++) key->left[v] = left[v];
	}
      }
    }
    key->compatible = last_incompatible <= 0; // if at most one incompatible
    // method, then no intermediate results must be stored

    if (!(finished=user_defined || *error==NOERROR)) {
      if (key->traditional) {
	// tabula rasa, since the (true) covariance or the nugget may
	// have failed; just to be sure...
	for (m=0; m<key->n_unimeth; m++) {
	  if (key->destruct[m]!=NULL) key->destruct[m](&(key->S[m]));
	  key->destruct[m]=NULL;
	}
	if (key->destructX!=NULL) key->destructX(&(key->SX));
	key->destructX=NULL;
	key->n_unimeth = 0;

	for (v=0; v<*ncov; v++) key->left[v]=true;
	// indicates which covariance function are left to be initialized

	// there is a first_method given by the programmer when defining
        // a cov model; this method is tried first; afterwards the 
        // standard list of methods is checked; however the first_method
	// should not be reconsidered!
	if (GENERAL_PRINTLEVEL>1) ErrorMessage(key->method[0], *error);
	act_number++; // within the preference list, act_number is initially -1
	if (preference_list[act_number]==first_method[0]) {
	  // the new(next or first) method out of the standard list might be 
	  // the first_method
	  act_number++;
	}
	key->method[0]=preference_list[act_number];
	// end of the list?
	finished=finished || (preference_list[act_number]==Nothing);
      } else { // not traditional
	assert(key->S[M]==NULL);
	act_number++;
	for (v=0; v<key->ncov; v++) {
	  if (key->left[v]) key->method[v] = preference_list[act_number];
	}
	finished=finished || (preference_list[act_number]==Nothing);
      }
    } //if !finished || err_occurred
  } // while !finished

  key->active = !*error;
  if ((*error && (GENERAL_PRINTLEVEL>0)) ||
      (!user_defined && GENERAL_PRINTLEVEL>1)) { 
    if (!user_defined && !key->traditional && *error)
      PRINTF("algorithm failed. Last error has been: ");
    if (*error && GENERAL_PRINTLEVEL>0) {
      ErrorMessage(Merr, *error);
      PRINTF("\n");
    }
  }
  if (!key->active) goto ErrorHandling;

  if (GENERAL_PRINTLEVEL>=4) PRINTF("Successfully initialised.\n");
  return;
  
  ErrorHandling:
  DeleteKey(keyNr);  
}

  
void DoPoissonRF(int *keyNr, Real *res, int *error)
{
  // not programmed yet
  assert(false);
}

void DoSimulateRF(int *keyNr, Real *res, int *error)
{
  Real  *part_result;
  long nx, i, endfor, m;
  key_type *key;

  part_result = NULL;
  key = &(KEY[*keyNr]);
  *error=NOERROR; 

  if ((*keyNr<0) || (*keyNr>=MAXKEYS)) {
    *error=ERRORKEYNROUTOFRANGE; goto ErrorHandling;
  }
  if (!key->active) {*error=ERRORNOTINITIALIZED; goto ErrorHandling;}

  if (key->n_unimeth==0) {
    for(i=0; i<key->totalpoints; i++) res[i] = key->mean;
    return;
  } 
  
  GetRNGstate();
  if (!key->compatible) {
    if ((part_result=(Real *) malloc(sizeof(Real) * key->totalpoints))
	==NULL) {*error=ERRORMEMORYALLOCATION; goto ErrorHandling;}
  }
  for(m=0; m<key->n_unimeth; m++) {
    if (GENERAL_PRINTLEVEL>=7)
      PRINTF("DoSimu %d %s\n",m, METHODNAMES[key->unimeth[m]]);
    if (do_compatiblemethod[key->unimeth[m]]!=NULL) {
      if (GENERAL_PRINTLEVEL>=7) PRINTF("compatible: add=%d\n",m!=0);
      do_compatiblemethod[key->unimeth[m]] (key, m!=0, m, res ); 
    } else {
      if (m==0) {
	if (GENERAL_PRINTLEVEL>=7) PRINTF("incomp\n");
	do_incompatiblemethod[key->unimeth[m]](key, m, res );
      }
      else {
	if (GENERAL_PRINTLEVEL>=7) PRINTF("extra %d\n",key->compatible);
	do_incompatiblemethod[key->unimeth[m]](key, m, part_result
					       );
	for (i=0; i<key->totalpoints; i++) {
	  res[i] += part_result[i];
	}
      }
    }
  } // for m
  
  if (key->mean!=0.0) {
    for(i=0; i<key->totalpoints; i++) res[i] += key->mean;
  }  
  
  if (part_result!=NULL) free(part_result);
  if (! (key->active = GENERAL_STORING && key->storing)) DeleteKey(keyNr);
  PutRNGstate();
  return;
  
 ErrorHandling: 
  PutRNGstate();
  if (GENERAL_PRINTLEVEL>0) ErrorMessage(Nothing,*error);
  if (part_result!=NULL) free(part_result);
  key->active = false;
  DeleteKey(keyNr);
  return;
}

void SimulateRF(Real *x, Real *T, int *dim,int *lx, int *grid, 
		int *Time, 
		int *covnr, Real *ParamList, int *nParam,
		Real *mean,
		int *ncov, int *anisotropy, int *op,
		int *method, int *distr,
		int *keyNr,
		Real *res, int *error)
{
  // grid   : boolean; if true then x,y,z are supposed to have three values: 
  //                      start, end, step
  //                   if false then (x,y,z) is a list of points where the 
  //                      value of the random field should be simulated
  // x,y,z  : see grid; if z==NULL it is assumed that the field is at most 
  //                    2 dimensional
  //                    if y==z==NULL, the random field is 1 dimensional
  // covnr  : nr of covariance function by getCovFct,
  // ParamList: vector with 5 to 8 elements
  //          see RFsimu.h, the definition of MEAN..KAPPAII, for the meaning 
  //          of the components
  // method : if <0 then programme chooses automatically for the method
  //          if ==-1 it starts with the preferred ("first") method for the 
  //                  specific covariance fct
  //            else starts search with previous successful method (if any)
  //          otherwise see list `SimulationType' in RFsimu.h
  // keyNR  : label where intermediate results should be stored
  // res : array/grid of simulation results

  InitSimulateRF(x, T, dim, lx, grid, Time, covnr, ParamList, nParam,
                 mean, ncov, anisotropy, op, method, distr, keyNr, error);
  if (*error==NOERROR) {
    DoSimulateRF(keyNr,res,error);
  } else {    
    if (GENERAL_PRINTLEVEL>0) PRINTF(" ** All methods have failed. **\n");
  } 
}