---

optimizer.h

// Emacs will be in -*- Mode: c++ -*-
//
// ************ DO NOT REMOVE THIS BANNER ****************
//
//  Nicolas Di Cesare <Nicolas.Dicesare@free.fr>
//  http://acm.emath.fr/~dicesare
//
//********************************************************
//
//  Optimization Method and Problem abstract classes
//
//********************************************************
#ifndef _optimizer_h
#define _optimizer_h

#include <iostream>
#include <iomanip>

#include <costfunction.h>
#include <criteria.h>

#include <ql/handle.hpp>

//using QuantLib::Handle;
using QuantLib::DotProduct;
using std::setw;

template <class V> class OptimizationMethod;
template <class V> class OptimizationProblem;


template <class V>
class OptimizationMethod {
public:
  typedef double value_type;
protected:
  V initialValue_;
  int iterationNumber_;
  OptimizationEndCriteria endCriteria_;
  int functionEvaluation_, gradientEvaluation_;
  value_type functionValue_, squaredNorm_;
  V x_, searchDirection_;
public:
  explicit OptimizationMethod() 
    : iterationNumber_(0),
      functionEvaluation_(0), gradientEvaluation_(0), 
      functionValue_(1), squaredNorm_(1)  {}
  virtual ~OptimizationMethod() {}

  inline void setInitialValue(const V& initialValue) 
  {
    iterationNumber_ = 0;
    initialValue_ = initialValue; 
    x_ = initialValue; 
    searchDirection_ = V(x_.size());
  } 
  inline void setEndCriteria(const OptimizationEndCriteria& endCriteria) { endCriteria_ = endCriteria;}

  inline int& iterationNumber() { return iterationNumber_;}
  inline OptimizationEndCriteria& endCriteria() { return endCriteria_;}
  inline int & functionEvaluation() { return functionEvaluation_; }
  inline int & gradientEvaluation() { return gradientEvaluation_; }
  inline value_type& functionValue() { return functionValue_; }
  inline value_type& gradientNormValue() { return squaredNorm_; }
  V& x() { return x_;}
  V& searchDirection() { return searchDirection_;}
  
  virtual void Minimize(OptimizationProblem<V> & P) = 0;
};



template <class V>
class OptimizationProblemOutput {
  std::ostream& os_;
public:
  typedef double value_type;

  OptimizationProblemOutput() : os_(std::cout) {}
  OptimizationProblemOutput(std::ostream& os) : os_(os) {}
  virtual ~OptimizationProblemOutput() {}

  virtual void init(OptimizationMethod<V>& method) 
  {
#ifdef DISPLAY_OPTI
    os_ << "Initial value : " << method.x() << std::endl;
    os_ << std::endl
    << " fen : function evaluation number" << std::endl
    << " gen : gradient evaluation number" << std::endl
    << "| iter | function value   |     gradient norm     |  line-search step     |  fen  |  gen  |" << std::endl
    << "|------|------------------|-----------------------|-----------------------|-------|-------|" << std::endl;
#endif
  }

  virtual void save(int iterationNumber, 
            value_type function, 
            value_type normGradient,
            value_type lineSearchStep,
            OptimizationMethod<V>& method)
  {
    method.functionValue() = function; method.gradientNormValue() = normGradient;
#ifdef DISPLAY_OPTI
    os_.setf(std::ios::scientific,std::ios::floatfield);
    os_ << "|" << setw(4) << iterationNumber
    << "  |" << setw(15) << function 
    << "   |" << setw(19) << normGradient
    << "    |" << setw(19) << lineSearchStep
    << "    |" << setw(5) << method.functionEvaluation() 
    << "  |" << setw(5) << method.gradientEvaluation() 
    << "  |" << std::endl;
    os_.flush();
#endif
  }
};

template <class V>
class OptimizationProblem {
public:
  typedef double value_type;
protected:
  CostFunction<V> & costFunction_;
  OptimizationMethod<V> & method_;
  QuantLib::Handle<OptimizationProblemOutput<V> > opo_;
public:
  OptimizationProblem(CostFunction<V> & f,           // Function and it gradient vector
              OptimizationMethod<V> & meth)  // Optimization method
    : costFunction_(f), method_(meth), 
      opo_(QuantLib::Handle<OptimizationProblemOutput<V> >(new OptimizationProblemOutput<V>()))
  {
    opo_->init(method_);
  }
  OptimizationProblem(CostFunction<V> & f,           // Function and it gradient vector
              OptimizationMethod<V> & meth,          // Optimization method
              QuantLib::Handle<OptimizationProblemOutput<V> >& opo)
    : costFunction_(f), method_(meth), opo_(opo)
  {
    opo_->init(method_);
  }

  ~OptimizationProblem() {}

  value_type value(const V& x)
  { 
    method_.functionEvaluation()++; 
    return costFunction_.value(x); 
  }

  void firstDerivative(V& grad_f, const V& x)
  { 
    method_.gradientEvaluation()++; 
    costFunction_.firstDerivative(grad_f, x); 
  }

  value_type valueAndFirstDerivative(V& grad_f, const V& x)
  { 
    method_.functionEvaluation()++; method_.gradientEvaluation()++; 
    return costFunction_.valueAndFirstDerivative(grad_f, x);
  }

  OptimizationMethod<V> & optimisationMethod() { return method_;}

  void Minimize()
  { 
    method_.Minimize( *this ); 
  }

  V& minimumValue() { return method_.x();}

  void Save(int iterationNumber, 
        value_type function, 
        value_type normGradient,
        value_type lineSearchStep,
        OptimizationMethod<V>& method)
  {
    opo_->save(iterationNumber, function, normGradient, lineSearchStep, method);
  }
};

#endif

---