estfunc.cpp

Go to the documentation of this file.

/*!
  \file   estfunc.cpp
  \author Klaus Holst
  \date   Jan 2005
  
  \brief  Demonstration of estimating function. Optimization via simulated annealing.
    
*/

#include <cmath>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <ctime>
#include <cstring>
#include <vector>
#include <map>
#include "linalg.h"       // various (matrix)functions
#include "models.h"
#include "util.h"
#include "debug.h"
#include "gnuplot.h"
#include "siman.h"
#include "qscr.h"


// Simulation and global options
bool extraoutput = false;
long N = 200;
double d = 1.0;
long q = 5;
long P = 20;
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 = 100;


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 square(double x) {
  return x*x;
}

double identity(double x) {
  return x;
}

double mysqrt(double x) {
  return sqrt(x);
}

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


void setsimplex(matrix& p) {
  p.newsize(4,3);
    
  p.set(1,0,0);
  p.set(1,0,1);
  p.set(1.5,0,2);

  p.set(1.5,1,0);
  p.set(1,1,1);
  p.set(1,1,2);

  p.set(1.5,2,0);
  p.set(1.5,2,1);
  p.set(1,2,2);

  p.set(1,3,0);
  p.set(1,3,1);
  p.set(1.5,3,2);
}




matrix anneal(PREESTMODEL *M, SimAnneal &sa, double epsilon, unsigned steps, double starttemp, matrix startpoint) {

//  (*M).set_iter_Mbar(1); // Calculate A(theta) next time estfunc is evaluated.
//  double estfunctrueval = call_estfunc(par);
//  string st = "estfunc(" + string(par.trans().row2string(0)) + ") = " + numStr(estfunctrueval);
//  qmvaddstr(15,55, st.c_str());

//  (*M).set_iter_Mbar(1); // Calculate A(theta) next time estfunc is evaluated.
//  matrix pt(3); pt.set(0.966078,0); pt.set(0.478924,1); pt.set(0.556969,2);
//  estfunctrueval = call_estfunc(pt);
//  st = "estfunc(" + string(pt.trans().row2string(0)) + ") = " + numStr(estfunctrueval);
//  qmvaddstr(16,55, st.c_str());


  (*M).set_iter_Mbar(1); // Calculate A(theta) next time estfunc is evaluated.
  sa.calcyb(); // Do that now.
  matrix bestpoint = startpoint;
  //  std::cerr << "bestpoint = "; bestpoint.stdprint();
  //  std::cerr << "bestvalue = " << sa.get_yb() << std::endl;
  
  bestvalue = sa.funceval(bestpoint);

  //  std::cerr << "bestvalue = " << bestvalue << std::endl;
  double curtt = starttemp;
  for (unsigned i=0; i<nosteps; i++) {
    (*M).set_iter_Mbar(0); // Do NOT calculate A(theta) the next couple of times estfunc is evaluated.
    sa.set_tt(curtt);
    sa.set_simplex( newsimplex(bestpoint) );   
    sa.optimize();
    (*M).set_iter_Mbar(1);    // Set iter_Mbar = true and call calcyb to update A(theta).
    qmvaddstr(25,7, "Calculating yb (and hence updating A(theta))...");
    sa.calcyb();
    if (sa.get_yb()<bestvalue) {
      bestvalue = sa.get_yb();
      bestpoint = sa.get_pb();
    }
    //    std::cerr << "bestvalue = " << bestvalue << std::endl;
    qmove(25,1); qclrtoeol();
    qmove(18,1); qclrtoeol();
    curtt -= epsilon;
  }
  
  return bestpoint;
}


void parseArguments(int argc, char* argv[]) {
  char syntax[] = 
    "Program demonstrating 'Prediction-Based Estimating Functions'\n"
    "Syntax:\n"
    "  diffusion <options>\n"
    "Options:\n"
    " -q[i]       (The lag-length in the prediction)\n"
    " -N[i]       (The number of simulations used)\n"
    " -d[f]       (The time-distance between the observations)\n"
    " -P[i]       (Points per Delta (time-distance))\n"
    " -f[s]       (Output-file = s)\n"
    " -l[i]       (Summation-start in the expresion of Mbar)\n"
    " -s[i]       (Seed)\n"
    " -i[i]       (Number of optimization-iterations)\n"
    " -v          (Ekstra output)\n"
    "\n";
  
  if((argc>1)&&( strcmp(argv[1],"-h")==0)) { 
    std::cout << syntax;
    exit(0); 
  }
  
  for(int i=1; i<argc; ++i) { 
    char* arg = argv[i];
    std::string str(arg);
    std::string strval(str.substr(2,str.length()));
    char option = arg[1];
    safevector<float> paramlist = getparm(strval);     
    switch (option)
      {
      case 'q': 
        q = static_cast <long> (paramlist[0]);
        break;      
      case 'N':
        N = static_cast <long> (paramlist[0]);
        break; 
      case 'd': 
        d = static_cast <double> (paramlist[0]);
        break;      
      case 'P':
        P = static_cast <long> (paramlist[0]);
        break; 
      case 'l': 
        l = static_cast <long> (paramlist[0]);
        break;      
      case 's':
        s = static_cast <long> (paramlist[0]);
        break; 
      case 'f':
        prefix = strval;
        break;
        //     case 's':
        //       inwidth = static_cast <unsigned> (paramlist[0]);
        //       inheight = static_cast <unsigned> (paramlist[1]);
      case 'v':
        extraoutput = true;
        break;
      default:      
        throw std::string("Unknown option -") + option + ".\n" + syntax;
      }
  }
}


