levrosen.cpp

Go to the documentation of this file.

#include <iostream>
#include <cmath>
#include "lm.h"

using namespace std;

void misfit(double *p, double *x, int m, int n, void *data)
  {

    x[0] = 100 * pow(p[1] - pow(p[0], 2), 2) + pow(1 - p[0], 2);
    x[1] = 100 * pow(p[1] - pow(p[0], 2), 2) + pow(1 - p[0], 2);
    x[2] = 100 * pow(p[1] - pow(p[0], 2), 2) + pow(1 - p[0], 2);
  }

int main()
  {
    const int nparams = 2;
    const int n = 3;
    const int maxiter = 100;
    double p[nparams], lb[nparams], ub[nparams]; // 6 is max(2, 3, 5)
    double *x; // 16 is max(2, 3, 5, 6, 16)
    x = new double[n];
    double opts[LM_OPTS_SZ], info[LM_INFO_SZ];

    opts[0] = LM_INIT_MU;
    opts[1] = 1E-15;
    opts[2] = 1E-15;
    opts[3] = 1E-20;
    opts[4] = LM_DIFF_DELTA; // relevant only if the finite difference jacobian version is used

    for (int i = 0; i < nparams; ++i)
      {
        p[i] = -1;
        lb[i] = -10000;
        ub[i] = 10000;
      }

    for (int i = 0; i < n; i++)
      x[i] = 0.0;

    double ret = dlevmar_bc_dif(misfit, p, x, nparams, n, lb, ub, maxiter,
        opts, info, NULL, NULL, NULL); // no jacobian
    cout << "Levenberg-Marquardt returned " << ret << " in " << info[5]
        << "iter, reason " << info[6] << endl;

    cout << "p: ";
    for (int i = 0; i < nparams; ++i)
      cout << p[i] << " ";
    cout << endl;
    cout << endl << " Minimization info:" << endl;
    for (int i = 0; i < LM_INFO_SZ; ++i)
      cout << info[i] << " ";
    cout << endl;

    return 0;
  }