estfuncR.cpp

Go to the documentation of this file.

/*!
  \file   estfuncR.cpp
  \author Klaus Holst
  \date   Jan 2005
  
  \brief  The interface to R.
    
*/

#include <map>
#include <string>
#include <cstdlib>
#include <cmath>
#include "linalg.h"
#include "linrdist.h"
#include "qscr.h"
#include "models.h"
#include "svd.h"
#include "debug.h"
#include "util.h"
// #include "siman.h"


// Simulation and global options
bool extraoutput = false;
long N = 100;
double d = 1.0;
long q = 5;
long P = 100;
std::string prefix = "out";
long s = 12345;
long l=10;
//matrix par(3);

// Optimization options:
double start_tt = 0.05;
double stop_tt = 0.0;
unsigned aiter = 200;
unsigned nosteps = 3;

PREESTMODEL *curModel;
PREpointerFn *curF;

double bestvalue;
  
double call_estfunc(matrix const& x) {
  double val = callMemberFunction(*curModel,*curF)(x);
  return val;
}
  

double cubic(double x) {
  return x*x*x;
}

double cubicabs(double x) {
  return x*x*x;
}

double square(double x) {
  return x*x;
}
 
double identity(double x) {
  return x;
}
  
double mysqrt(double x) {
  return sqrt(x);
}
  
double mysqrtabs(double x) {
  return sqrt(fabs(x));
}

double myabs(double x) {
  return fabs(x);
}

double myexp(double x) {
  return exp(x);
}

double mylogabs(double x) {
  return log(fabs(x));
}

double myx4(double x) {
  double x2 = x*x;
  return x2*x2;
}


double estfunc(double *par) {
  unsigned n = (*curModel).getdim();
    matrix x(n);
    for (unsigned i=0; i<n; i++)
      x.set(par[i],i);
    // Evaluating Estimating Function in x
    qinit();    
    double val = call_estfunc(x);
    qend();    
    return(val);
}



// ************************* Export to R ****************************




std::map<std::string, int> modelslist; 

extern "C" {


  void estfunceval(char **model, double *indata, int *dimdata, double *par, int *dimpar,
               int *NN, double *dd, int *PP,
                   int *qq, int *ll, int *ss, int *innerfunc, int *diminnerfunc, int *calcAstar, 
                   int *debug,
                   double *res) {

    PREpointerFn f_sumH, f_estfunc;

    f_sumH = &PREESTMODEL::twonorm_sumH;
    f_estfunc = &PREESTMODEL::twonorm_estfunc;


    N = *NN;
    d = *dd;
    q = *qq;
    P = *PP;
    s = *ss;
    l = *ll;

    // inner functions:
    // 0: identity x
    // 1: |x|
    // 2: x^2    
    // 3: x^3
    // 4: x^4
    // 5: sqrt(abs())
    // 6: exp
    // 7: log(abs())
    // 8: |x|^3
    safevector <double (*)(double)> myfunctions;    
    for (int i=0; i<*diminnerfunc; i++) {
      switch(innerfunc[i]) {
      case 0:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = x" << std::endl;
        }
        myfunctions.push_back(&identity);
        break;
      case 1:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = |x|" << std::endl;
        }
        myfunctions.push_back(&myabs);
        break;
      case 2:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = x^2" << std::endl;
        }
        myfunctions.push_back(&square);
        break;
      case 3:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = x^3" << std::endl;
        }
        myfunctions.push_back(&cubic);
        break;
      case 4:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = x^4" << std::endl;
        }
        myfunctions.push_back(&myx4);
        break;
      case 5:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = sqrt(|x|)" << std::endl;
        }
        myfunctions.push_back(&mysqrtabs);
        break;
      case 6:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = exp(x)" << std::endl;
        }

        myfunctions.push_back(&myexp);
        break;
      case 7:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = log(|x|)" << std::endl;
        }

        myfunctions.push_back(&mylogabs);
        break;
      case 8:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = |x|^3" << std::endl;
        }
        myfunctions.push_back(&cubicabs);
        break;
      default:
        if (*debug) {
          std::cerr << "f" << i+1 << "(x) = |x|" << std::endl;
        }
        myfunctions.push_back(&myabs);
        break;
      }
    }
    

    enum Models { Error = 0, SVCIR_MODEL, SVEOU_MODEL };
    modelslist["svcir"] = SVCIR_MODEL; 
    modelslist["sveou"] = SVEOU_MODEL;     

    SVCIR mSVCIR(myfunctions,d,P,N,N,q,q+1,l,s);      
    SVEOU mSVEOU(myfunctions,d,P,N,N,q,q+1,l,s);  

    switch (modelslist[std::string(*model)]) { 
    case SVCIR_MODEL: 
      if (*debug) {
        std::cerr << "Using SVCIR-model" << std::endl;
      }
      curModel = &mSVCIR;
      break; 
    case SVEOU_MODEL:     
      if (*debug) {
        std::cerr << "Using SVEOU-model" << std::endl;
      }
      curModel = &mSVEOU;
      break;      
    default:
      std::cerr << std::endl << "You specified an unknown model" << std::endl;
      break;
    }   

    curF = &f_estfunc;

    if (modelslist[std::string(*model)]!=Error) {

      matrix datamatrix(*dimdata);
      for (unsigned i=0; i<*dimdata; i++) 
        datamatrix.set(indata[i],i);
      (*curModel).setdata(datamatrix);

      matrix p(*dimpar);
      for (unsigned i=0; i<*dimpar; i++)
        p.set(par[i],i);     

      if ((*debug)!=0) {
        (*curModel).printdebug = 1;
      }


      if (*calcAstar) {
        curF = &f_estfunc;
        (*curModel).set_sumH_est(p.get(1)); // theta-estimate
        (*curModel).set_iter_Mbar(1); // Make sure Astar is calculated
      }
      else {
        curF = &f_sumH; 
      }
      double val = call_estfunc(p);

      *res = val;
    }
    
  }


  void simdata(char **model, double *par, int *dimpar, int *NN, double *dd, int*PP, int *ss,  double *data) {

    N = *NN;   
    d = *dd;
    P = *PP;
    s = *ss;

    matrix p(*dimpar);
    for (unsigned i=0; i<*dimpar; i++)
      p.set(par[i],i);     

    safevector <double (*)(double)> myfunctions;
    myfunctions.push_back(&square);


    enum Models { Error = 0, SVCIR_MODEL, SVEOU_MODEL };
    modelslist["svcir"] = SVCIR_MODEL; 
    modelslist["sveou"] = SVEOU_MODEL;     

    SVCIR mSVCIR(myfunctions,d,P,N,N,q,q+1,l,s);      
    SVEOU mSVEOU(myfunctions,d,P,N,N,q,q+1,l,s);

    switch (modelslist[std::string(*model)]) { 
    case SVCIR_MODEL: 
      curModel = &mSVCIR;
      break; 
    case SVEOU_MODEL:     
      curModel = &mSVEOU;
      break;      
    default:
      std::cerr << std::endl << "You specified an unknown model";
      break;
    }   

    if (modelslist[std::string(*model)]!=Error) {      
      (*curModel).setpar(p);
      matrix d = (*curModel).simulate(N);
      for (int i=0; i<N; i++)
        data[i] = d.get(i);
    }
  
  }


}

---