lmdemo.c

a copylefted C/C++ implementation of the Levenberg-Marquardt non-linear least squares algorithm

/////////////////////////////////////////////////////////////////////////////////
// 
//  Demonstration driver program for the Levenberg - Marquardt minimization
//  algorithm
//  Copyright (C) 2004-05  Manolis Lourakis (lourakis at ics forth gr)
//  Institute of Computer Science, Foundation for Research & Technology - Hellas
//  Heraklion, Crete, Greece.
//
//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
/////////////////////////////////////////////////////////////////////////////////

/******************************************************************************** 
 * Levenberg-Marquardt minimization demo driver. Only the double precision versions
 * are tested here. See the Meyer case for an example of verifying the Jacobian 
 ********************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

#include "lm.h"

#ifndef LM_DBL_PREC
#error Demo program assumes that levmar has been compiled with double precision, see LM_DBL_PREC!
#endif


/* Sample functions to be minimized with LM and their Jacobians.
 * More test functions at http://www.csit.fsu.edu/~burkardt/f_src/test_nls/test_nls.html
 * Check also the CUTE problems collection at ftp://ftp.numerical.rl.ac.uk/pub/cute/;
 * CUTE is searchable through http://numawww.mathematik.tu-darmstadt.de:8081/opti/select.html
 * CUTE problems can also be solved through the AMPL web interface at http://www.ampl.com/TRYAMPL/startup.html
 *
 * Nonlinear optimization models in AMPL can be found at http://www.princeton.edu/~rvdb/ampl/nlmodels/
 */

#define ROSD 105.0

/* Rosenbrock function, global minimum at (1, 1) */
void ros(double *p, double *x, int m, int n, void *data)
{
register int i;

  for(i=0; i<n; ++i)
    x[i]=((1.0-p[0])*(1.0-p[0]) + ROSD*(p[1]-p[0]*p[0])*(p[1]-p[0]*p[0]));
}

void jacros(double *p, double *jac, int m, int n, void *data)
{
register int i, j;

  for(i=j=0; i<n; ++i){
    jac[j++]=(-2 + 2*p[0]-4*ROSD*(p[1]-p[0]*p[0])*p[0]);
    jac[j++]=(2*ROSD*(p[1]-p[0]*p[0]));
  }
}


#define MODROSLAM 1E02
/* Modified Rosenbrock problem, global minimum at (1, 1) */
void modros(double *p, double *x, int m, int n, void *data)
{
register int i;

  for(i=0; i<n; i+=3){
    x[i]=10*(p[1]-p[0]*p[0]);
	  x[i+1]=1.0-p[0];
	  x[i+2]=MODROSLAM;
  }
}

void jacmodros(double *p, double *jac, int m, int n, void *data)
{
register int i, j;

  for(i=j=0; i<n; i+=3){
    jac[j++]=-20.0*p[0];
	  jac[j++]=10.0;

	  jac[j++]=-1.0;
	  jac[j++]=0.0;

	  jac[j++]=0.0;
	  jac[j++]=0.0;
  }
}


/* Powell's function, minimum at (0, 0) */
void powell(double *p, double *x, int m, int n, void *data)
{
register int i;

  for(i=0; i<n; i+=2){
    x[i]=p[0];
    x[i+1]=10.0*p[0]/(p[0]+0.1) + 2*p[1]*p[1];
  }
}

void jacpowell(double *p, double *jac, int m, int n, void *data)
{
register int i, j;

  for(i=j=0; i<n; i+=2){
    jac[j++]=1.0;
    jac[j++]=0.0;

    jac[j++]=1.0/((p[0]+0.1)*(p[0]+0.1));
    jac[j++]=4.0*p[1];
  }
}

/* Wood's function, minimum at (1, 1, 1, 1) */
void wood(double *p, double *x, int m, int n, void *data)
{
register int i;

  for(i=0; i<n; i+=6){
    x[i]=10.0*(p[1] - p[0]*p[0]);
    x[i+1]=1.0 - p[0];
    x[i+2]=sqrt(90.0)*(p[3] - p[2]*p[2]);
    x[i+3]=1.0 - p[2];
    x[i+4]=sqrt(10.0)*(p[1]+p[3] - 2.0);
    x[i+5]=(p[1] - p[3])/sqrt(10.0);
  }
}

/* Meyer's (reformulated) problem, minimum at (2.48, 6.18, 3.45) */
void meyer(double *p, double *x, int m, int n, void *data)
{
register int i;
double ui;

	for(i=0; i<n; ++i){
		ui=0.45+0.05*i;
		x[i]=p[0]*exp(10.0*p[1]/(ui+p[2]) - 13.0);
	}
}

