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 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezusing std::sqrt; 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathint fcn_chkder(const VectorXd &x, VectorXd &fvec, MatrixXd &fjac, int iflag) 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* subroutine fcn for chkder example. */ 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(15 == fvec.size()); 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(3 == x.size()); 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3, tmp4; 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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, 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iflag == 0) 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (iflag != 2) 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i=0; i<15; i++) { 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16-i-1; 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = tmp1; 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i >= 8) tmp3 = tmp2; 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else { 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i = 0; i < 15; i++) { 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16-i-1; 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* error introduced into next statement for illustration. */ 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* corrected statement should read tmp3 = tmp1 . */ 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = tmp2; 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i >= 8) tmp3 = tmp2; 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4=tmp4*tmp4; 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = -1.; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = tmp1*tmp2/tmp4; 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = tmp1*tmp3/tmp4; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testChkder() 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m=15, n=3; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n), fvec(m), xp, fvecp(m), err; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd fjac(m,n); 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXi ipvt; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following values should be suitable for */ 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* checking the jacobian matrix. */ 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x << 9.2e-1, 1.3e-1, 5.4e-1; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::chkder(x, fvec, fjac, xp, fvecp, 1, err); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fcn_chkder(x, fvec, fjac, 1); 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fcn_chkder(x, fvec, fjac, 2); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fcn_chkder(xp, fvecp, fjac, 1); 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::chkder(x, fvec, fjac, xp, fvecp, 2, err); 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvecp -= fvec; 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check those 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd fvec_ref(m), fvecp_ref(m), err_ref(m); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec_ref << 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -1.181606, -1.429655, -1.606344, 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -1.745269, -1.840654, -1.921586, 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -1.984141, -2.022537, -2.468977, 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -2.827562, -3.473582, -4.437612, 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -6.047662, -9.267761, -18.91806; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvecp_ref << 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -7.724666e-09, -3.432406e-09, -2.034843e-10, 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2.313685e-09, 4.331078e-09, 5.984096e-09, 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 7.363281e-09, 8.53147e-09, 1.488591e-08, 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2.33585e-08, 3.522012e-08, 5.301255e-08, 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 8.26666e-08, 1.419747e-07, 3.19899e-07; 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath err_ref << 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.1141397, 0.09943516, 0.09674474, 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.09980447, 0.1073116, 0.1220445, 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.1526814, 1, 1, 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1, 1, 1, 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1, 1, 1; 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fvec, fvec_ref); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fvecp, fvecp_ref); 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(err, err_ref); 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Generic functor 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _Scalar, int NX=Dynamic, int NY=Dynamic> 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct Functor 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _Scalar Scalar; 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath InputsAtCompileTime = NX, 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ValuesAtCompileTime = NY 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,InputsAtCompileTime,1> InputType; 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,ValuesAtCompileTime,1> ValueType; 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m_inputs, m_values; 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {} 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {} 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int inputs() const { return m_inputs; } 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int values() const { return m_values; } 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // you should define that in the subclass : 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// void operator() (const InputType& x, ValueType* v, JacobianType* _j=0) const; 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lmder_functor : Functor<double> 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmder_functor(void): Functor<double>(3,15) {} 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) const 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3; 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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, 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i < values(); i++) 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &x, MatrixXd &fjac) const 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3, tmp4; 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i = 0; i < values(); i++) 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = -1; 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = tmp1*tmp2/tmp4; 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = tmp1*tmp3/tmp4; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmder1() 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n=3, info; 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x; 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmder_functor functor; 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmder_functor> lm(functor); 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.lmder1(x); 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmder() 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m=15, n=3; 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fnorm, covfac; 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmder_functor functor; 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmder_functor> lm(functor); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return values 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = lm.fvec.blueNorm(); 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fnorm, 0.09063596); 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check covariance 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath covfac = fnorm*fnorm/(m-n); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov_ref(n,n); 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov_ref << 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.0001531202, 0.002869941, -0.002656662, 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.002869941, 0.09480935, -0.09098995, 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.002656662, -0.09098995, 0.08778727; 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << fjac*covfac << std::endl; 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov; 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov = covfac*lm.fjac.topLeftCorner<n,n>(); 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX( cov, cov_ref); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: why isn't this allowed ? : 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct hybrj_functor : Functor<double> 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrj_functor(void) : Functor<double>(9,9) {} 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double temp, temp1, temp2; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = x.size(); 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==n); 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k = 0; k < n; k++) 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = (3. - 2.*x[k])*x[k]; 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp1 = 0.; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k) temp1 = x[k-1]; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp2 = 0.; 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != n-1) temp2 = x[k+1]; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[k] = temp - temp1 - 2.*temp2 + 1.; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &x, MatrixXd &fjac) 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = x.size(); 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==n); 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==n); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k = 0; k < n; k++) 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int j = 0; j < n; j++) 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(k,j) = 0.; 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(k,k) = 3.- 4.*x[k]; 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k) fjac(k,k-1) = -1.; 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != n-1) fjac(k,k+1) = -2.; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrj1() 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=9; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrj_functor functor; 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrj_functor> solver(functor); 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.hybrj1(x); 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 11); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.njev, 1); 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// check x 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.5706545, -0.6816283, -0.7017325, 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.7042129, -0.701369, -0.6918656, 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.665792, -0.5960342, -0.4164121; 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrj() 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=9; 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrj_functor functor; 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrj_functor> solver(functor); 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.diag.setConstant(n, 1.); 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.useExternalScaling = true; 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.solve(x); 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 11); 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.njev, 1); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// check x 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.5706545, -0.6816283, -0.7017325, 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.7042129, -0.701369, -0.6918656, 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.665792, -0.5960342, -0.4164121; 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct hybrd_functor : Functor<double> 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrd_functor(void) : Functor<double>(9,9) {} 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) const 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double temp, temp1, temp2; 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n = x.size(); 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==n); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int k=0; k < n; k++) 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = (3. - 2.*x[k])*x[k]; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp1 = 0.; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k) temp1 = x[k-1]; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp2 = 0.; 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != n-1) temp2 = x[k+1]; 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[k] = temp - temp1 - 2.*temp2 + 1.; 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrd1() 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n=9, info; 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough solution. */ 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrd_functor functor; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrd_functor> solver(functor); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.hybrd1(x); 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 20); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << -0.5706545, -0.6816283, -0.7017325, -0.7042129, -0.701369, -0.6918656, -0.665792, -0.5960342, -0.4164121; 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testHybrd() 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=9; 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x; 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, -1.); 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hybrd_functor functor; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HybridNonLinearSolver<hybrd_functor> solver(functor); 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.parameters.nb_of_subdiagonals = 1; 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.parameters.nb_of_superdiagonals = 1; 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.diag.setConstant(n, 1.); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solver.useExternalScaling = true; 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = solver.solveNumericalDiff(x); 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(solver.nfev, 14); 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(solver.fvec.blueNorm(), 1.192636e-08); 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.5706545, -0.6816283, -0.7017325, 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.7042129, -0.701369, -0.6918656, 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.665792, -0.5960342, -0.4164121; 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lmstr_functor : Functor<double> 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmstr_functor(void) : Functor<double>(3,15) {} 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* subroutine fcn for lmstr1 example. */ 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3; 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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, 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(15==fvec.size()); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(3==x.size()); 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (int i=0; i<15; i++) 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &x, VectorXd &jac_row, VectorXd::Index rownb) 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(x.size()==3); 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(jac_row.size()==x.size()); 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1, tmp2, tmp3, tmp4; 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i = rownb-2; 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 16 - i - 1; 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = (i>=8)? tmp2 : tmp1; 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac_row[0] = -1; 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac_row[1] = tmp1*tmp2/tmp4; 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath jac_row[2] = tmp1*tmp3/tmp4; 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmstr1() 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmstr_functor functor; 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmstr_functor> lm(functor); 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.lmstr1(x); 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.blueNorm(), 0.09063596); 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695 ; 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmstr() 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fnorm; 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmstr_functor functor; 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lmstr_functor> lm(functor); 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimizeOptimumStorage(x); 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return values 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = lm.fvec.blueNorm(); 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fnorm, 0.09063596); 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lmdif_functor : Functor<double> 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmdif_functor(void) : Functor<double>(3,15) {} 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &x, VectorXd &fvec) const 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double tmp1,tmp2,tmp3; 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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, 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(x.size()==3); 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==15); 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (i=0; i<15; i++) 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp1 = i+1; 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp2 = 15 - i; 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp3 = tmp1; 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i >= 8) tmp3 = tmp2; 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmdif1() 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n), fvec(15); 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmdif_functor functor; 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath DenseIndex nfev; 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = LevenbergMarquardt<lmdif_functor>::lmdif1(functor, x, &nfev); 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(nfev, 26); 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath functor(x, fvec); 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.0824106, 1.1330366, 2.3436947; 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testLmdif() 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int m=15, n=3; 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fnorm, covfac; 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* the following starting values provide a rough fit. */ 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x.setConstant(n, 1.); 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lmdif_functor functor; 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath NumericalDiff<lmdif_functor> numDiff(functor); 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<NumericalDiff<lmdif_functor> > lm(numDiff); 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return values 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 26); 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fnorm = lm.fvec.blueNorm(); 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(fnorm, 0.09063596); 608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x_ref(n); 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x_ref << 0.08241058, 1.133037, 2.343695; 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x, x_ref); 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check covariance 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath covfac = fnorm*fnorm/(m-n); 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::covar(lm.fjac, lm.permutation.indices()); // TODO : move this as a function of lm 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov_ref(n,n); 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov_ref << 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.0001531202, 0.002869942, -0.002656662, 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 0.002869942, 0.09480937, -0.09098997, 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -0.002656662, -0.09098997, 0.08778729; 623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// std::cout << fjac*covfac << std::endl; 625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixXd cov; 627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cov = covfac*lm.fjac.topLeftCorner<n,n>(); 628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX( cov, cov_ref); 629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // TODO: why isn't this allowed ? : 630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner<n,n>() , cov_ref); 631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct chwirut2_functor : Functor<double> 634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chwirut2_functor(void) : Functor<double>(3,54) {} 636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_x[54]; 637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_y[54]; 638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int i; 641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==54); 644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(i=0; i<54; i++) { 645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = m_x[i]; 646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; 647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==54); 654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<54; i++) { 656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = m_x[i]; 657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double factor = 1./(b[1]+b[2]*x); 658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(-b[0]*x); 659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = -x*e*factor; 660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = -e*factor*factor; 661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = -x*e*factor*factor; 662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml 670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistChwirut2(void) 671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 676c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 677c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 678c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 679c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 680c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.1, 0.01, 0.02; 681c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 682c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath chwirut2_functor functor; 683c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<chwirut2_functor> lm(functor); 684c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 10); 689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 690c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 691c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 692c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 693c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 694c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 695c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 696c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 697c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 698c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 699c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 700c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.15, 0.008, 0.010; 701c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 702c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 703c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 704c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 705c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 706c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 707c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 708c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 709c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 7); 710c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 6); 711c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 712c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.1304802941E+02); 713c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 714c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.6657666537E-01); 715c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.1653291286E-03); 716c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 1.2150007096E-02); 717c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 718c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 719c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 720c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct misra1a_functor : Functor<double> 721c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 722c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1a_functor(void) : Functor<double>(2,14) {} 723c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_x[14]; 724c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_y[14]; 725c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 726c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 727c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 728c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==14); 729c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 730c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; 731c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 732c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 733c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 734c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 735c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 736c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 737c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==14); 738c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==2); 739c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 740c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = (1.-exp(-b[1]*m_x[i])); 741c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); 742c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 743c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 744c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 745c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 746c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 747c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 748c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 749c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml 750c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMisra1a(void) 751c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 752c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=2; 753c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 754c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 755c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 756c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 757c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 758c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 759c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 760c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 500., 0.0001; 761c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 762c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1a_functor functor; 763c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<misra1a_functor> lm(functor); 764c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 765c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 766c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 767c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 768c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 19); 769c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 15); 770c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 771c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 772c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 773c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 774c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 775c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 776c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 777c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 778c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 779c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 250., 0.0005; 780c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 781c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 782c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 783c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 784c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 785c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 5); 786c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 4); 787c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 788c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.2455138894E-01); 789c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 790c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.3894212918E+02); 791c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.5015643181E-04); 792c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 793c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 794c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct hahn1_functor : Functor<double> 795c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 796c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hahn1_functor(void) : Functor<double>(7,236) {} 797c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double m_x[236]; 798c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 799c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 8007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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 , 8017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 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 , 8027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez13.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 }; 803c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 804c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 805c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 806c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 807c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==236); 808c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<236; i++) { 809c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=m_x[i], xx=x*x, xxx=xx*x; 810c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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]; 811c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 812c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 813c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 814c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 815c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 816c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 817c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 818c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==236); 819c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==7); 820c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<236; i++) { 821c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=m_x[i], xx=x*x, xxx=xx*x; 822c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 823c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.*fact; 824c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = x*fact; 825c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = xx*fact; 826c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = xxx*fact; 827c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 828c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = x*fact; 829c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,5) = xx*fact; 830c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,6) = xxx*fact; 831c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 832c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 833c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 834c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 8357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 , 8367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez282.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 , 8377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez141.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}; 838c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 839c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml 840c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistHahn1(void) 841c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 842c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=7; 843c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 844c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 845c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 846c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 847c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 848c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 849c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 850c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 10., -1., .05, -.00001, -.05, .001, -.000001; 851c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 852c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hahn1_functor functor; 853c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<hahn1_functor> lm(functor); 854c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 855c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 856c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 857c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 858c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 11); 859c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 10); 860c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 861c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 862c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 863c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.0776351733E+00); 864c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1],-1.2269296921E-01); 865c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 4.0863750610E-03); 866c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 867c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 868c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 2.4053735503E-04); 869c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6],-1.2314450199E-07); 870c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 871c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 872c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 873c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 874c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; 875c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 876c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 877c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 878c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 879c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 880c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 11); 881c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 10); 882c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 883c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.5324382854E+00); 884c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 885c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 886c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 887c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 888c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 889c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4],-5.7609940901E-03); 890c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 891c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 892c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 893c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 894c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 895c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct misra1d_functor : Functor<double> 896c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 897c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1d_functor(void) : Functor<double>(2,14) {} 898c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[14]; 899c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[14]; 900c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 901c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 902c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 903c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==14); 904c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 905c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; 906c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 907c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 908c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 909c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 910c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 911c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 912c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==14); 913c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==2); 914c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<14; i++) { 915c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double den = 1.+b[1]*x[i]; 916c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = b[1]*x[i] / den; 917c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; 918c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 919c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 920c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 921c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 922c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 923c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 924c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 925c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml 926c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMisra1d(void) 927c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 928c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=2; 929c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 930c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 931c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 932c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 933c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 934c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 935c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 936c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 500., 0.0001; 937c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 938c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath misra1d_functor functor; 939c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<misra1d_functor> lm(functor); 940c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 941c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 942c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 943c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 3); 944c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 9); 945c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 7); 946c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 947c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 948c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 949c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 950c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 951c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 952c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 953c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 954c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 955c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 450., 0.0003; 956c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 957c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 958c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 959c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 960c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 961c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 4); 962c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 3); 963c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 964c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6419295283E-02); 965c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 966c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 4.3736970754E+02); 967c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 3.0227324449E-04); 968c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 969c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 970c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 971c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct lanczos1_functor : Functor<double> 972c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 973c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lanczos1_functor(void) : Functor<double>(6,24) {} 974c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[24]; 975c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[24]; 976c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 977c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 978c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==6); 979c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==24); 980c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<24; i++) 981c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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]; 982c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 983c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 984c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 985c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 986c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==6); 987c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==24); 988c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==6); 989c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<24; i++) { 990c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = exp(-b[1]*x[i]); 991c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); 992c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = exp(-b[3]*x[i]); 993c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); 994c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = exp(-b[5]*x[i]); 995c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); 996c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 997c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 998c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 999c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1000c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1001c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1002c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1003c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml 1004c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistLanczos1(void) 1005c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1006c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=6; 1007c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1008c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1009c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1010c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1011c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1012c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1013c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1014c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; 1015c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1016c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lanczos1_functor functor; 1017c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<lanczos1_functor> lm(functor); 1018c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1019c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1020c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1021c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1022c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 79); 1023c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 72); 1024c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1025c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.430899764097e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 1026c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1027c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1028c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1029c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1030c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1031c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1032c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1033c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1034c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1035c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1036c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1037c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; 1038c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1039c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1040c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1041c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1042c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1043c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 9); 1044c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 1045c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1046c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.428595533845e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats 1047c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1048c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 9.5100000027E-02); 1049c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.0000000001E+00); 1050c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 8.6070000013E-01); 1051c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 3.0000000002E+00); 1052c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 1.5575999998E+00); 1053c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 5.0000000001E+00); 1054c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1055c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1056c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1057c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct rat42_functor : Functor<double> 1058c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1059c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat42_functor(void) : Functor<double>(3,9) {} 1060c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[9]; 1061c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[9]; 1062c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1063c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1064c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1065c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==9); 1066c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<9; i++) { 1067c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; 1068c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1069c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1070c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1071c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1072c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1073c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1074c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1075c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==9); 1076c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1077c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<9; i++) { 1078c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(b[1]-b[2]*x[i]); 1079c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1./(1.+e); 1080c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); 1081c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); 1082c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1083c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1084c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1085c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1086c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; 1087c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; 1088c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1089c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml 1090c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistRat42(void) 1091c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1092c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1093c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1094c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1095c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1096c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1097c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1098c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1099c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 100., 1., 0.1; 1101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat42_functor functor; 1103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<rat42_functor> lm(functor); 1104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 10); 1109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 1110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 75., 2.5, 0.07; 1121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 6); 1127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 5); 1128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.0565229338E+00); 1130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 7.2462237576E+01); 1132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 2.6180768402E+00); 1133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 6.7359200066E-02); 1134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct MGH10_functor : Functor<double> 1137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH10_functor(void) : Functor<double>(3,16) {} 1139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[16]; 1140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[16]; 1141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==16); 1145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<16; i++) 1146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; 1147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==16); 1153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<16; i++) { 1155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double factor = 1./(x[i]+b[2]); 1156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(b[1]*factor); 1157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = e; 1158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*factor*e; 1159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = -b[1]*b[0]*factor*factor*e; 1160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml 1168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMGH10(void) 1169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 2., 400000., 25000.; 1179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH10_functor functor; 1181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<MGH10_functor> lm(functor); 1182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 284 ); 1187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 249 ); 1188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.02, 4000., 250.; 1199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 3); 1204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 126); 1205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 116); 1206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7945855171E+01); 1208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 5.6096364710E-03); 1210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 6.1813463463E+03); 1211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 3.4522363462E+02); 1212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct BoxBOD_functor : Functor<double> 1216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BoxBOD_functor(void) : Functor<double>(2,6) {} 1218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[6]; 1219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[6] = { 109., 149., 149., 191., 213., 224. }; 1222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 1223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==6); 1224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<6; i++) 1225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; 1226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==2); 1231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==6); 1232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==2); 1233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<6; i++) { 1234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(-b[1]*x[i]); 1235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.-e; 1236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*x[i]*e; 1237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; 1242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml 1244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistBoxBOD(void) 1245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=2; 1247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1., 1.; 1255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath BoxBOD_functor functor; 1257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<BoxBOD_functor> lm(functor); 1258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.factor = 10.; 1261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 31); 1266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 25); 1267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 100., 0.75; 1277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = NumTraits<double>::epsilon(); 1280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = NumTraits<double>::epsilon(); 1281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 15 ); 1286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 14 ); 1287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.1680088766E+03); 1289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 2.1380940889E+02); 1291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.4723748542E-01); 1292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct MGH17_functor : Functor<double> 1295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH17_functor(void) : Functor<double>(5,33) {} 1297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[33]; 1298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[33]; 1299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==5); 1302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==33); 1303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<33; i++) 1304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; 1305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==5); 1310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==33); 1311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==5); 1312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<33; i++) { 1313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.; 1314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = exp(-b[3]*x[i]); 1315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = exp(-b[4]*x[i]); 1316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); 1317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); 1318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml 1326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMGH17(void) 1327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=5; 1329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 50., 150., -100., 1., 2.; 1337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH17_functor functor; 1339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<MGH17_functor> lm(functor); 1340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = NumTraits<double>::epsilon(); 1341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = NumTraits<double>::epsilon(); 1342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.maxfev = 1000; 1343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 2); 1347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 602 ); 1348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 545 ); 1349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; 1362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 18); 1369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 15); 1370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.4648946975E-05); 1372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 3.7541005211E-01); 1374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.9358469127E+00); 1375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], -1.4646871366E+00); 1376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2867534640E-02); 1377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 2.2122699662E-02); 1378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct MGH09_functor : Functor<double> 1381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH09_functor(void) : Functor<double>(4,11) {} 1383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double _x[11]; 1384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[11]; 1385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==11); 1389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<11; i++) { 1390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = _x[i], xx=x*x; 1391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; 1392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==11); 1399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==4); 1400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<11; i++) { 1401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x = _x[i], xx=x*x; 1402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double factor = 1./(xx+x*b[2]+b[3]); 1403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = (xx+x*b[1]) * factor; 1404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = b[0]*x* factor; 1405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; 1406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; 1407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml 1415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistMGH09(void) 1416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=4; 1418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 25., 39, 41.5, 39.; 1426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MGH09_functor functor; 1428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<MGH09_functor> lm(functor); 1429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.maxfev = 1000; 1430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 490 ); 1435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 376 ); 1436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 1440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 1441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 1442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 1443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 0.25, 0.39, 0.415, 0.39; 1448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 18); 1455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 16); 1456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 3.0750560385E-04); 1458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 1460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 1461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 1462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 1463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct Bennett5_functor : Functor<double> 1468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Bennett5_functor(void) : Functor<double>(3,154) {} 1470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[154]; 1471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[154]; 1472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==154); 1476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<154; i++) 1477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; 1478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==154); 1484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<154; i++) { 1486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = pow(b[1]+x[i],-1./b[2]); 1487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = e; 1488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); 1489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; 1490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 14947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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, 14957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez11.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 }; 14967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezconst 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 14977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez,-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 , 14987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; 1499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml 1501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistBennett5(void) 1502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< -2000., 50., 0.8; 1512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Bennett5_functor functor; 1514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<Bennett5_functor> lm(functor); 1515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.maxfev = 1000; 1516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 758); 1521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 744); 1522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], -2.5235058043E+03); 1526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 4.6736564644E+01); 1527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 9.3218483193E-01); 1528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< -1500., 45., 0.85; 1532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 203); 1539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 192); 1540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.2404744073E-04); 1542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 1544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); 1545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); 1546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct thurber_functor : Functor<double> 1549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath thurber_functor(void) : Functor<double>(7,37) {} 1551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double _x[37]; 1552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double _y[37]; 1553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; 1556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 1557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==37); 1558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<37; i++) { 1559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=_x[i], xx=x*x, xxx=xx*x; 1560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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]; 1561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==7); 1567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==37); 1568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==7); 1569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<37; i++) { 1570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double x=_x[i], xx=x*x, xxx=xx*x; 1571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); 1572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = 1.*fact; 1573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = x*fact; 1574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = xx*fact; 1575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = xxx*fact; 1576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; 1577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,4) = x*fact; 1578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,5) = xx*fact; 1579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,6) = xxx*fact; 1580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 1586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml 1588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistThurber(void) 1589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=7; 1591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; 1599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath thurber_functor functor; 1601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<thurber_functor> lm(functor); 1602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1608c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 39); 1609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 36); 1610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1623c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1624c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; 1625c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1626c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1627c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E4*NumTraits<double>::epsilon(); 1628c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E4*NumTraits<double>::epsilon(); 1629c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1630c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1631c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1632c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1633c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 29); 1634c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 28); 1635c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1636c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 5.6427082397E+03); 1637c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1638c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.2881396800E+03); 1639c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 1.4910792535E+03); 1640c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 5.8323836877E+02); 1641c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 7.5416644291E+01); 1642c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[4], 9.6629502864E-01); 1643c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[5], 3.9797285797E-01); 1644c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[6], 4.9727297349E-02); 1645c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1646c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1647c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct rat43_functor : Functor<double> 1648c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1649c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat43_functor(void) : Functor<double>(4,15) {} 1650c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[15]; 1651c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[15]; 1652c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1653c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1654c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1655c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==15); 1656c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<15; i++) 1657c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; 1658c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1659c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1660c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1661c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1662c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==4); 1663c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==15); 1664c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==4); 1665c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<15; i++) { 1666c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(b[1]-b[2]*x[i]); 1667c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double power = -1./b[3]; 1668c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = pow(1.+e, power); 1669c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); 1670c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); 1671c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); 1672c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1673c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1674c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1675c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1676c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathconst double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; 1677c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1678c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1679c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml 1680c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistRat43(void) 1681c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1682c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=4; 1683c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1684c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1685c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1686c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1687c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1688c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1689c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1690c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 100., 10., 1., 1.; 1691c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1692c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath rat43_functor functor; 1693c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<rat43_functor> lm(functor); 1694c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E6*NumTraits<double>::epsilon(); 1695c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E6*NumTraits<double>::epsilon(); 1696c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1697c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1698c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1699c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1700c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 27); 1701c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 20); 1702c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1703c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1704c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1705c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1706c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1707c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1708c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1709c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1710c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1711c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1712c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1713c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 700., 5., 0.75, 1.3; 1714c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1715c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.resetParameters(); 1716c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.ftol = 1.E5*NumTraits<double>::epsilon(); 1717c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lm.parameters.xtol = 1.E5*NumTraits<double>::epsilon(); 1718c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1719c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1720c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1721c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1722c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 9); 1723c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 8); 1724c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1725c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 8.7864049080E+03); 1726c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1727c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 6.9964151270E+02); 1728c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 5.2771253025E+00); 1729c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 7.5962938329E-01); 1730c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[3], 1.2792483859E+00); 1731c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1732c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1733c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1734c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1735c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct eckerle4_functor : Functor<double> 1736c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1737c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eckerle4_functor(void) : Functor<double>(3,35) {} 1738c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double x[35]; 1739c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static const double y[35]; 1740c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int operator()(const VectorXd &b, VectorXd &fvec) 1741c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1742c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1743c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fvec.size()==35); 1744c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<35; i++) 1745c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; 1746c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1747c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1748c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int df(const VectorXd &b, MatrixXd &fjac) 1749c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1750c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(b.size()==3); 1751c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.rows()==35); 1752c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(fjac.cols()==3); 1753c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<35; i++) { 1754c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double b12 = b[1]*b[1]; 1755c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); 1756c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,0) = e / b[1]; 1757c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; 1758c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; 1759c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1760c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return 0; 1761c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 1762c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 1763c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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}; 1764c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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 }; 1765c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1766c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml 1767c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid testNistEckerle4(void) 1768c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1769c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const int n=3; 1770c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int info; 1771c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1772c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorXd x(n); 1773c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1774c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1775c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * First try 1776c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1777c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1., 10., 500.; 1778c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1779c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eckerle4_functor functor; 1780c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LevenbergMarquardt<eckerle4_functor> lm(functor); 1781c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1782c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1783c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1784c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1785c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 18); 1786c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 15); 1787c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1788c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1789c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1790c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.5543827178); 1791c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 4.0888321754); 1792c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1793c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1794c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* 1795c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Second try 1796c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 1797c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath x<< 1.5, 5., 450.; 1798c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // do the computation 1799c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath info = lm.minimize(x); 1800c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1801c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check return value 1802c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(info, 1); 1803c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.nfev, 7); 1804c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_EQUAL(lm.njev, 6); 1805c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check norm^2 1806c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(lm.fvec.squaredNorm(), 1.4635887487E-03); 1807c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // check x 1808c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[0], 1.5543827178); 1809c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[1], 4.0888321754); 1810c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VERIFY_IS_APPROX(x[2], 4.5154121844E+02); 1811c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1812c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1813c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid test_NonLinearOptimization() 1814c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 1815c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Tests using the examples provided by (c)minpack 1816c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_1*/(testChkder()); 1817c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_1*/(testLmder1()); 1818c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_1*/(testLmder()); 1819c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrj1()); 1820c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrj()); 1821c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrd1()); 1822c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_2*/(testHybrd()); 1823c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmstr1()); 1824c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmstr()); 1825c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmdif1()); 1826c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_3*/(testLmdif()); 1827c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1828c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // NIST tests, level of difficulty = "Lower" 1829c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_4*/(testNistMisra1a()); 1830c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_4*/(testNistChwirut2()); 1831c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1832c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // NIST tests, level of difficulty = "Average" 1833c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_5*/(testNistHahn1()); 1834c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_6*/(testNistMisra1d()); 18357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST/*_7*/(testNistMGH17()); 18367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST/*_8*/(testNistLanczos1()); 1837c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// // NIST tests, level of difficulty = "Higher" 1839c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_9*/(testNistRat42()); 18407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST/*_10*/(testNistMGH10()); 1841c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_11*/(testNistBoxBOD()); 18427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// CALL_SUBTEST/*_12*/(testNistMGH09()); 1843c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_13*/(testNistBennett5()); 1844c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_14*/(testNistThurber()); 1845c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_15*/(testNistRat43()); 1846c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CALL_SUBTEST/*_16*/(testNistEckerle4()); 1847c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 1848c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1849c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/* 1850c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Can be useful for debugging... 1851c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev : %d, %d\n", info, lm.nfev); 1852c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev, njev : %d, %d, %d\n", info, solver.nfev, solver.njev); 1853c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev : %d, %d\n", info, solver.nfev); 1854c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[0] : %.32g\n", x[0]); 1855c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[1] : %.32g\n", x[1]); 1856c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[2] : %.32g\n", x[2]); 1857c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("x[3] : %.32g\n", x[3]); 1858c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("fvec.blueNorm() : %.32g\n", solver.fvec.blueNorm()); 1859c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("fvec.blueNorm() : %.32g\n", lm.fvec.blueNorm()); 1860c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1861c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("info, nfev, njev : %d, %d, %d\n", info, lm.nfev, lm.njev); 1862c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath printf("fvec.squaredNorm() : %.13g\n", lm.fvec.squaredNorm()); 1863c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x << std::endl; 1864c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout.precision(9); 1865c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[0] << std::endl; 1866c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[1] << std::endl; 1867c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[2] << std::endl; 1868c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::cout << x[3] << std::endl; 1869c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath*/ 1870c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 1871