# Using the Lourakis constrained Levenberg-Marquardt library

Here is a very quickly prepared example of how to use the Lourakis library to do constrained minimisation in n-dimensions. This shows simple fitting of a second degree polynomial to a handful of observed points. The box constraint is used and it is set so that it does not contain the maximum likelihood point -- this shows off the constraint feature.

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

#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,
{
pmodel->a=p;
pmodel->b=p;
pmodel->c=p;

std::cout<<">> Interim point: "
<<p<<","
<<p<<","
<<p<<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;

// 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,
3,
obs.size(),
&lowerbounds,
&upperbounds,
1000,
opts,
info,
NULL,
NULL,
reinterpret_cast<void*>(&model)
);

std::cout<<"The final result is:"
<<ic<<","
<<ic<<","
<<ic<<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);

}
```