void jacmeyer(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double ui, tmp;

  for(i=j=0; i<n; ++i){
	  ui=0.45+0.05*i;
	  tmp=exp(10.0*p[1]/(ui+p[2]) - 13.0);

	  jac[j++]=tmp;
	  jac[j++]=10.0*p[0]*tmp/(ui+p[2]);
	  jac[j++]=-10.0*p[0]*p[1]*tmp/((ui+p[2])*(ui+p[2]));
  }
}

/* helical valley function, minimum at (1.0, 0.0, 0.0) */
#ifndef M_PI
#define M_PI   3.14159265358979323846  /* pi */
#endif

void helval(double *p, double *x, int m, int n, void *data)
{
double theta;

  if(p[0]<0.0)
     theta=atan(p[1]/p[0])/(2.0*M_PI) + 0.5;
  else if(0.0<p[0])
     theta=atan(p[1]/p[0])/(2.0*M_PI);
  else 
    theta=(p[1]>=0)? 0.25 : -0.25;

  x[0]=10.0*(p[2] - 10.0*theta);
  x[1]=10.0*(sqrt(p[0]*p[0] + p[1]*p[1]) - 1.0);
  x[2]=p[2];
}

void jachelval(double *p, double *jac, int m, int n, void *data)
{
register int i=0;
double tmp;

  tmp=p[0]*p[0] + p[1]*p[1];

  jac[i++]=50.0*p[1]/(M_PI*tmp);
  jac[i++]=-50.0*p[0]/(M_PI*tmp);
  jac[i++]=10.0;

  jac[i++]=10.0*p[0]/sqrt(tmp);
  jac[i++]=10.0*p[1]/sqrt(tmp);
  jac[i++]=0.0;

  jac[i++]=0.0;
  jac[i++]=0.0;
  jac[i++]=1.0;
}

/* Boggs - Tolle problem 3 (linearly constrained), minimum at (-0.76744, 0.25581, 0.62791, -0.11628, 0.25581)
 * constr1: p[0] + 3*p[1] = 0;
 * constr2: p[2] + p[3] - 2*p[4] = 0;
 * constr3: p[1] - p[4] = 0;
 */
void bt3(double *p, double *x, int m, int n, void *data)
{
register int i;
double t1, t2, t3, t4;

  t1=p[0]-p[1];
  t2=p[1]+p[2]-2.0;
  t3=p[3]-1.0;
  t4=p[4]-1.0;

  for(i=0; i<n; ++i)
    x[i]=t1*t1 + t2*t2 + t3*t3 + t4*t4;
}

void jacbt3(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double t1, t2, t3, t4;

  t1=p[0]-p[1];
  t2=p[1]+p[2]-2.0;
  t3=p[3]-1.0;
  t4=p[4]-1.0;

  for(i=j=0; i<n; ++i){
    jac[j++]=2.0*t1;
    jac[j++]=2.0*(t2-t1);
    jac[j++]=2.0*t2;
    jac[j++]=2.0*t3;
    jac[j++]=2.0*t4;
  }
}

/* Hock - Schittkowski problem 28 (linearly constrained), minimum at (0.5, -0.5, 0.5)
 * constr1: p[0] + 2*p[1] + 3*p[2] = 1;
 */
void hs28(double *p, double *x, int m, int n, void *data)
{
register int i;
double t1, t2;

  t1=p[0]+p[1];
  t2=p[1]+p[2];

  for(i=0; i<n; ++i)
    x[i]=t1*t1 + t2*t2;
}

void jachs28(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double t1, t2;

  t1=p[0]+p[1];
  t2=p[1]+p[2];

  for(i=j=0; i<n; ++i){
    jac[j++]=2.0*t1;
    jac[j++]=2.0*(t1+t2);
    jac[j++]=2.0*t2;
  }
}

/* Hock - Schittkowski problem 48 (linearly constrained), minimum at (1.0, 1.0, 1.0, 1.0, 1.0)
 * constr1: sum {i in 0..4} p[i] = 5;
 * constr2: p[2] - 2*(p[3]+p[4]) = -3;
 */
void hs48(double *p, double *x, int m, int n, void *data)
{
register int i;
double t1, t2, t3;

  t1=p[0]-1.0;
  t2=p[1]-p[2];
  t3=p[3]-p[4];

  for(i=0; i<n; ++i)
    x[i]=t1*t1 + t2*t2 + t3*t3;
}

void jachs48(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double t1, t2, t3;

  t1=p[0]-1.0;
  t2=p[1]-p[2];
  t3=p[3]-p[4];

  for(i=j=0; i<n; ++i){
    jac[j++]=2.0*t1;
    jac[j++]=2.0*t2;
    jac[j++]=-2.0*t2;
    jac[j++]=2.0*t3;
    jac[j++]=-2.0*t3;
  }
}

/* Hock - Schittkowski problem 51 (linearly constrained), minimum at (1.0, 1.0, 1.0, 1.0, 1.0)
 * constr1: p[0] + 3*p[1] = 4;
 * constr2: p[2] + p[3] - 2*p[4] = 0;
 * constr3: p[1] - p[4] = 0;
 */
