///////////////////////////////////////////////////////////////////////
//functorEx.cpp
#include <iostream>
#include <cmath>

using std::cout; using std::endl;

struct Functor{
  int nCalls;
  Functor(){nCalls = 0;}

  int operator() (double x){
    nCalls++;
    cout << "Calling from ()" << endl;
    return floor(x);
  }
};

int main(){
  Functor f;
  cout << f(2.5) << endl; 
  cout << f.nCalls << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//pointToFunc.cpp
#include <iostream>
#include <cmath>

using std::cout; using std::endl;

struct Z{
  int func(double x){
    cout << "Calling from func" << endl;
    return ceil(x);
  }
};

struct Functor{
  int (*pF)(double x);

  Functor(int (*pf)(double x)){
    pF = pf;
  }

  int operator()(double x){
    cout << "Calling from ()" << endl;
    return pF(x);
  }
};

int func1(double x){
  cout << "Calling from func1" << endl;
  return floor(x);
}

int func2(double x){
  cout << "Calling from func2" << endl;
  return ceil(x);
}
  
int main(){
  Functor f1(&func1);
  cout << f1(2.5) << endl;
  Functor f2(&func2);
  cout << f2(2.5) << endl;
  Z z;
  //Functor f3(&z.func);
  //cout << f3(2.5) << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//pointToMemFunc.cpp
#include <iostream>
#include <cmath>

using std::cout; using std::endl;

class Z{

public:
  int func1(double x){
    cout << "Calling from Z func1" << endl;
    return floor(x);
  }

 int func2(double x){
   cout << "Calling from Z func2" << endl;
   return ceil(x);
  }
};

struct Functor{
  Z* pZ;
  int (Z::*pF)(double x);
  
  Functor(Z* pz, int (Z::*pf)(double x)){
    pZ = pz; pF = pf;
  }

  int operator() (double x) const {
    return (pZ->*pF)(x);
  }
};

int main(){
  Z z;
  Functor f1(&z, &Z::func1);
  cout << f1(2.5) << endl;
  Functor f2(&z, &Z::func2);
  cout << f2(2.5) << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//funcTemp.h
template <class T> class Functor{

  int (T::*pF)(double);
  T* pT;

public:
  Functor(){}
  Functor(T* pt, int (T::*pf)(double)){pT = pt; pF = pf;}
  int operator() (double u) const {
    return (pT->*pF)(u);
  }
};

//functorDriver.cpp
#include <iostream>
#include <cmath>
#include "funcTemp.h"

using std::cout; using std::endl;

struct Y{

  int func(double x){
    cout << "Calling from Y func" << endl;
    return ceil(x);
  }
};

class Z{

public:
  int func(double x){
    cout << "Calling from Z func" << endl;
    return floor(x);
  }
};

int main(){
  Y y; Z z;
  Functor<Y> f1(&y, &Y::func);
  cout << f1(2.5) << endl;
  Functor<Z> f2(&z, &Z::func);
  cout << f2(2.5) << endl;
  return 0;
}
///////////////////////////////////////////////////////////////////////
//optim.h
#ifndef OPTIM_H
#define OPTIM_H

class Functor;//forward declaration

class Optim{
  int itMax;
  double delta;

public:

  Optim(int ItMax = 38, double delta = 0.);
  
  double golden(Functor& f, double low, double up) const;
  double ranGolden(Functor& f, double a, double b, int nSearch, 
        unsigned long seed) const;
 };
#endif

//optim.cpp
#include <cmath>
#include <limits>
#include <algorithm>
#include "optim.h"
#include "functor.h"
#include "ranGen.h"

const double rup = (-1. + sqrt(5.))/2;
const double rlow = 1 - rup;

Optim::Optim(int ItMax, double Delta){
    itMax = ItMax;
    delta = Delta;
}

double Optim::golden(Functor& f, double low, double up) const {
  int niter = 0;
  double minTempLast = low, minTemp = up;
  double flow = f(low + rlow*(up - low));
  double fup = f(low + rup*(up - low));
  double change = std::numeric_limits<double>::infinity();

  while(change > delta*fabs(minTempLast) && niter < itMax){
    minTempLast = minTemp;
    minTemp = (low + up)/2.;
    if(flow < fup){
      up = low + rup*(up - low), fup = flow;
      flow = f(low + rlow*(up - low));
    }
    else{
      low = low + rlow*(up - low), flow = fup;
      fup = f(low + rup*(up - low));
    }
    change = fabs(minTemp - minTempLast);
    niter++;
  }
  return minTemp;
}

double Optim::ranGolden(Functor& f, double a, double b, int nSearch, 
        unsigned long seed) const {
  double newMinTemp, minTemp = b;
  double up, low;

  //instantiate ranGen object 
  ranGen RNG;
  RNG.setSeed(seed);
  //pointer for random uniforms
  double* pu = new double[2];
  
  //random search begins
  for(int i = 0; i < nSearch ; i++){
    //generate uniforms for the search interval
    RNG.ranUnif(2, pu);

    up = a + (b - a)*std::max(pu[0], pu[1]); 
    low = a + (b - a)*std::min(pu[0],pu[1]);
    
    newMinTemp = golden(f, low, up);

    //compare local to global minimizer
    if(f(newMinTemp) < f(minTemp)){
      minTemp = newMinTemp;
    }
  }
  delete[] pu;
  return minTemp;
}

//functor.h
#ifndef FUNCTOR_H
#define FUNCTOR_H

struct Functor{
  double (*pF)(double theta);

  Functor(double (*pf)(double theta)){
    pF = pf;
  }

  double operator()(double theta) const {
     return pF(theta);
  }
};
#endif

//optDriver.cpp
#include <iostream>
#include <cstdlib>
#include <iomanip>
#include "optim.h"
#include "functor.h"

double fun(double theta){
  return (-theta*(1. - theta));
}

int main(int argc, char** argv){
  int itMax = atoi(argv[1]);
  double delta = atof(argv[2]);
  Functor f(&fun);
  Optim Opt = Optim(itMax, delta); 
  double fMin = Opt.golden(f, 0., 1.);
  std::cout << "The approximate minimizer is " 
            << std::setprecision(16) << fMin << std::endl;

  return 0;
}
///////////////////////////////////////////////////////////////////////
//To recreate the Section 8.3.2 Example, do the following:
// i) take the matrix.h, matrix.cpp, vector.h and vector.cpp files and 
//     change the extensions to cc and hh. The include directives in
//     the files must be altered to reflect this
// ii) use the template version of the Optim class (below) with hh file extesion  
// iii) use one (e.g., the one with derivatives) of the template functor structs with
//      hh file extension

//linss.hh
#ifndef LINSS_H
#define LINSS_H

#include "vector.hh"
#include "matrix.hh"

class Linss{
  int n;
  Vector x, y;
  Matrix lhsMat(double lam) const;

public:
  Linss(const Vector& x, const Vector& Y);
  Vector smooth(double lam, double* GCV) const;
  double gcv(double lLam) const;
  Vector smooth(double* lam, double* GCV, double a = -5., 
        double b = 5., int nGridPts = 100) const; 
}; 
#endif

//linss.cc
#include "linss.hh"
#include "optim.hh"
#include "functor.hh"
#include <iostream>
#include <cmath>

Linss::Linss(const Vector& X, const Vector& Y){
  x = Vector(X);
  y = Vector(Y);
  n = x.getnRows();
}

Vector Linss::smooth(double lam, double* GCV) const {
  Matrix G(n, n, 0.);
  Matrix Hinv = lhsMat(lam);
  Vector fit = Hinv.bandChol(y, G, 1);

  //compute the GCV criterion
  double temp = 1/(G[n-1][n-1]*G[n-1][n-1]);
  double tr = temp;
  *GCV = (y[n-1] - fit[n-1])*(y[n-1] - fit[n-1]);
  for(int i = (n-2); i >= 0; i--){    
    temp = 1/(G[i][i]*G[i][i])*(1 + G[i+1][i]*G[i+1][i]*temp);
    *GCV += (y[i] - fit[i])*(y[i] - fit[i]);
    tr += temp;
  }
  *GCV /= (1. - tr/((double)n))*(1. - tr/((double)n))*(double)n;

  return fit;
}

double Linss::gcv(double lLam) const {
  double GCV = 0;
  Vector fit = smooth(pow(10, lLam), &GCV);
  return GCV;
}
 
Matrix Linss::lhsMat(double lam) const {
  double** pH = new double*[n];
  pH[0] = new double[n];
  pH[n - 1] = new double[n];
  pH[0][0] = 1. + (lam/(x[1] - x[0]));
  pH[0][1] = -lam/(x[1] - x[0]);
  pH[n - 1][n - 2] = -lam/(x[n - 1] - x[n - 2]);
  pH[n - 1][n - 1] = 1. + (lam/(x[n - 1] - x[n - 2]));
  for(int i = 1; i < (n - 1); i++){
    pH[i] = new double[n];
    pH[i][i - 1] = -lam/(x[i] - x[i - 1]);
    pH[i][i] = 1 + (lam/(x[i] - x[i - 1])) 
      + (lam/(x[i + 1] - x[i]));
    pH[i][i + 1] = -lam/(x[i + 1] - x[i]);
  }
  Matrix H(n, n, pH);

  for(int i = 0; i < n; i++)
    delete[] pH[i];
  delete[] pH;

  return H;
}

Vector Linss::smooth(double* lam, double* GCV, double a, double b, 
        int nGridPts) const {
  Functor<const Linss> F(this, &Linss::gcv);
  Optim<const Linss> Opt; 
  //initial grid search
  double gcvMin = Opt.gridSearch(F, a, b, nGridPts);

  //golden section in the neighborhood of the minimizer
  double halfRange = (b - a)/((double)nGridPts);
  double low = gcvMin - halfRange;
  double up = gcvMin + halfRange;
  gcvMin = Opt.golden(F, low, up);
  *lam = pow(10, gcvMin);
  //now compute the fit at the GCV optimal lambda
  Vector fit = smooth(*lam, GCV);
  return fit;
}
  
//linssDriver.cc
#include "R.h" 
#include "vector.hh"
#include "linss.hh"

extern "C" {
  
  void linss(int* n, double* lam, double* px, double* py, 
        double* pfit, double* GCV){

    Vector y(*n, py);
    Vector x(*n, px);
  
    Linss S(x, y);
    Vector fit = S.smooth(lam, GCV);
    for(int i = 0; i < *n; i++)
      pfit[i] = fit[i];

    Rprintf("The value of GCV at lambdaOpt = %f is %f\n", *lam, 
	    *GCV);
  }
}
///////////////////////////////////////////////////////////////////////
double Optim::newton(Functor& f, double a, double b, 
        double start) const {
  int niter = 0;
  double minTemp = start, minTempLast = b;
 
  while(fabs(minTemp - minTempLast) > delta*fabs(minTempLast) 
        && niter < itMax){
    minTempLast = minTemp;
    try{
    if(f(minTempLast, 2) == 0)
      throw std::overflow_error("Division by zero: ");    
    minTemp = minTempLast - f(minTempLast, 1)/f(minTempLast, 2);
    }

    catch(std::overflow_error e){

      std::cout << e.what() << std::endl;
      std::cout << "Enter another starting value (q to quit)" 
		<< std::endl;
      std::cin >> minTemp;
      minTempLast = b;
      continue;
    }
    
    try{
      if(minTemp > b || minTemp < a)
	throw std::out_of_range("Value outside [a, b]: ");    
    }

    catch(std::out_of_range e){
      std::cout << e.what() << std::endl;

      std::cout << "Enter another starting value (q to quit)" 
		<< std::endl;
      std::cin >> minTemp;
      minTempLast = b;
      continue;
    }    
    niter++;
  }
  return minTemp;
}
///////////////////////////////////////////////////////////////////////
//functor.h
#ifndef FUNCTOR_H
#define FUNCTOR_H

struct Functor{
  double (*pF)(double theta);
  double (*pF1)(double theta);
  double (*pF2)(double theta);

  Functor(double (*pf)(double theta), double (*pf1)(double theta), 
        double (*pf2)(double theta)){
    pF = pf; pF1 = pf1; pF2 = pf2;
  }

  Functor(double (*pf)(double theta)){
    pF = pf;
  }

  double operator()(double theta, int i) const {
    if(i == 0) return pF(theta);
    if(i == 1) return pF1(theta);
    return pF2(theta);
  }

  double operator()(double theta) const {
    return pF(theta);
  }
};
#endif
///////////////////////////////////////////////////////////////////////
//optDriver.cpp
#include <iostream>
#include <cstdlib>
#include "optim.h"
#include "functor.h"

double fun(double theta){
  return (-theta*(1. - theta));
}

double fun1(double theta){
  return (-1. + 2.*theta);
}

double fun2(double theta){
  return 2.;
}

int main(int argc, char** argv){
  int itMax = atoi(argv[1]);
  double delta = atof(argv[2]); 
  double start = atof(argv[3]); 
  Functor f(&fun, &fun1, &fun2);
  Optim Opt = Optim(itMax, delta); 

  double fMin = Opt.newton(f, 0., 1., start);
  double fMinGold = Opt.golden(f, 0., 1.);
 
  std::cout << "The Newton approximate minimizer is " 
            << fMin << std::endl;  
  std::cout << "The golden section approximate minimizer is " 
            << fMinGold << std::endl;

  return 0;
}
///////////////////////////////////////////////////////////////////////    
//optDriverTrig.cpp
#include <iostream>
#include <cstdlib>
#include <cmath>
#include "optim.h"
#include "functor.h"

double pi = 2*acos(0.); 

double fun(double theta){
  return theta*sin(4*pi*theta*theta);
}

double fun1(double theta){
  double temp = sin(4*pi*theta*theta) +
    8*pi*theta*theta*cos(4*pi*theta*theta);
  return temp;
}

double fun2(double theta){
  double temp = 24*pi*theta*cos(4*pi*theta*theta) -
    64*pi*pi*theta*theta*fun(theta);
  return temp;
}

int main(int argc, char** argv){

  int itMax = atoi(argv[1]);
  int nSearch = atoi(argv[2]);
  int seed = atoi(argv[3]);
  double delta = atof(argv[4]); 
  Functor f(&fun, &fun1, &fun2);
  Optim Opt = Optim(itMax, delta); 
  int nFail = 0;

  double fMin = Opt.ranNewton(f, 0., 1., nSearch, seed, &nFail);
  std::cout << "The approximate minimizer is " 
            << fMin << std::endl;
   std::cout << "The number of failures is " << nFail << std::endl;

  return 0;
}
 
opt : optDriverTrig.o optim.o ranGen.o
	g++ -Wall optDriverTrig.o optim.o ranGen.o -o opt

ranGen.o : ranGen.cpp ranGen.h 
	g++ -Wall -c ranGen.cpp

optim.o : optim.cpp optim.h ranGen.h functor.h
	g++ -Wall -c optim.cpp

optDriverTrig.o : optDriverTrig.cpp optim.h functor.h
	g++ -Wall -c optDriverTrig.cpp
///////////////////////////////////////////////////////////////////////    
//To recreate the mle example use the template version of Optim below

//functorIn.h
#ifndef FUNCTORIN_H
#define FUNCTORIN_H

struct FunctorIn{
  double (*pF )(double x, double theta);
  double (*pF1)(double x, double theta);
  double (*pF2)(double x, double theta);

  FunctorIn(){}
  FunctorIn(double (*pf)(double x, double theta), 
        double (*pf1)(double x, double theta), 
        double (*pf2)(double x, double theta)){
    pF = pf; pF1 = pf1; pF2 = pf2;
  }

  double operator()(double x, double theta, int i) const {
    if(i == 0) return pF(x, theta);
    if(i == 1) return pF1(x, theta);
    return pF2(x, theta);
  }
};
#endif

//like.h
#ifndef LIKE_H
#define LIKE_H
#include "functorIn.h"

class Like{
  int n;
  FunctorIn f;
  double* pxO;
  double l(double theta) const; 
  double lp(double theta) const; 
  double lpp(double theta) const;

public:

  Like(int N, char* fname, FunctorIn& F);
       
  double mle(double a, double b, double start, int itMax = 38, 
	double delta = 0.) const ;
};
#endif

//like.cpp
#include "like.h"
#include "functor.h"
#include "optim.h"
#include "fileIn.h"
     
Like::Like(int N, char* fname, FunctorIn& F){
  n = N; f = F;
  pxO = new double[n];
  fileIn<double> inObj(fname);
  inObj.read(n, pxO);  
}
  
double Like::l(double theta) const {
  double logL = 0;
  for(int i = 0; i < n; i++)
    logL += f(pxO[i], theta, 0);
  return (-logL);
}
    
double Like::lp(double theta) const {
  double logLp = 0;
  for(int i = 0; i < n; i++)
    logLp += f(pxO[i], theta, 1);
  return (-logLp);
}
      
double Like::lpp(double theta) const {
  double logLpp = 0;
  for(int i = 0; i < n; i++)
    logLpp += f(pxO[i], theta, 2);
  return (-logLpp);
}

double Like::mle(double a, double b, double start, int itMax, 
	double delta) const {
  int ier = 0;
  Functor<const Like> logL(this, &Like::l, &Like::lp, &Like::lpp);
  Optim<const Like> Opt(itMax, delta);
  double thetaHat = Opt.newton(logL, a, b, start, &ier);
  return thetaHat;
}

//functor.h
#ifndef FUNCTOR_H
#define FUNCTOR_H

template <class T> 
class Functor{

  double (T::*pF)(double);
  double (T::*pF1)(double);
  double (T::*pF2)(double);
  T* pT;

public:
  Functor(){}

  Functor(T* pt, double (T::*pf)(double)){
    pT = pt; pF = pf;
  }

 Functor(T* pt, double (T::*pf)(double), double (T::*pf1)(double),
        double (T::*pf2)(double)){
    pT = pt; pF = pf; pF1 = pf1; pF2 = pf2;
  }

 double operator() (double u) const {
    return (pT->*pF)(u);
 }

  double operator() (double u, int i) const {
    if(i == 0) return (pT->*pF)(u);
    if(i == 1) return (pT->*pF1)(u);
    return (pT->*pF2)(u);
  }
};
#endif
/////////////////////////////class Optim///////////////////////////////
//optim.h
#ifndef OPTIM_H
#define OPTIM_H

class Functor;//forward declaration

class Optim{
  int itMax;
  double delta;

public:

  Optim(int ItMax = 38, double delta = 0.);

  double golden(Functor& f, double low, double up) const;
  double ranGolden(Functor& f, double a, double b, int nSearch, 
	unsigned long seed) const;


  double newton(Functor& f, double a, double b, double start) const;
  double newton(Functor& f, double a, double b, double start, 
	int* ier) const;
  double ranNewton(Functor& f, double a, double b, int nSearch, 
        unsigned long seed, int* nFail) const;

  double ranSearch(Functor&F, double a, double b, 
	unsigned long seed) const;
  double sAnn(Functor& F, double a, double b, double d, 
	double rho, unsigned  long seed) const;
};
#endif

//optim.cpp
#include <iostream>
#include <cstdlib>
#include <cmath>
#include <stdexcept>
#include <algorithm>
#include "optim.h"
#include "functor.h"
#include "ranGen.h"

const double rup = (-1. + sqrt(5.))/2;
const double rlow = 1 - rup;

Optim::Optim(int ItMax, double Delta){
    itMax = ItMax;
    delta = Delta;
}

double Optim::golden(Functor& f, double low, double up) const {
  int niter = 0;
  double minTempLast = low, minTemp = up;
  double flow = f(low + rlow*(up - low));
  double fup = f(low + rup*(up - low));
  double change = std::numeric_limits<double>::infinity();

  while(change > delta*fabs(minTempLast) && niter < itMax){
    minTempLast = minTemp;
    minTemp = (low + up)/2.;
    if(flow < fup){
      up = low + rup*(up - low), fup = flow;
      flow = f(low + rlow*(up - low));
    }
    else{
      low = low + rlow*(up - low), flow = fup;
      fup = f(low + rup*(up - low));
    }
    change = fabs(minTemp - minTempLast);
    niter++;
  }
  return minTemp;
}

double Optim::ranGolden(Functor& f, double a, double b, int nSearch, 
        unsigned long seed) const {
  double newMinTemp, minTemp = b;
  double up, low;

  //instantiate ranGen object 
  ranGen RNG;
  RNG.setSeed(seed);
  //pointer for random uniforms
  double* pu = new double[2];
  
  //random search begins
  for(int i = 0; i < nSearch ; i++){
    //generate uniforms for the search interval
    RNG.ranUnif(2, pu);

    up = a + (b - a)*std::max(pu[0], pu[1]); 
    low = a + (b - a)*std::min(pu[0],pu[1]);
    
    newMinTemp = golden(f, low, up);

    //compare local to global minimizer
    if(f(newMinTemp, 0) < f(minTemp, 0)){
      minTemp = newMinTemp;
    }
  }
  delete[] pu;
  return minTemp;
}

double Optim::newton(Functor& f, double a, double b, 
	double start) const {
  int niter = 0;
  double minTemp = start, minTempLast = b;
 
  while(fabs(minTemp - minTempLast) > delta*fabs(minTempLast) 
        && niter < itMax){
    minTempLast = minTemp;
    try{
    if(f(minTempLast, 2) == 0)
      throw std::overflow_error("Division by zero: ");    
    minTemp = minTempLast - f(minTempLast, 1)/f(minTempLast, 2);
    }

    catch(std::overflow_error e){
      std::cout << e.what() << std::endl;

      std::cout << "Enter another staring value (q to quit)" 
		<< std::endl;
      std::cin >> minTemp;
      minTempLast = b;
      continue;
    }
    
    try{
      if(minTemp > b || minTemp < a)
	throw std::out_of_range("Value outside [a, b]: ");    
    }

    catch(std::out_of_range e){
      std::cout << e.what() << std::endl;

      std::cout << "Enter another staring value (q to quit)" 
		<< std::endl;
      std::cin >> minTemp;
      minTempLast = b;
      continue;
    }    
    niter++;
  }
  return minTemp;
}

double Optim::newton(Functor& f, double a, double b, double start, 
	int* ier) const {
  int niter = 0;
  *ier = 0;
  double minTemp = start, minTempLast = b;
 
  while(fabs(minTemp - minTempLast) > delta*fabs(minTempLast) 
        && niter < itMax){
    minTempLast = minTemp;
    if(f(minTempLast, 2) != 0)
      minTemp = minTempLast - f(minTempLast, 1)/f(minTempLast, 2);
    else{
      *ier = 1;
      break;
    }
    if(minTemp > b || minTemp < a){
      *ier = 2;
      break;
    }
    niter++;
  }
  return minTemp;
}

double Optim::ranNewton(Functor& f, double a, double b, 
	int nSearch, unsigned long seed, int* nFail) const {
  double start, newMinTemp, minTemp = b;
  int ier = 0;
  *nFail = 0;

  //instantiate ranGen object 
  ranGen RNG;
  RNG.setSeed(seed);
  //pointer for random uniforms
  double* pu = new double;

  //random search begins
  for(int i = 0; i < nSearch ; i++){
    //generate a uniform for the starting value
    RNG.ranUnif(1, pu);
    start = a + (b - a)*(*pu);
    newMinTemp = newton(f, a, b, start, &ier);

    //compare local to global minimizer
    if(f(newMinTemp, 0) < f(minTemp, 0) && ier == 0)
      minTemp = newMinTemp;
 
    else if(ier != 0) 
      ++*nFail;
  }
  delete pu;
  return minTemp;
}

double Optim::ranSearch(Functor& f, double a, double b, 
	unsigned long seed) const {
  double minTemp, U;
  int niter = 0;
  //instantiate ranGen object 
  ranGen RNG;
  RNG.setSeed(seed);
  //pointer for random uniform
  double* pu = new double;
  RNG.ranUnif(1, pu); 
  minTemp = a + (b - a)*(*pu);

  //random search begins
  while(niter < itMax){
    RNG.ranUnif(1, pu);
    U = a + (b - a)*(*pu);
    if(f(U) < f(minTemp))
      minTemp = U;

    niter++;
  }
  delete pu;
  return minTemp;
}

double Optim::sAnn(Functor& f, double a, double b, double d, 
        double rho, unsigned long seed) const {
  double minTemp, U, low, up, prob;
  int niter = 0;

  ranGen RNG;
  RNG.setSeed(seed);

  double* pu = new double;
  RNG.ranUnif(1, pu); 
 
  //random search begins
  double T = 1.;
  while(niter < itMax){
    T *= rho;
    RNG.ranUnif(1, pu);
    low = std::max(minTemp - d, a);
    up = std::min(minTemp + d, b);
    U = low + (up - low)*(*pu);
    //always accept if the objective function is decreased
    if(f(U) < f(minTemp)){
      minTemp = U;
    }
    else{//the annealing step
      prob = exp((f(minTemp) - f(U))/T);
      RNG.ranUnif(1, pu);
      if(*pu < prob)
	minTemp = U;
    }
    niter++;
  }
  delete pu;
  return minTemp;      
}
//////////////////class Optim: template version////////////////////////
//optim.h
#ifndef OPTIM_H
#define OPTIM_H
#include <cstdlib>
#include <cmath>
#include <iostream>
#include "functor.h"

template <class T> class Optim{
  int itMax;
  double delta;

public:

  Optim(int ItMax = 38, double Delta = 0.){
    itMax = ItMax; delta = Delta;
  }

  double golden(Functor<T>& f, double low, double up) const;
  double gridSearch(Functor<T>& f, double low, double up, 
	int nGridPts) const;
  double newton(Functor<T>& f, double a, double b, 
        double start, int* ier);
};

const double rup = (-1. + sqrt(5.))/2;
const double rlow = 1 - rup;

template<class T> double Optim<T>::golden(Functor<T>& f, double low, 
        double up) const {
  int niter = 0;
  double minTempLast = low, minTemp = up;
  double flow = f(low + rlow*(up - low));
  double fup = f(low + rup*(up - low));
  double change = std::numeric_limits<double>::infinity();

  while(change > delta*fabs(minTempLast) && niter < itMax){
    minTempLast = minTemp;
    minTemp = (low + up)/2.;
    if(flow < fup){
      up = low + rup*(up - low), fup = flow;
      flow = f(low + rlow*(up - low));
    }
    else{
      low = low + rlow*(up - low), flow = fup;
      fup = f(low + rup*(up - low));
    }
    change = fabs(minTemp - minTempLast);
    niter++;
  }
  return minTemp;
}

template<class T> double Optim<T>::gridSearch(Functor<T>& f, 
	double low, double up, int nGridPts) const {
  double minTemp, fTemp, gridVal, range = (up - low);
  double minF = std::numeric_limits<double>::infinity();
  for(int i = 0; i < nGridPts; i++){
    gridVal = low + range*((double)i)/((double)nGridPts);
    fTemp = f(gridVal);
    if(fTemp < minF){
      minF = fTemp;
      minTemp = gridVal;
    }
  }
  return minTemp;
}

template <class T> double Optim<T>::newton(Functor<T>& f, double a, 
        double b, double start, int* ier){
  int niter = 0;
  *ier = 0;
  double minTemp = start, minTempLast = b; 

  while(fabs(minTemp - minTempLast) > delta*fabs(minTempLast) 
        && niter < itMax){
    minTempLast = minTemp;
    if(f(minTempLast, 2) != 0)
      minTemp = minTempLast - f(minTempLast, 1)/f(minTempLast, 2);
    else{
      *ier = 1;
      break;
    }
    if(minTemp > b || minTemp < a){
      *ier = 2;
      std::cout  << minTemp << std::endl;
      break;
    }
    niter++;
  }
  if(*ier == 0){
    return minTemp;
  }
  else
    return start;
}
#endif