Using the Lourakis constrained Levenberg-Marquardt library

// Bojan Nikolic <bojan@bnikolic.co.uk>
//
// Example of how to use the Lourakis L-M library. Note this library
// is licensed under the GPL.

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

#include <boost/assign/list_of.hpp>

#include "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;

  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 Lourakis library and C++
/// code. It is hard coded for this particular problem. See BNMin1 for
/// a more general approach.
void LMHelper(double *p,
              double *hx,
              int m,
              int n,
              void *adata)
{
  PolyModel *pmodel=reinterpret_cast<PolyModel*>(adata);
  pmodel->a=p[0];
  pmodel->b=p[1];
  pmodel->c=p[2];

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


  std::vector<double> scratch;

  (*pmodel)(scratch);

  for (size_t i=0; i<n; ++i)
  {
    hx[i]=scratch[i];
  }
}

void minimise(PolyModel &model,
              std::vector<double> &obs,
              std::vector<double> &lowerbounds,
              std::vector<double> &upperbounds)
{
  /// Options to the solver, see the manual of the library
  double opts[]={0.1, 0.0001, 0.0001, 0.0001};

  /// Output information from the solver
  double info[9];

  // This is the starting point of the minimisation. The finate
  // difference doesn't work too great if it is started from 0,0,0
  double ic[] = {0,0,1.0};

  dlevmar_bc_dif(LMHelper,
                 ic,
                 &obs[0],
                 3,
                 obs.size(),
                 &lowerbounds[0],
                 &upperbounds[0],
                 1000,
                 opts,
                 info,
                 NULL,
                 NULL,
                 reinterpret_cast<void*>(&model)
                 );

  std::cout<<"The final result is:"
           <<ic[0]<<","
           <<ic[1]<<","
           <<ic[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.9)(10)(10);

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

}