Skip to content
Snippets Groups Projects
minimize.h 2.9 KiB
Newer Older
Mark Woolrich's avatar
Mark Woolrich committed
/*  minimize
 
    Tim Behrens, FMRIB Image Analysis Group

    Copyright (C) 1999-2000 University of Oxford  */

/*  CCOPYRIGHT  */

#if !defined(minimize_h)
#define minimize_h

#include <string>
#include <iostream>
#include <fstream>
David Flitney's avatar
David Flitney committed
//#include <unistd.h>
Mark Woolrich's avatar
Mark Woolrich committed
#include <vector>
#include <algorithm>
#include "newmatap.h"
#include "newmatio.h"
Mark Woolrich's avatar
Mark Woolrich committed


#define WANT_STREAM
#define WANT_MATH

using namespace MISCMATHS;
using namespace NEWMAT;
using namespace std;
///////////////////////////////////////////////////////

//fminsearch.m

namespace MISCMATHS {

class pair_comparer
{
public:
  bool operator()(const pair<float,ColumnVector>& p1,const pair<float,ColumnVector>& p2) const
  {
    return p1.first < p2.first;
  }    
};

Mark Woolrich's avatar
Mark Woolrich committed
  class EvalFunction;
  class gEvalFunction;

float diff1(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff derivative

float diff2(const ColumnVector& x, const EvalFunction& func, int i,float h,int errorord=4);// finite diff 2nd derivative

float diff2(const ColumnVector& x, const EvalFunction& func, int i,int j,float h,int errorord=4);// finite diff cross derivative

ReturnMatrix gradient(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff derivative vector 

Tim Behrens's avatar
Tim Behrens committed
ReturnMatrix hessian(const ColumnVector& x, const EvalFunction& func,float h,int errorord=4);// finite diff hessian
Tim Behrens's avatar
Tim Behrens committed
void minsearch(ColumnVector& x, const EvalFunction& func, ColumnVector& paramstovary);
Mark Woolrich's avatar
Mark Woolrich committed
void scg(ColumnVector& x, const gEvalFunction& func, ColumnVector& paramstovary, float tol = 0.0000001, float eps=1e-16, int niters=500);
Mark Woolrich's avatar
Mark Woolrich committed
class EvalFunction
Mark Woolrich's avatar
Mark Woolrich committed
{//Function where gradient is not analytic (or you are too lazy to work it out) (required for fminsearch)
Mark Woolrich's avatar
Mark Woolrich committed
public:
  EvalFunction(){}
Tim Behrens's avatar
Tim Behrens committed
  virtual float evaluate(const ColumnVector& x) const = 0; //evaluate the function
Mark Woolrich's avatar
Mark Woolrich committed
  virtual ~EvalFunction(){};

Mark Woolrich's avatar
Mark Woolrich committed
  virtual void minimize(ColumnVector& x)
  {
Tim Behrens's avatar
Tim Behrens committed
    ColumnVector paramstovary(x.Nrows());
    paramstovary = 1;
    minsearch(x,*this,paramstovary);
Tim Behrens's avatar
Tim Behrens committed
  virtual void minimize(ColumnVector& x, ColumnVector& paramstovary)
Tim Behrens's avatar
Tim Behrens committed
    minsearch(x,*this,paramstovary);
Mark Woolrich's avatar
Mark Woolrich committed
private:
Mark Woolrich's avatar
Mark Woolrich committed
  const EvalFunction& operator=(EvalFunction& par);
  EvalFunction(const EvalFunction&);
};

Tim Behrens's avatar
Tim Behrens committed
class gEvalFunction : public EvalFunction
{//Function where gradient is analytic (required for scg)
public:
  gEvalFunction() : EvalFunction(){}
  // evaluate is inherited from EvalFunction
Mark Woolrich's avatar
Mark Woolrich committed
  
  virtual ReturnMatrix g_evaluate(const ColumnVector& x) const = 0; //evaluate the gradient
Tim Behrens's avatar
Tim Behrens committed
  virtual ~gEvalFunction(){};

Mark Woolrich's avatar
Mark Woolrich committed
  virtual void minimize(ColumnVector& x)
  {
Tim Behrens's avatar
Tim Behrens committed
    ColumnVector paramstovary(x.Nrows());
    paramstovary = 1;
    scg(x,*this,paramstovary);
Tim Behrens's avatar
Tim Behrens committed
  virtual void minimize(ColumnVector& x, ColumnVector& paramstovary)
Tim Behrens's avatar
Tim Behrens committed
    scg(x,*this,paramstovary);
Tim Behrens's avatar
Tim Behrens committed
private:
    
  const gEvalFunction& operator=(gEvalFunction& par);
  gEvalFunction(const gEvalFunction&);
};

}