Constrained Levenberg-Marquardt using the InMin library

// Bojan Nikolic <bojan@bnikolic.co.uk>
//
// Example of constrained Levenberg-Marquardt optimisiation using the
// inmin library http://www.bnikolic.co.uk/inmin/inmin-library.html

#include <vector>
#include <cmath>
#include <iostream>

#include <boost/assign/list_of.hpp>

#include "inmin_lm.h"

/** A toy second-degree polynomial model
 */
struct PolyModel {

  /// The model parameters
  double a, b, c;

  /// These are the observation points, i.e., the points at which the
  /// model should be evaluated
  std::vector<double> x;

  // The observed points
  std::vector<double> obs;

  PolyModel(const std::vector<double> &x):
    a(1.0),
    b(2.0),
    c(3.0),
    x(x)
  {
  }

  void operator() (std::vector<double> &res)
  {
    const size_t n=x.size();
    res.resize(n);
    for(size_t i=0; i<n; ++i)
      res[i]=a+b*x[i]+c*pow(x[i],2);
  }
};


/// This function interfaces between the imnin library and simple C++
/// code. It is hard coded for this particular problem.
void LMHelper(const double *x,
              double *res,
              void *data)
{
  PolyModel *pmodel=reinterpret_cast<PolyModel*>(data);
  pmodel->a=x[0];
  pmodel->b=x[1];
  pmodel->c=x[2];

  std::cout<<">> Interim point: "
           <<x[0]<<","
           <<x[1]<<","
           <<x[2]<<std::endl;


  std::vector<double> scratch;

  (*pmodel)(scratch);

  for (size_t i=0; i<pmodel->obs.size(); ++i)
  {
    res[i]=scratch[i]-pmodel->obs[i];
  }
}

void minimise(PolyModel &model,
              std::vector<double> &obs,
              std::vector<double> &lowerbounds,
              std::vector<double> &upperbounds)
{

  inmin_lm_in in;
  memset(&in, 0, sizeof(inmin_lm_in));
  in.N=3;
  in.M=obs.size();
  in.f=LMHelper;

  in.ftol=1e-8;
  in.xtol=1e-8;
  in.gtol=1e-8;
  in.maxfev=100;

  model.obs=obs;

  double box[6];
  for (size_t i=0; i<3; ++i)
  {
     box[i*2]=lowerbounds[i];
     box[i*2+1]=upperbounds[i];
  }
  in.box=box;

  std::vector<double> res(3);
  inmin_lm_out out;
  out.x=&res[0];
  out.covar=NULL;

  double ic[] = {0.0,0.0,1.0};

  inmin_lm_run(&in,
               ic,
               &out,
               &model);

  std::cout<<"The final result is:"
           <<res[0]<<","
           <<res[1]<<","
           <<res[2]<<std::endl;

};




int main(void)
{
  // This is the list of values of x at which the unknown function is
  // observed
  std::vector<double> x = boost::assign::list_of
    (-10)(-8)(-5)(1)(4)(10)(40);

  PolyModel model(x);

  // The observed values corresponding to the x coordinates above. For
  // this problem simply create mock data using the model itself
  std::vector<double> obs;
  model(obs);

  std::cout<<"Simulated observation: ";
  for (size_t i=0; i<obs.size(); ++i)
  {
    std::cout<<obs[i]<<",";
  }
  std::cout<<std::endl;

  std::vector<double> lower = boost::assign::list_of
    (-10)(-10)(-10);

  // Make the upper bound of the first model parameter smaller than
  // 1.0 to show off the constrained minimisation
  std::vector<double> upper = boost::assign::list_of
    (0.5)(10)(10);


  // Run the minimiser
  minimise(model,
           obs,
           lower,
           upper);

}