/*
 Authors 
 Martin Schlather, schlath@hsu-hh.de 

 Copyright (C) 2001 -- 2006 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.  
*/

/* TO DO
  * noch fast alles auf key->x programmiert statt s->x
  * time component: aenderungen siehe TBM !
*/

#include <math.h>
#include <assert.h>
#include "RFsimu.h"
#include "MPPFcts.h"

double MPP_RADIUS =0.0;
double MPP_APPROXZERO  =0.001;
double ADDMPP_REALISATIONS = 100.0; // logically an integer, implementation 
                                  //  allows for positive real ! 

void SetMPP(int *action, double *approxzero, double *realisations, 
	    double *radius)
{
  if (*action) {
    if (*approxzero<=0.0) {
      if (GENERAL_PRINTLEVEL>0) 
	PRINTF("\nWARNING! `approx.zero' not a positive number. Ignored.\n");
    } else  MPP_APPROXZERO = *approxzero;
   if (*realisations<=0.0) {
     if (GENERAL_PRINTLEVEL>0) 
       PRINTF("\nWARNING! `realisation' not a positive number. Ignored.\n");
   } else ADDMPP_REALISATIONS = *realisations;
   MPP_RADIUS = *radius;
   if ((MPP_RADIUS>0.0) && (GENERAL_PRINTLEVEL>1))
     PRINTF("\nWARNING! Positive `addradius' may result in simulation results that are difficult to detect.\n");
  } else {
    *approxzero = MPP_APPROXZERO;
    *realisations =  ADDMPP_REALISATIONS;
    *radius = MPP_RADIUS;
  }
}

void mpp_destruct(void **S) {
  // mpp_storage in RFsimu.h !!
  if (*S!=NULL) {
    free(*S);
    *S = NULL;
  }
}

void MPP_NULL(mpp_storage* s) {
  int i;
  s->integral = s->integralsq = s->integralpos = s->factor = s->maxheight = 
    s->effectiveRadius = s->effectivearea = s->addradius = RF_NAN;
  for (i=0; i<MAXDIM; i++) s->min[i] = s->length[i] = RF_NAN;
  for (i=0; i<6; i++) s->c[i] = RF_NAN;
  s->MppFct = NULL;
  s->dim = 0;  // otherwise error possible in InternalGetKeyInfo
}

