/*
 Authors
 Martin Schlather, schlather@math.uni-mannheim.de

 Simulation of a random field by Cholesky or SVD decomposition

 Copyright (C) 2001 -- 2014 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 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 <math.h>
#include <stdio.h>  
#include <stdlib.h>
 
#include "RF.h"
#include "Covariance.h"
#include <R_ext/Lapack.h>
//#include <R_ext/Linpack.h>


int check_specificGauss(cov_model *cov) {
#define nsel 4
  cov_model
    *key = cov->key,
    *next=cov->sub[0];
  int err ; // taken[MAX DIM],
  //  direct_param *gp  = &(GLOBAL.direct); //
  
  ROLE_ASSERT(ROLE_GAUSS);
   
  if (cov->tsdim != cov->xdimprev || cov->tsdim != cov->xdimown) 
    return ERRORDIM;
  if (CovList[next->nr].Specific == MISMATCH)
    SERR1("specific method for '%s' not known", NAME(next));

  //    APMI(cov); //assert(false);
  if (key == NULL) {
#define SPEC_TYPES 4
    Types type[SPEC_TYPES] = {PosDefType, PosDefType, NegDefType, TrendType};
    int i,
      iso[SPEC_TYPES] = {SYMMETRIC, SYMMETRIC, SYMMETRIC, CARTESIAN_COORD},
      dom[SPEC_TYPES] = {XONLY, KERNEL, XONLY, XONLY};    
    
      // APMI(cov->calling);

    for (i=0; i<SPEC_TYPES; i++) {
      if ((err = CHECK(next, cov->tsdim,  cov->tsdim,
		       type[i], dom[i], iso[i],
		       SUBMODEL_DEP, ROLE_COV)) == NOERROR) break;
    }
    if (err != NOERROR) return err;
    if (next->pref[Specific] == PREF_NONE) return ERRORPREFNONE;
  } else {
    if ((err = CHECK(key, cov->tsdim,  cov->tsdim, ProcessType,
		     XONLY, cov->isoown,
		     SUBMODEL_DEP, ROLE_GAUSS)) != NOERROR) {
      return err;
    }

    //APMI(key);

  }
  cov_model *sub = cov->key == NULL ? next : key;
  setbackward(cov, sub);
  cov->vdim2[0] = sub->vdim2[0];
  cov->vdim2[1] = sub->vdim2[1];

  //PMI(cov);

  return NOERROR;
}


int struct_specificGauss(cov_model *cov, cov_model VARIABLE_IS_NOT_USED **newmodel) {
  cov_model 
    *next = cov->sub[0];
  location_type *loc = cov->prevloc;
   int err;
  //   PMI(cov); //assert(false);
  
  if (next->pref[Specific] == PREF_NONE) {
    return ERRORPREFNONE;
  }

  ROLE_ASSERT_GAUSS;

  if (cov->key != NULL) COV_DELETE(&(cov->key));
  if ((err = covcpy(&(cov->key), next)) != NOERROR) return err;
  if ((err = CHECK(cov->key, next->tsdim, next->xdimprev, next->typus, 
		   next->domprev, next->isoprev, next->vdim2,
		   next->role))!= NOERROR) {
    //PMI(cov->key); XERR(err);
    //printf("specific ok\n");
    // crash();
    return err;
  }

  assert(CovList[next->nr].Specific >= 0);


  cov->key->nr = CovList[cov->key->nr].Specific ;
  cov->key->role = ROLE_GAUSS;
  cov->key->typus = ProcessType;

  //PMI(cov->key);
  // APMI(cov->key);
  if ((err = STRUCT(cov->key, NULL)) != NOERROR) {
    return err;
  }

  //   APMI(cov->key);

  if ((err = CHECK(cov->key, loc->timespacedim, cov->xdimown, ProcessType,
		   XONLY, CARTESIAN_COORD, cov->vdim2, ROLE_GAUSS)) != NOERROR) {
    //PMI(cov->key); XERR(err);
    //printf("specific ok\n");
    // crash();
    return err;
  }

  //APMI(cov->key);
  //  if (next->nr == PLUS) AERR(err);
 
  return err;
}



int init_specificGauss(cov_model *cov, gen_storage *S) {
  cov_model *key = cov->key;
  int err;

  assert(key != NULL);

  if (cov->role == ROLE_COV) {
    return NOERROR;
  }

  ROLE_ASSERT_GAUSS;

  cov->method = Specific;
  if ((err = INIT(key, 0, S)) != NOERROR) return err;

  key->simu.active = true;
  cov->fieldreturn = true;
  cov->origrf = false;
  cov->rf = key->rf;

  return err;
}


void do_specificGauss(cov_model *cov, gen_storage *S) {  
  cov_model *key = cov->key;
  location_type *loc = Loc(cov);
  bool loggauss = GLOBAL.gauss.loggauss;
  double *res = cov->rf;

  assert(key != NULL);
  DO(key, S);
  if (loggauss) {
    long i, vdimtot = loc->totalpoints * cov->vdim2[0];
    for (i=0; i<vdimtot; i++) res[i] = exp(res[i]);
  }
}