1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Thomas Capricelli <orzel@freehackers.org> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <stdio.h> 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "main.h" 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <unsupported/Eigen/NonLinearOptimization> 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This disables some useless Warnings on MSVC. 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// It is intended to be done for this test only. 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include <Eigen/src/Core/util/DisableStupidWarnings.h> 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag) 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* subroutine fcn for chkder example. */ 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(15 == fvec.size()); 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(3 == x.size()); 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3, tmp4; 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iflag == 0) 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iflag != 2) 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i=0; i<15; i++) { 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16-i-1; 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = tmp1; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i >= 8) tmp3 = tmp2; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else { 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i = 0; i < 15; i++) { 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16-i-1; 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* error introduced into next statement for illustration. */ 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* corrected statement should read tmp3 = tmp1 . */ 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = tmp2; 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i >= 8) tmp3 = tmp2; 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4=tmp4*tmp4; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = -1.; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = tmp1*tmp2/tmp4; 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = tmp1*tmp3/tmp4; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testChkder() 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m=15, n=3; 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n), fvec(m), xp, fvecp(m), err; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd fjac(m,n); 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXi ipvt; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following values should be suitable for */ 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* checking the jacobian matrix. */ 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x << 9.2e-1, 1.3e-1, 5.4e-1; 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::chkder(x, fvec, fjac, xp, fvecp, 1, err); 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fcn_chkder(x, fvec, fjac, 1); 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fcn_chkder(x, fvec, fjac, 2); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fcn_chkder(xp, fvecp, fjac, 1); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::chkder(x, fvec, fjac, xp, fvecp, 2, err); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvecp -= fvec; 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check those 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd fvec_ref(m), fvecp_ref(m), err_ref(m); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec_ref << 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -1.181606, -1.429655, -1.606344, 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -1.745269, -1.840654, -1.921586, 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -1.984141, -2.022537, -2.468977, 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -2.827562, -3.473582, -4.437612, 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -6.047662, -9.267761, -18.91806; 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvecp_ref << 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -7.724666e-09, -3.432406e-09, -2.034843e-10, 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2.313685e-09, 4.331078e-09, 5.984096e-09, 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 7.363281e-09, 8.53147e-09, 1.488591e-08, 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2.33585e-08, 3.522012e-08, 5.301255e-08, 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 8.26666e-08, 1.419747e-07, 3.19899e-07; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath err_ref << 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.1141397, 0.09943516, 0.09674474, 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.09980447, 0.1073116, 0.1220445, 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.1526814, 1, 1, 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1, 1, 1, 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1, 1, 1; 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fvec, fvec_ref); 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fvecp, fvecp_ref); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(err, err_ref); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Generic functor 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _Scalar, int NX=Dynamic, int NY=Dynamic> 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct Functor 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath InputsAtCompileTime = NX, 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ValuesAtCompileTime = NY 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m_inputs, m_values; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {} 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int inputs() const { return m_inputs; } 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int values() const { return m_values; } 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // you should define that in the subclass : 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// void operator() (const InputType& x, ValueType* v, JacobianType* _j=0) const; 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lmder_functor : Functor<double> 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmder_functor(void): Functor<double>(3,15) {} 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) const 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3; 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[15] = {1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i < values(); i++) 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &x, MatrixXd &fjac) const 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3, tmp4; 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i < values(); i++) 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = -1; 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = tmp1*tmp2/tmp4; 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = tmp1*tmp3/tmp4; 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmder1() 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n=3, info; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x; 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmder_functor functor; 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmder_functor> lm(functor); 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.lmder1(x); 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmder() 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m=15, n=3; 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fnorm, covfac; 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmder_functor functor; 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmder_functor> lm(functor); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return values 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = lm.fvec.blueNorm(); 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fnorm, 0.09063596); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check covariance 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath covfac = fnorm*fnorm/(m-n); 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov_ref(n,n); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov_ref << 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.0001531202, 0.002869941, -0.002656662, 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.002869941, 0.09480935, -0.09098995, 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.002656662, -0.09098995, 0.08778727; 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << fjac*covfac << std::endl; 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov; 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov = covfac*lm.fjac.topLeftCorner<n,n>(); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX( cov, cov_ref); 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: why isn't this allowed ? : 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct hybrj_functor : Functor<double> 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrj_functor(void) : Functor<double>(9,9) {} 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double temp, temp1, temp2; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = x.size(); 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==n); 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k = 0; k < n; k++) 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = (3. - 2.*x[k])*x[k]; 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp1 = 0.; 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k) temp1 = x[k-1]; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp2 = 0.; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != n-1) temp2 = x[k+1]; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[k] = temp - temp1 - 2.*temp2 + 1.; 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &x, MatrixXd &fjac) 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = x.size(); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==n); 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==n); 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k = 0; k < n; k++) 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j = 0; j < n; j++) 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(k,j) = 0.; 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(k,k) = 3.- 4.*x[k]; 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k) fjac(k,k-1) = -1.; 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != n-1) fjac(k,k+1) = -2.; 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrj1() 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=9; 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrj_functor functor; 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrj_functor> solver(functor); 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.hybrj1(x); 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 11); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.njev, 1); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// check x 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.5706545, -0.6816283, -0.7017325, 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.7042129, -0.701369, -0.6918656, 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.665792, -0.5960342, -0.4164121; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrj() 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=9; 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrj_functor functor; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrj_functor> solver(functor); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.diag.setConstant(n, 1.); 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.useExternalScaling = true; 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.solve(x); 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 11); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.njev, 1); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// check x 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.5706545, -0.6816283, -0.7017325, 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.7042129, -0.701369, -0.6918656, 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.665792, -0.5960342, -0.4164121; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct hybrd_functor : Functor<double> 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrd_functor(void) : Functor<double>(9,9) {} 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) const 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double temp, temp1, temp2; 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = x.size(); 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==n); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k=0; k < n; k++) 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = (3. - 2.*x[k])*x[k]; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp1 = 0.; 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k) temp1 = x[k-1]; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp2 = 0.; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != n-1) temp2 = x[k+1]; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[k] = temp - temp1 - 2.*temp2 + 1.; 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrd1() 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n=9, info; 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough solution. */ 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrd_functor functor; 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrd_functor> solver(functor); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.hybrd1(x); 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 20); 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrd() 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=9; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x; 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrd_functor functor; 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrd_functor> solver(functor); 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.parameters.nb_of_subdiagonals = 1; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.parameters.nb_of_superdiagonals = 1; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.diag.setConstant(n, 1.); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.useExternalScaling = true; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.solveNumericalDiff(x); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 14); 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.5706545, -0.6816283, -0.7017325, 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.7042129, -0.701369, -0.6918656, 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.665792, -0.5960342, -0.4164121; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lmstr_functor : Functor<double> 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmstr_functor(void) : Functor<double>(3,15) {} 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* subroutine fcn for lmstr1 example. */ 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3; 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[15]={1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(15==fvec.size()); 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(3==x.size()); 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i=0; i<15; i++) 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb) 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(x.size()==3); 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(jac_row.size()==x.size()); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3, tmp4; 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i = rownb-2; 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac_row[0] = -1; 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac_row[1] = tmp1*tmp2/tmp4; 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac_row[2] = tmp1*tmp3/tmp4; 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmstr1() 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmstr_functor functor; 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmstr_functor> lm(functor); 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.lmstr1(x); 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695 ; 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmstr() 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fnorm; 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmstr_functor functor; 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmstr_functor> lm(functor); 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimizeOptimumStorage(x); 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return values 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = lm.fvec.blueNorm(); 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fnorm, 0.09063596); 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lmdif_functor : Functor<double> 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmdif_functor(void) : Functor<double>(3,15) {} 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) const 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1,tmp2,tmp3; 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[15]={1.4e-1,1.8e-1,2.2e-1,2.5e-1,2.9e-1,3.2e-1,3.5e-1,3.9e-1, 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(x.size()==3); 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==15); 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i=0; i<15; i++) 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 15 - i; 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = tmp1; 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i >= 8) tmp3 = tmp2; 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmdif1() 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n), fvec(15); 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmdif_functor functor; 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex nfev; 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(nfev, 26); 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath functor(x, fvec); 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.0824106, 1.1330366, 2.3436947; 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmdif() 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m=15, n=3; 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fnorm, covfac; 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmdif_functor functor; 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff<lmdif_functor> numDiff(functor); 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return values 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 26); 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = lm.fvec.blueNorm(); 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fnorm, 0.09063596); 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check covariance 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath covfac = fnorm*fnorm/(m-n); 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov_ref(n,n); 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov_ref << 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.0001531202, 0.002869942, -0.002656662, 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.002869942, 0.09480937, -0.09098997, 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.002656662, -0.09098997, 0.08778729; 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << fjac*covfac << std::endl; 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov; 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov = covfac*lm.fjac.topLeftCorner<n,n>(); 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX( cov, cov_ref); 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: why isn't this allowed ? : 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct chwirut2_functor : Functor<double> 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chwirut2_functor(void) : Functor<double>(3,54) {} 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_x[54]; 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_y[54]; 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==54); 642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(i=0; i<54; i++) { 643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = m_x[i]; 644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==54); 652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<54; i++) { 654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = m_x[i]; 655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double factor = 1./(b[1]+b[2]*x); 656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(-b[0]*x); 657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = -x*e*factor; 658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = -e*factor*factor; 659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = -x*e*factor*factor; 660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double chwirut2_functor::m_x[54] = { 0.500E0, 1.000E0, 1.750E0, 3.750E0, 5.750E0, 0.875E0, 2.250E0, 3.250E0, 5.250E0, 0.750E0, 1.750E0, 2.750E0, 4.750E0, 0.625E0, 1.250E0, 2.250E0, 4.250E0, .500E0, 3.000E0, .750E0, 3.000E0, 1.500E0, 6.000E0, 3.000E0, 6.000E0, 1.500E0, 3.000E0, .500E0, 2.000E0, 4.000E0, .750E0, 2.000E0, 5.000E0, .750E0, 2.250E0, 3.750E0, 5.750E0, 3.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .750E0, 2.500E0, 4.000E0, .500E0, 6.000E0, 3.000E0, .500E0, 2.750E0, .500E0, 1.750E0}; 665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double chwirut2_functor::m_y[54] = { 92.9000E0 ,57.1000E0 ,31.0500E0 ,11.5875E0 ,8.0250E0 ,63.6000E0 ,21.4000E0 ,14.2500E0 ,8.4750E0 ,63.8000E0 ,26.8000E0 ,16.4625E0 ,7.1250E0 ,67.3000E0 ,41.0000E0 ,21.1500E0 ,8.1750E0 ,81.5000E0 ,13.1200E0 ,59.9000E0 ,14.6200E0 ,32.9000E0 ,5.4400E0 ,12.5600E0 ,5.4400E0 ,32.0000E0 ,13.9500E0 ,75.8000E0 ,20.0000E0 ,10.4200E0 ,59.5000E0 ,21.6700E0 ,8.5500E0 ,62.0000E0 ,20.2000E0 ,7.7600E0 ,3.7500E0 ,11.8100E0 ,54.7000E0 ,23.7000E0 ,11.5500E0 ,61.3000E0 ,17.7000E0 ,8.7400E0 ,59.2000E0 ,16.3000E0 ,8.6200E0 ,81.0000E0 ,4.8700E0 ,14.6200E0 ,81.7000E0 ,17.1700E0 ,81.3000E0 ,28.9000E0 }; 666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistChwirut2(void) 669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 676c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 677c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 678c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.1, 0.01, 0.02; 679c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 680c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chwirut2_functor functor; 681c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<chwirut2_functor> lm(functor); 682c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 683c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 684c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 10); 687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 690c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 691c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 692c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 693c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 694c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 695c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 696c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 697c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 698c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.15, 0.008, 0.010; 699c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 700c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 701c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 702c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 703c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 704c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 705c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 706c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 707c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 7); 708c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 6); 709c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 710c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 711c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 712c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 713c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 714c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 715c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 716c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 717c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 718c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct misra1a_functor : Functor<double> 719c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 720c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1a_functor(void) : Functor<double>(2,14) {} 721c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_x[14]; 722c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_y[14]; 723c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 724c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 725c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 726c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==14); 727c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 728c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 729c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 730c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 731c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 732c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 733c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 734c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 735c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==14); 736c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==2); 737c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 738c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 739c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 740c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 741c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 742c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 743c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 744c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double misra1a_functor::m_x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; 745c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double misra1a_functor::m_y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; 746c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 747c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 748c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMisra1a(void) 749c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 750c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=2; 751c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 752c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 753c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 754c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 755c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 756c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 757c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 758c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 500., 0.0001; 759c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 760c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1a_functor functor; 761c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<misra1a_functor> lm(functor); 762c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 763c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 764c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 765c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 766c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 19); 767c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 15); 768c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 769c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 770c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 771c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 772c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 773c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 774c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 775c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 776c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 777c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 250., 0.0005; 778c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 779c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 780c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 781c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 782c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 783c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 5); 784c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 4); 785c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 786c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 787c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 788c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 789c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 790c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 791c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 792c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct hahn1_functor : Functor<double> 793c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 794c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hahn1_functor(void) : Functor<double>(7,236) {} 795c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_x[236]; 796c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 797c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 798c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0 , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0 , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0 , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 , 16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0 , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0 , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 , 13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 }; 799c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 800c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 801c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 802c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 803c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==236); 804c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<236; i++) { 805c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=m_x[i], xx=x*x, xxx=xx*x; 806c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - m_y[i]; 807c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 808c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 809c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 810c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 811c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 812c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 813c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 814c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==236); 815c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==7); 816c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<236; i++) { 817c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=m_x[i], xx=x*x, xxx=xx*x; 818c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 819c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.*fact; 820c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = x*fact; 821c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = xx*fact; 822c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = xxx*fact; 823c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 824c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = x*fact; 825c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,5) = xx*fact; 826c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,6) = xxx*fact; 827c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 828c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 829c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 830c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 831c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , 282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , 141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0}; 832c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 833c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 834c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistHahn1(void) 835c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 836c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=7; 837c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 838c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 839c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 840c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 841c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 842c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 843c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 844c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 845c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 846c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hahn1_functor functor; 847c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<hahn1_functor> lm(functor); 848c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 849c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 850c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 851c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 852c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 11); 853c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 10); 854c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 855c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 856c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 857c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 858c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 859c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 860c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 861c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 862c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 863c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 864c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 865c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 866c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 867c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 868c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 869c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 870c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 871c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 872c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 873c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 874c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 11); 875c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 10); 876c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 877c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 878c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 879c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 880c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 881c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 882c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 883c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 884c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 885c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 886c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 887c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 888c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 889c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct misra1d_functor : Functor<double> 890c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 891c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1d_functor(void) : Functor<double>(2,14) {} 892c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[14]; 893c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[14]; 894c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 895c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 896c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 897c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==14); 898c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 899c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 900c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 901c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 902c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 903c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 904c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 905c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 906c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==14); 907c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==2); 908c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 909c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double den = 1.+b[1]*x[i]; 910c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = b[1]*x[i] / den; 911c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 912c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 913c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 914c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 915c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 916c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double misra1d_functor::x[14] = { 77.6E0, 114.9E0, 141.1E0, 190.8E0, 239.9E0, 289.0E0, 332.8E0, 378.4E0, 434.8E0, 477.3E0, 536.8E0, 593.1E0, 689.1E0, 760.0E0}; 917c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double misra1d_functor::y[14] = { 10.07E0, 14.73E0, 17.94E0, 23.93E0, 29.61E0, 35.18E0, 40.02E0, 44.82E0, 50.76E0, 55.05E0, 61.01E0, 66.40E0, 75.47E0, 81.78E0}; 918c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 919c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 920c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMisra1d(void) 921c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 922c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=2; 923c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 924c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 925c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 926c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 927c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 928c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 929c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 930c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 500., 0.0001; 931c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 932c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1d_functor functor; 933c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<misra1d_functor> lm(functor); 934c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 935c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 936c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 937c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 3); 938c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 9); 939c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 7); 940c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 941c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 942c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 943c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 944c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 945c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 946c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 947c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 948c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 949c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 450., 0.0003; 950c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 951c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 952c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 953c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 954c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 955c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 4); 956c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 3); 957c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 958c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 959c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 960c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 961c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 962c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 963c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 964c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 965c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lanczos1_functor : Functor<double> 966c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 967c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lanczos1_functor(void) : Functor<double>(6,24) {} 968c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[24]; 969c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[24]; 970c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 971c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 972c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==6); 973c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==24); 974c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<24; i++) 975c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*exp(-b[1]*x[i]) + b[2]*exp(-b[3]*x[i]) + b[4]*exp(-b[5]*x[i]) - y[i]; 976c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 977c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 978c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 979c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 980c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==6); 981c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==24); 982c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==6); 983c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<24; i++) { 984c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = exp(-b[1]*x[i]); 985c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 986c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = exp(-b[3]*x[i]); 987c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 988c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = exp(-b[5]*x[i]); 989c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 990c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 991c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 992c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 993c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 994c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double lanczos1_functor::x[24] = { 0.000000000000E+00, 5.000000000000E-02, 1.000000000000E-01, 1.500000000000E-01, 2.000000000000E-01, 2.500000000000E-01, 3.000000000000E-01, 3.500000000000E-01, 4.000000000000E-01, 4.500000000000E-01, 5.000000000000E-01, 5.500000000000E-01, 6.000000000000E-01, 6.500000000000E-01, 7.000000000000E-01, 7.500000000000E-01, 8.000000000000E-01, 8.500000000000E-01, 9.000000000000E-01, 9.500000000000E-01, 1.000000000000E+00, 1.050000000000E+00, 1.100000000000E+00, 1.150000000000E+00 }; 995c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double lanczos1_functor::y[24] = { 2.513400000000E+00 ,2.044333373291E+00 ,1.668404436564E+00 ,1.366418021208E+00 ,1.123232487372E+00 ,9.268897180037E-01 ,7.679338563728E-01 ,6.388775523106E-01 ,5.337835317402E-01 ,4.479363617347E-01 ,3.775847884350E-01 ,3.197393199326E-01 ,2.720130773746E-01 ,2.324965529032E-01 ,1.996589546065E-01 ,1.722704126914E-01 ,1.493405660168E-01 ,1.300700206922E-01 ,1.138119324644E-01 ,1.000415587559E-01 ,8.833209084540E-02 ,7.833544019350E-02 ,6.976693743449E-02 ,6.239312536719E-02 }; 996c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 997c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 998c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistLanczos1(void) 999c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1000c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=6; 1001c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1002c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1003c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1004c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1005c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1006c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1007c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1008c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 1009c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1010c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lanczos1_functor functor; 1011c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lanczos1_functor> lm(functor); 1012c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1013c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1014c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1015c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1016c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 79); 1017c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 72); 1018c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1019c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 1020c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1021c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1022c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1023c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1024c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1025c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1026c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1027c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1028c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1029c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1030c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1031c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 1032c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1033c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1034c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1035c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1036c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1037c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 9); 1038c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 1039c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1040c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 1041c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1042c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1043c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1044c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1045c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1046c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1047c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1048c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1049c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1050c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1051c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct rat42_functor : Functor<double> 1052c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1053c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat42_functor(void) : Functor<double>(3,9) {} 1054c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[9]; 1055c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[9]; 1056c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1057c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1058c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1059c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==9); 1060c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<9; i++) { 1061c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 1062c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1063c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1064c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1065c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1066c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1067c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1068c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1069c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==9); 1070c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1071c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<9; i++) { 1072c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(b[1]-b[2]*x[i]); 1073c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1./(1.+e); 1074c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 1075c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 1076c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1077c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1078c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1079c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1080c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 1081c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 1082c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1083c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 1084c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistRat42(void) 1085c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1086c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1087c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1088c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1089c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1090c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1091c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1092c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1093c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1094c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 100., 1., 0.1; 1095c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1096c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat42_functor functor; 1097c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<rat42_functor> lm(functor); 1098c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1099c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 10); 1103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 1104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 75., 2.5, 0.07; 1115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 1121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 1122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct MGH10_functor : Functor<double> 1131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH10_functor(void) : Functor<double>(3,16) {} 1133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[16]; 1134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[16]; 1135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==16); 1139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<16; i++) 1140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; 1141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==16); 1147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<16; i++) { 1149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double factor = 1./(x[i]+b[2]); 1150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(b[1]*factor); 1151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = e; 1152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*factor*e; 1153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = -b[1]*b[0]*factor*factor*e; 1154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double MGH10_functor::x[16] = { 5.000000E+01, 5.500000E+01, 6.000000E+01, 6.500000E+01, 7.000000E+01, 7.500000E+01, 8.000000E+01, 8.500000E+01, 9.000000E+01, 9.500000E+01, 1.000000E+02, 1.050000E+02, 1.100000E+02, 1.150000E+02, 1.200000E+02, 1.250000E+02 }; 1159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double MGH10_functor::y[16] = { 3.478000E+04, 2.861000E+04, 2.365000E+04, 1.963000E+04, 1.637000E+04, 1.372000E+04, 1.154000E+04, 9.744000E+03, 8.261000E+03, 7.030000E+03, 6.005000E+03, 5.147000E+03, 4.427000E+03, 3.820000E+03, 3.307000E+03, 2.872000E+03 }; 1160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml 1162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMGH10(void) 1163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 2., 400000., 25000.; 1173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH10_functor functor; 1175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<MGH10_functor> lm(functor); 1176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 284 ); 1181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 249 ); 1182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.02, 4000., 250.; 1193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 3); 1198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 126); 1199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 116); 1200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct BoxBOD_functor : Functor<double> 1210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BoxBOD_functor(void) : Functor<double>(2,6) {} 1212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[6]; 1213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 1216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 1217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==6); 1218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<6; i++) 1219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 1220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 1225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==6); 1226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==2); 1227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<6; i++) { 1228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(-b[1]*x[i]); 1229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.-e; 1230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*x[i]*e; 1231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 1236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 1238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistBoxBOD(void) 1239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=2; 1241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1., 1.; 1249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BoxBOD_functor functor; 1251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<BoxBOD_functor> lm(functor); 1252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.factor = 10.; 1255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 31); 1260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 25); 1261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 100., 0.75; 1271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = NumTraits<double>::epsilon(); 1274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = NumTraits<double>::epsilon(); 1275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 15 ); 1280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 14 ); 1281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct MGH17_functor : Functor<double> 1289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH17_functor(void) : Functor<double>(5,33) {} 1291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[33]; 1292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[33]; 1293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==5); 1296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==33); 1297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<33; i++) 1298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 1299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==5); 1304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==33); 1305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==5); 1306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<33; i++) { 1307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.; 1308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = exp(-b[3]*x[i]); 1309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = exp(-b[4]*x[i]); 1310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 1311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 1312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double MGH17_functor::x[33] = { 0.000000E+00, 1.000000E+01, 2.000000E+01, 3.000000E+01, 4.000000E+01, 5.000000E+01, 6.000000E+01, 7.000000E+01, 8.000000E+01, 9.000000E+01, 1.000000E+02, 1.100000E+02, 1.200000E+02, 1.300000E+02, 1.400000E+02, 1.500000E+02, 1.600000E+02, 1.700000E+02, 1.800000E+02, 1.900000E+02, 2.000000E+02, 2.100000E+02, 2.200000E+02, 2.300000E+02, 2.400000E+02, 2.500000E+02, 2.600000E+02, 2.700000E+02, 2.800000E+02, 2.900000E+02, 3.000000E+02, 3.100000E+02, 3.200000E+02 }; 1317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double MGH17_functor::y[33] = { 8.440000E-01, 9.080000E-01, 9.320000E-01, 9.360000E-01, 9.250000E-01, 9.080000E-01, 8.810000E-01, 8.500000E-01, 8.180000E-01, 7.840000E-01, 7.510000E-01, 7.180000E-01, 6.850000E-01, 6.580000E-01, 6.280000E-01, 6.030000E-01, 5.800000E-01, 5.580000E-01, 5.380000E-01, 5.220000E-01, 5.060000E-01, 4.900000E-01, 4.780000E-01, 4.670000E-01, 4.570000E-01, 4.480000E-01, 4.380000E-01, 4.310000E-01, 4.240000E-01, 4.200000E-01, 4.140000E-01, 4.110000E-01, 4.060000E-01 }; 1318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 1320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMGH17(void) 1321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=5; 1323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 50., 150., -100., 1., 2.; 1331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH17_functor functor; 1333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<MGH17_functor> lm(functor); 1334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = NumTraits<double>::epsilon(); 1335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = NumTraits<double>::epsilon(); 1336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.maxfev = 1000; 1337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 602 ); 1342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 545 ); 1343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 1356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 18); 1363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 15); 1364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct MGH09_functor : Functor<double> 1375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH09_functor(void) : Functor<double>(4,11) {} 1377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double _x[11]; 1378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[11]; 1379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==11); 1383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<11; i++) { 1384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = _x[i], xx=x*x; 1385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 1386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==11); 1393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==4); 1394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<11; i++) { 1395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = _x[i], xx=x*x; 1396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double factor = 1./(xx+x*b[2]+b[3]); 1397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = (xx+x*b[1]) * factor; 1398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*x* factor; 1399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 1400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 1401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double MGH09_functor::_x[11] = { 4., 2., 1., 5.E-1 , 2.5E-01, 1.670000E-01, 1.250000E-01, 1.E-01, 8.330000E-02, 7.140000E-02, 6.250000E-02 }; 1406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double MGH09_functor::y[11] = { 1.957000E-01, 1.947000E-01, 1.735000E-01, 1.600000E-01, 8.440000E-02, 6.270000E-02, 4.560000E-02, 3.420000E-02, 3.230000E-02, 2.350000E-02, 2.460000E-02 }; 1407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 1409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMGH09(void) 1410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=4; 1412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 25., 39, 41.5, 39.; 1420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH09_functor functor; 1422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<MGH09_functor> lm(functor); 1423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.maxfev = 1000; 1424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 490 ); 1429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 376 ); 1430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 1434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 1435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 1436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 1437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.25, 0.39, 0.415, 0.39; 1442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 18); 1449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 16); 1450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 1454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 1455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 1456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 1457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct Bennett5_functor : Functor<double> 1462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Bennett5_functor(void) : Functor<double>(3,154) {} 1464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[154]; 1465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[154]; 1466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==154); 1470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<154; i++) 1471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 1472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==154); 1478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<154; i++) { 1480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = pow(b[1]+x[i],-1./b[2]); 1481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = e; 1482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 1483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 1484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, 11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 }; 1489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 ,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 ,-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 1490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 1492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistBennett5(void) 1493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< -2000., 50., 0.8; 1503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Bennett5_functor functor; 1505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<Bennett5_functor> lm(functor); 1506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.maxfev = 1000; 1507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 758); 1512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 744); 1513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 1517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 1518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 1519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< -1500., 45., 0.85; 1523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 203); 1530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 192); 1531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 1535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 1536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 1537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct thurber_functor : Functor<double> 1540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath thurber_functor(void) : Functor<double>(7,37) {} 1542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double _x[37]; 1543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double _y[37]; 1544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 1547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 1548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==37); 1549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<37; i++) { 1550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=_x[i], xx=x*x, xxx=xx*x; 1551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) / (1.+b[4]*x+b[5]*xx+b[6]*xxx) - _y[i]; 1552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 1558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==37); 1559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==7); 1560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<37; i++) { 1561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=_x[i], xx=x*x, xxx=xx*x; 1562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 1563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.*fact; 1564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = x*fact; 1565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = xx*fact; 1566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = xxx*fact; 1567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 1568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = x*fact; 1569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,5) = xx*fact; 1570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,6) = xxx*fact; 1571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double thurber_functor::_x[37] = { -3.067E0, -2.981E0, -2.921E0, -2.912E0, -2.840E0, -2.797E0, -2.702E0, -2.699E0, -2.633E0, -2.481E0, -2.363E0, -2.322E0, -1.501E0, -1.460E0, -1.274E0, -1.212E0, -1.100E0, -1.046E0, -0.915E0, -0.714E0, -0.566E0, -0.545E0, -0.400E0, -0.309E0, -0.109E0, -0.103E0, 0.010E0, 0.119E0, 0.377E0, 0.790E0, 0.963E0, 1.006E0, 1.115E0, 1.572E0, 1.841E0, 2.047E0, 2.200E0 }; 1576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double thurber_functor::_y[37] = { 80.574E0, 84.248E0, 87.264E0, 87.195E0, 89.076E0, 89.608E0, 89.868E0, 90.101E0, 92.405E0, 95.854E0, 100.696E0, 101.060E0, 401.672E0, 390.724E0, 567.534E0, 635.316E0, 733.054E0, 759.087E0, 894.206E0, 990.785E0, 1090.109E0, 1080.914E0, 1122.643E0, 1178.351E0, 1260.531E0, 1273.514E0, 1288.339E0, 1327.543E0, 1353.863E0, 1414.509E0, 1425.208E0, 1421.384E0, 1442.962E0, 1464.350E0, 1468.705E0, 1447.894E0, 1457.628E0}; 1577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 1579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistThurber(void) 1580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=7; 1582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 1590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath thurber_functor functor; 1592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<thurber_functor> lm(functor); 1593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 39); 1600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 36); 1601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 1616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 29); 1625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 28); 1626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct rat43_functor : Functor<double> 1639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat43_functor(void) : Functor<double>(4,15) {} 1641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[15]; 1642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[15]; 1643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==15); 1647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<15; i++) 1648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 1649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==15); 1655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==4); 1656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<15; i++) { 1657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(b[1]-b[2]*x[i]); 1658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double power = -1./b[3]; 1659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = pow(1.+e, power); 1660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 1661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 1662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 1663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 1668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat43_functor::y[15] = { 16.08, 33.83, 65.80, 97.20, 191.55, 326.20, 386.87, 520.53, 590.03, 651.92, 724.93, 699.56, 689.96, 637.56, 717.41 }; 1669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 1671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistRat43(void) 1672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=4; 1674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1676c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1677c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1678c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1679c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1680c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1681c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 100., 10., 1., 1.; 1682c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1683c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat43_functor functor; 1684c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<rat43_functor> lm(functor); 1685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1690c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1691c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 27); 1692c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 20); 1693c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1694c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1695c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1696c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1697c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1698c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1699c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1700c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1701c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1702c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1703c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1704c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 700., 5., 0.75, 1.3; 1705c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1706c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1707c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon(); 1708c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon(); 1709c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1710c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1711c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1712c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1713c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 9); 1714c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 1715c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1716c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1717c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1718c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1719c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1720c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1721c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1722c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1723c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1724c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1725c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1726c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct eckerle4_functor : Functor<double> 1727c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1728c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eckerle4_functor(void) : Functor<double>(3,35) {} 1729c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[35]; 1730c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[35]; 1731c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1732c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1733c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1734c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==35); 1735c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<35; i++) 1736c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 1737c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1738c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1739c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1740c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1741c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1742c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==35); 1743c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1744c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<35; i++) { 1745c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double b12 = b[1]*b[1]; 1746c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 1747c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = e / b[1]; 1748c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 1749c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 1750c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1751c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1752c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1753c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1754c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double eckerle4_functor::x[35] = { 400.0, 405.0, 410.0, 415.0, 420.0, 425.0, 430.0, 435.0, 436.5, 438.0, 439.5, 441.0, 442.5, 444.0, 445.5, 447.0, 448.5, 450.0, 451.5, 453.0, 454.5, 456.0, 457.5, 459.0, 460.5, 462.0, 463.5, 465.0, 470.0, 475.0, 480.0, 485.0, 490.0, 495.0, 500.0}; 1755c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double eckerle4_functor::y[35] = { 0.0001575, 0.0001699, 0.0002350, 0.0003102, 0.0004917, 0.0008710, 0.0017418, 0.0046400, 0.0065895, 0.0097302, 0.0149002, 0.0237310, 0.0401683, 0.0712559, 0.1264458, 0.2073413, 0.2902366, 0.3445623, 0.3698049, 0.3668534, 0.3106727, 0.2078154, 0.1164354, 0.0616764, 0.0337200, 0.0194023, 0.0117831, 0.0074357, 0.0022732, 0.0008800, 0.0004579, 0.0002345, 0.0001586, 0.0001143, 0.0000710 }; 1756c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1757c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 1758c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistEckerle4(void) 1759c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1760c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1761c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1762c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1763c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1764c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1765c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1766c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1767c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1768c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1., 10., 500.; 1769c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1770c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eckerle4_functor functor; 1771c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<eckerle4_functor> lm(functor); 1772c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1773c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1774c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1775c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1776c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 18); 1777c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 15); 1778c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1779c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1780c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1781c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.5543827178); 1782c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 4.0888321754); 1783c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1784c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1785c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1786c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1787c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1788c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1.5, 5., 450.; 1789c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1790c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1791c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1792c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1793c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1794c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 7); 1795c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 6); 1796c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1797c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1798c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1799c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.5543827178); 1800c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 4.0888321754); 1801c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1802c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1803c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1804c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_NonLinearOptimization() 1805c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1806c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Tests using the examples provided by (c)minpack 1807c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_1*/(testChkder()); 1808c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_1*/(testLmder1()); 1809c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_1*/(testLmder()); 1810c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrj1()); 1811c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrj()); 1812c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrd1()); 1813c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrd()); 1814c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmstr1()); 1815c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmstr()); 1816c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmdif1()); 1817c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmdif()); 1818c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1819c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // NIST tests, level of difficulty = "Lower" 1820c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_4*/(testNistMisra1a()); 1821c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_4*/(testNistChwirut2()); 1822c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1823c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // NIST tests, level of difficulty = "Average" 1824c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_5*/(testNistHahn1()); 1825c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_6*/(testNistMisra1d()); 1826c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_7*/(testNistMGH17()); 1827c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_8*/(testNistLanczos1()); 1828c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1829c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // NIST tests, level of difficulty = "Higher" 1830c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_9*/(testNistRat42()); 1831c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_10*/(testNistMGH10()); 1832c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_11*/(testNistBoxBOD()); 1833c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_12*/(testNistMGH09()); 1834c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_13*/(testNistBennett5()); 1835c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_14*/(testNistThurber()); 1836c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_15*/(testNistRat43()); 1837c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_16*/(testNistEckerle4()); 1838c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1839c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1840c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* 1841c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Can be useful for debugging... 1842c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev : %d, %d\n", info, lm.nfev); 1843c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev); 1844c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev : %d, %d\n", info, solver.nfev); 1845c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[0] : %.32g\n", x[0]); 1846c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[1] : %.32g\n", x[1]); 1847c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[2] : %.32g\n", x[2]); 1848c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[3] : %.32g\n", x[3]); 1849c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm()); 1850c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm()); 1851c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1852c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev); 1853c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm()); 1854c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x << std::endl; 1855c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout.precision(9); 1856c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[0] << std::endl; 1857c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[1] << std::endl; 1858c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[2] << std::endl; 1859c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[3] << std::endl; 1860c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath*/ 1861c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1862