int init_mpp(key_type * key, int m) {
  methodvalue_type *meth; 
  int error, d, i, v, reduceddim, actcov;
  double max[MAXDIM];  
  mpp_storage *s;
  covinfo_type *kc;

  error = NOERROR;
  meth = &(key->meth[m]);
  SET_DESTRUCT(mpp_destruct, m);
  assert(meth->S==NULL);
  if ((meth->S=malloc(sizeof(mpp_storage)))==0){
    error=ERRORMEMORYALLOCATION; goto ErrorHandling;
  }
  s = (mpp_storage*) meth->S;
  MPP_NULL(s);
  kc = NULL;
  for (actcov=v=0; v<key->ncov; v++) {
    ERRORMODELNUMBER = v;	
    kc = &(key->cov[v]);
    if ((kc->method==AdditiveMpp) && kc->left
        && (kc->param[VARIANCE]>0)) {
      cov_fct *cov;
      assert(kc->nr>=0 && kc->nr<currentNrCov);
      meth->covlist[actcov] = v;
      cov = &(CovList[kc->nr]);
      if (key->Time && !kc->genuine_last_dimension && 
	  (cov->type==SPACEISOTROPIC)) {// unclear whether this is too cautious
	  error= ERRORWITHOUTTIME; goto ErrorHandling;}
      if ((error=cov->check(kc->param, kc->reduceddim, kc->dim,
			    AdditiveMpp)) != NOERROR) {
	ERRORMODELNUMBER = v;	
	goto ErrorHandling;
      }
      if (cov->add_mpp_scl==NULL) {error=ERRORNOTDEFINED; goto ErrorHandling;}
      if (cov->type==ISOHYPERMODEL || cov->type==ANISOHYPERMODEL) {
	error=ERRORHYPERNOTALLOWED; goto ErrorHandling;}
      if ((v<key->ncov-1 &&  kc->op)) {
	// no multiplicative models are allowed
	  error=ERRORNOMULTIPLICATION; goto ErrorHandling;}
      if ((error=Transform2NoGrid(key, v)) != NOERROR) goto ErrorHandling;
      kc->left=false;
      actcov++; 
      // do not pile up more than 1 elementary covariance function.
      // creating max-stable fields may cause difficulties! 
      break;
    }
  }
  ERRORMODELNUMBER = -1;	
  if (actcov==0) {
    if (error==NOERROR) error=NOERROR_ENDOFLIST;
    goto ErrorHandling;
  } else assert(actcov==1);

  assert(kc != NULL);
  if (v>key->ncov) {
    error=ERRORFAILED; goto ErrorHandling;
  }
  meth->actcov = actcov;
  // kc = &(key->cov[0]); -- passt nicht waere v
  // determine the minimum area where random field values are to be generated
  reduceddim = kc->reduceddim;
  if (kc->simugrid) {
    for (d=0; d<key->timespacedim; d++) {
//      printf("%d, %d %d %f %f\n", d, kc->idx[d],
//	     XSTARTD[d], kc->xsimugr[XSTARTD[d]], s->integralpos);
      s->min[kc->idx[d]]   = kc->xsimugr[XSTARTD[d]];
      s->length[kc->idx[d]]= kc->xsimugr[XSTEPD[d]] * 
	(double) (key->length[d] - 1);
      // x[XSTARTD[d]] and x[XSTEPD[d]] are assumed to be precise, but
      // not x[XENDD[d]] 
    }
  } else {
    int ix;
      for (d=0; d<reduceddim; d++) {
      s->min[d]=RF_INF; 
      max[d]=RF_NEGINF;
    }
    v = 0;
    kc = &(key->cov[meth->covlist[v]]);
    for (ix=i=0; i<key->totalpoints; i++, ix+=reduceddim) {
      for (d=0; d<reduceddim; d++) {
	if (kc->x[ix+d] < s->min[d]) s->min[d] = kc->x[ix+d];
	if (kc->x[ix+d] > max[d]) max[d] = kc->x[ix+d];
      }
    }
    for (d=0; d<reduceddim; d++) 
      s->length[d] = max[d] - s->min[d];
  }
  v = 0;
  kc = &(key->cov[meth->covlist[v]]);
  s->addradius = MPP_RADIUS;  // must be set BEFORE the next command!!
  CovList[kc->nr].add_mpp_scl(s, reduceddim, kc->param);
  s->MppFct = CovList[kc->nr].add_mpp_rnd;
  if ((MPP_RADIUS>0.0) && (GENERAL_PRINTLEVEL>=2))
    PRINTF("Note: window has been enlarged by fixed value (%f)\n",s->addradius);
 
  if (key->anisotropy) return NOERROR_REPEAT;
  else return 0;

 ErrorHandling:
  return error; 
}


