Constrained Levenberg-Marquardt using the InMin library

This is the equivalent program to the own shown in this old blog post, but implemented using the new InMin library. The program demonstrates Levenberg-Marquardt fitting of a function with a box constraint on the parameters.

```// 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);

}
```