matrix SimulateAndEstimate() { 

  std::string line = "";
  for (unsigned i=0; i<50; i++) 
    line += "-";
  line += "\n";


  qmvaddstr(3,3, "Prediction-Based Estimating Function in SV-Model.");

  PREpointerFn f_sumH, f_estfunc;

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


  PREpointerFn A, B;
  A = &PREESTMODEL::twonorm_sumH;
  B = &PREESTMODEL::twonorm_estfunc;

  
  matrix simdata, estimate(3); 

  std::string st;

  qmvaddstr(5,3,"Options");
  st = "N = " + numStr(N);  qmvaddstr(6,5, st.c_str());
  st = "d = " + numStr(d);  qmvaddstr(7,5, st.c_str());
  st = "q = " + numStr(q);  qmvaddstr(8,5, st.c_str());
  st = "P = " + numStr(P);  qmvaddstr(9,5, st.c_str());
  st = "l = " + numStr(l);  qmvaddstr(10,5, st.c_str());
  st = "s = " + numStr(s);  qmvaddstr(12,5, st.c_str());  
  qmvaddstr(13,3,"Optimization options:");
  st = "start_tt = " + numStr(start_tt);  qmvaddstr(14,5, st.c_str());
  st = "stop_tt = " + numStr(stop_tt);  qmvaddstr(15,5, st.c_str());
  st = "aiter = " + numStr(aiter);  qmvaddstr(16,5, st.c_str());
  

  safevector <double (*)(double)> myfunctions;
  //  myfunctions.push_back(&identity);
  myfunctions.push_back(&square);
  myfunctions.push_back(&cubic);
  myfunctions.push_back(&myabs);
  //  myfunctions.push_back(&mysqrt);



    //  SVCIR testmodel(myfunctions,d,P,N,200,q,q+1,l,s);
  SVEOU testmodel(myfunctions,d,P,N,N,q,q+1,l,s);
    
    //  SVCIR testmodel(myfunctions,d,P,N,q,q+1,l,s);
  curModel = &testmodel;




  
  par.set(10.0, 0); 
  par.set(5.0, 1); 
  par.set(5.0, 2);

  testmodel.setpar(par);
  simdata = testmodel.simulate(N-1);
  testmodel.setdata(simdata);



  //  par.set(12.94339, 0); 
  //  par.set(5.43643, 1); 
  //  par.set(5.34323, 2);

  //  curF = &f_estfunc;
  //  testmodel.set_sumH_est(5); // theta-estimate
  //  double estfunctrueval = call_estfunc(par);
  //  st = "estfunc(truepar) = " + numStr(estfunctrueval);
  //  std::cerr << st << estfunctrueval << std::endl;
  //  exit(0);


  matrix p;
  setsimplex(p);
  
  qmvaddstr(24,10, "Calculating initial estimate using twonorm_sumH (simplified est.func.)"); qclrtoeol();
 
  {
    curF = &f_sumH;
    
    Function F(call_estfunc, 3); // Define 3d-function

    SimAnneal sa(&F, p, 0.0, 100);

    sa.optimize();
    estimate = sa.get_pb(); // parameter (alpha,theta,sigma)
    bestvalue = sa.get_yb();
  } // F and sa destructed

  //  std::cerr << "estimate="; estimate.stdprint();
  //  std::cerr << "value=" << bestvalue << std::endl; 
  
  testmodel.set_sumH_est(estimate.get(1)); // theta-estimate
    
  
  //  testmodel.set_sumH_est(1000); // theta-estimate

  qmvaddstr(24,10, "Calculating estimate using twonorm_estfunc."); qclrtoeol();

  p = newsimplex(p.getrow(0).trans()); // Setting new simplex around current parameter-estimate.

  {
    curF = &f_estfunc;
    Function F(call_estfunc, 3); // Define 3d-function
    SimAnneal sa(&F, p, start_tt, aiter); // Create annealing-object 
    estimate = anneal(curModel, sa, start_tt/(double)nosteps, nosteps, start_tt, p.getrow(0).trans());  // parameter (alpha,theta,sigma)
    //    bestvalue = sa.get_yb();
  } // F and sa destructed
  
  
  return estimate;

}
 

int main(int argc, char* argv[]) {
  try {
    parseArguments(argc,argv);

    qinit();
    matrix estimate = SimulateAndEstimate();
    qend();

    std::cerr << std::endl << "Program terminated successfully." << std::endl;
    std::cerr << "Minimum value: " << bestvalue << std::endl; 
    std::cerr << "ESTIMATE: [";
    unsigned n=estimate.row();
    for (unsigned i=0; i<n-1; i++) {
      std::cerr << estimate.get(i) << ", ";
    }
    std::cerr << estimate.get(n-1) << "]" << std::endl;

    std::cerr << std::endl;
     
    
    return EXIT_SUCCESS;
  } catch(std::string const& err) {
    std::cerr << "exn: " << err << std::endl;
    return EXIT_FAILURE;
  } 
}

---