void do_addmpp(key_type *key, int m, double *res )
{ 
  methodvalue_type *meth; 
  covinfo_type *kc;
  int i, d, v, timespacedim;
  double lambda, //[MAXDIM], 
    min[MAXDIM], max[MAXDIM], 
    factor, average, poisson;
  long segment[MAXDIM+1];
  mpp_storage *s;
  mppmodel model;

  meth = &(key->meth[m]);
  assert(meth->actcov==1);
  {
    register long i;
    for (i=key->totalpoints - 1; i>=0; res[i--]=0.0);
  }

  s =  (mpp_storage*) meth->S;  
  v = 0;  // NOTE !
  
  switch (key->distribution) {
  case DISTR_GAUSS : 
    kc = &(key->cov[meth->covlist[v]]);
    lambda = ADDMPP_REALISATIONS * s->effectivearea;
    factor = sqrt(kc->param[VARIANCE] / 
		  (ADDMPP_REALISATIONS * s->integralsq));
    average = ADDMPP_REALISATIONS * s->integral;
    break;
  case DISTR_POISSON :
    // not done yet -- it is only a fragment
    // on the R level, it is excluded that the program runs 
    // into this part
    //
    // only actcov==1 is possible!!
    //
    // maybe nice to switch MPPFcts back to actcov=1
    kc = &(key->cov[meth->covlist[v]]);
    lambda = kc->param[VARIANCE] * s->effectivearea; 
    factor = average = 0.0 /* replace by true value !! */;  
    assert(false);
    if (key->mean!=0) {
      // do not compare mean against VARIANCE since
      // VARIANCEs may differ with v
      //      if (GENERAL_PRINTLEVEL>0)
      ErrorMessage(AdditiveMpp,ERRORVARMEAN); assert(false); // print it always
      
      assert(false); // for debugging -- delete this later on
      
      // finally, the value of MEAN is ignored, and that of VARIANCE
      // is used for simulating the random field
    }
    break;
  default : assert(false);
  }
  
  segment[0] = 1;  // only used if simugrid
  for (d=0; d<key->timespacedim; d++)  
      segment[d+1] = segment[d] * key->length[d];
  for (i=0; i<key->totalpoints; i++) res[i]=0.0;
  kc = &(key->cov[meth->covlist[v]]);
  timespacedim = kc->simugrid ? key->timespacedim : kc->reduceddim;
  for (poisson = rexp(1.0); poisson < lambda; poisson += rexp(1.0)) {   
    (s->MppFct)(s, min, max, &model);
    if (kc->simugrid) {
      int start[MAXDIM], end[MAXDIM], resindex, index[MAXDIM],
	  segmentdelta[MAXDIM], dimM1;
      double coord[MAXDIM], startcoord[MAXDIM];
	  
      resindex = 0;
      for (d=0; d<timespacedim; d++) {	 
        // determine rectangle of indices, where the mpp function
        // is different from zero
	if (kc->genuine_dim[d]) {
          if (min[kc->idx[d]]< kc->xsimugr[XSTARTD[d]]) {start[d]=0;}
	  else start[d] = (int) ((min[kc->idx[d]] - kc->xsimugr[XSTARTD[d]]) / 
				 kc->xsimugr[XSTEPD[d]]);
	  // "end[d] = 1 + *"  since inequalities are "* < end[d]" 
	  // not "* <= end[d]" !
	  end[d] = 1 + (int) ((max[kc->idx[d]] - kc->xsimugr[XSTARTD[d]]) /
			      kc->xsimugr[XSTEPD[d]]);
	  if (end[d] > key->length[d]) end[d] = key->length[d];	 
	} else {
	    start[d] = 0;
	    end[d] = key->length[d];
	}
	  
	// prepare coord starting value and the segment vectors for res
	index[d] = start[d];
	segmentdelta[d] = segment[d] * (end[d]-start[d]);
	resindex += segment[d] * start[d];
	coord[kc->idx[d]] = startcoord[d] = 
	    kc->xsimugr[XSTARTD[d]] + (double) start[d] * kc->xsimugr[XSTEPD[d]];
      }

      // add mpp function to res
      dimM1 = timespacedim -1;
      d = 0; // only needed for assert(d<dim)
      while (index[dimM1] < end[dimM1]) {
	assert(d<timespacedim); 
	res[resindex] += model(coord); 
	d=0;
	index[d]++;
	coord[kc->idx[d]] += kc->xsimugr[XSTEPD[d]];
	resindex++;
	while (d<dimM1 && index[d] >= end[d]) { 
	  // loop never entered if dim=1
	  index[d] = start[d];
	  coord[kc->idx[d]] = startcoord[d];
	  resindex -= segmentdelta[d];
	  d++; // if (d>=dim) break;
	  index[d]++;
	  coord[kc->idx[d]] += kc->xsimugr[XSTEPD[d]];
	  resindex += segment[d];
	}
      }
    } else { // !key->simugrid
      double y[MAXDIM];
      long j;
      // the following algorithm can greatly be improved !
      // but for ease, just the stupid algorithm
      for (j=i=0; i<key->totalpoints; i++) {
	for (d=0; d<timespacedim; d++, j++) y[d] = kc->x[j];
	res[i] += model(y); 
      }
    }
  }

  if (key->distribution==DISTR_GAUSS) {
    for (i=0; i<key->totalpoints;i++) {
      res[i] = factor * (res[i] - average);
    }
  }
}