void hs51(double *p, double *x, int m, int n, void *data)
{
register int i;
double t1, t2, t3, t4;

  t1=p[0]-p[1];
  t2=p[1]+p[2]-2.0;
  t3=p[3]-1.0;
  t4=p[4]-1.0;

  for(i=0; i<n; ++i)
    x[i]=t1*t1 + t2*t2 + t3*t3 + t4*t4;
}

void jachs51(double *p, double *jac, int m, int n, void *data)
{
register int i, j;
double t1, t2, t3, t4;

  t1=p[0]-p[1];
  t2=p[1]+p[2]-2.0;
  t3=p[3]-1.0;
  t4=p[4]-1.0;

  for(i=j=0; i<n; ++i){
    jac[j++]=2.0*t1;
    jac[j++]=2.0*(t2-t1);
    jac[j++]=2.0*t2;
    jac[j++]=2.0*t3;
    jac[j++]=2.0*t4;
  }
}

/* Hock - Schittkowski problem 01 (box constrained), minimum at (1.0, 1.0)
 * constr1: p[1]>=-1.5;
 */
void hs01(double *p, double *x, int m, int n, void *data)
{
double t;

  t=p[0]*p[0];
  x[0]=10.0*(p[1]-t);
  x[1]=1.0-p[0];
}

void jachs01(double *p, double *jac, int m, int n, void *data)
{
register int j=0;

  jac[j++]=-20.0*p[0];
  jac[j++]=10.0;

  jac[j++]=-1.0;
  jac[j++]=0.0;
}

/* Hock - Schittkowski MODIFIED problem 21 (box constrained), minimum at (2.0, 0.0)
 * constr1: 2 <= p[0] <=50;
 * constr2: -50 <= p[1] <=50;
 *
 * Original HS21 has the additional constraint 10*p[0] - p[1] >= 10; which is inactive
 * at the solution, so it is dropped here.
 */
void hs21(double *p, double *x, int m, int n, void *data)
{
  x[0]=p[0]/10.0;
  x[1]=p[1];
}

void jachs21(double *p, double *jac, int m, int n, void *data)
{
register int j=0;

  jac[j++]=0.1;
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=1.0;
}

/* Problem hatfldb (box constrained), minimum at (0.947214, 0.8, 0.64, 0.4096)
 * constri: p[i]>=0.0; (i=1..4)
 * constr5: p[1]<=0.8;
 */
void hatfldb(double *p, double *x, int m, int n, void *data)
{
register int i;

  x[0]=p[0]-1.0;

  for(i=1; i<m; ++i)
     x[i]=p[i-1]-sqrt(p[i]);
}

void jachatfldb(double *p, double *jac, int m, int n, void *data)
{
register int j=0;

  jac[j++]=1.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[1]);
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[2]);
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[3]);
}

/* Problem hatfldc (box constrained), minimum at (1.0, 1.0, 1.0, 1.0)
 * constri: p[i]>=0.0; (i=1..4)
 * constri+4: p[i]<=10.0; (i=1..4)
 */
void hatfldc(double *p, double *x, int m, int n, void *data)
{
register int i;

  x[0]=p[0]-1.0;

  for(i=1; i<m-1; ++i)
     x[i]=p[i-1]-sqrt(p[i]);

  x[m-1]=p[m-1]-1.0;
}

void jachatfldc(double *p, double *jac, int m, int n, void *data)
{
register int j=0;

  jac[j++]=1.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[1]);
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=1.0;
  jac[j++]=-0.5/sqrt(p[2]);
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=1.0;
}

/* Hock - Schittkowski (modified) problem 52 (box/linearly constrained), minimum at (-0.09, 0.03, 0.25, -0.19, 0.03)
 * constr1: p[0] + 3*p[1] = 0;
 * constr2: p[2] +   p[3] - 2*p[4] = 0;
 * constr3: p[1] -   p[4] = 0;
 *
 * To the above 3 constraints, we add the following 5:
 * constr4: -0.09 <= p[0];
 * constr5:   0.0 <= p[1] <= 0.3;
 * constr6:          p[2] <= 0.25;
 * constr7:  -0.2 <= p[3] <= 0.3;
 * constr8:   0.0 <= p[4] <= 0.3;
 *
 */
void modhs52(double *p, double *x, int m, int n, void *data)
{
  x[0]=4.0*p[0]-p[1];
  x[1]=p[1]+p[2]-2.0;
  x[2]=p[3]-1.0;
  x[3]=p[4]-1.0;
}

void jacmodhs52(double *p, double *jac, int m, int n, void *data)
{
register int j=0;

  jac[j++]=4.0;
  jac[j++]=-1.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=1.0;
  jac[j++]=1.0;
  jac[j++]=0.0;
  jac[j++]=0.0;

  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=0.0;
  jac[j++]=1.0;