From cc0fef980736302b92bd811b44afcbf523450757 Mon Sep 17 00:00:00 2001 From: Desire NUENTSA Date: Fri, 7 Dec 2012 15:48:21 +0100 Subject: [PATCH] Add tests for dense and sparse levenberg-Marquardt --- test/CMakeLists.txt | 8 + test/denseLM.cpp | 190 +++++ test/levenberg_marquardt.cpp | 1447 ++++++++++++++++++++++++++++++++++ test/sparseLM.cpp | 176 +++++ 4 files changed, 1821 insertions(+) create mode 100644 test/denseLM.cpp create mode 100644 test/levenberg_marquardt.cpp create mode 100644 test/sparseLM.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 59cd68cab..f4d9eac94 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -220,6 +220,10 @@ ei_add_test(conjugate_gradient) ei_add_test(bicgstab) ei_add_test(sparselu) +ei_add_test(levenberg_marquardt) + +# ei_add_test(denseLM) + if(UMFPACK_FOUND) ei_add_test(umfpack_support "" "${UMFPACK_ALL_LIBS}") endif() @@ -244,6 +248,10 @@ if(SPQR_FOUND AND CHOLMOD_FOUND) ei_add_test(spqr_support "" "${SPQR_ALL_LIBS}") endif() +# if(SPQR_FOUND AND CHOLMOD_FOUND) +# ei_add_test(sparseLM "" "${SPQR_ALL_LIBS}") +# endif() + string(TOLOWER "${CMAKE_CXX_COMPILER}" cmake_cxx_compiler_tolower) if(cmake_cxx_compiler_tolower MATCHES "qcc") set(CXX_IS_QCC "ON") diff --git a/test/denseLM.cpp b/test/denseLM.cpp new file mode 100644 index 000000000..0aa736ea3 --- /dev/null +++ b/test/denseLM.cpp @@ -0,0 +1,190 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa +// Copyright (C) 2012 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include +#include +#include + +#include "main.h" +#include +using namespace std; +using namespace Eigen; + +template +struct DenseLM : DenseFunctor +{ + typedef DenseFunctor Base; + typedef typename Base::JacobianType JacobianType; + typedef Matrix VectorType; + + DenseLM(int n, int m) : DenseFunctor(n,m) + { } + + VectorType model(const VectorType& uv, VectorType& x) + { + VectorType y; // Should change to use expression template + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + eigen_assert(x.size() == m); + y.setZero(m); + int half = n/2; + VectorBlock u(uv, 0, half); + VectorBlock v(uv, half, half); + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + y(j) += u(i)*std::exp(-(x(j)-i)*(x(j)-i)/(v(i)*v(i))); + } + return y; + + } + void initPoints(VectorType& uv_ref, VectorType& x) + { + m_x = x; + m_y = this->model(uv_ref, x); + } + + int operator()(const VectorType& uv, VectorType& fvec) + { + + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + eigen_assert(fvec.size() == m); + int half = n/2; + VectorBlock u(uv, 0, half); + VectorBlock v(uv, half, half); + for (int j = 0; j < m; j++) + { + fvec(j) = m_y(j); + for (int i = 0; i < half; i++) + { + fvec(j) -= u(i) *std::exp(-(m_x(j)-i)*(m_x(j)-i)/(v(i)*v(i))); + } + } + + return 0; + } + int df(const VectorType& uv, JacobianType& fjac) + { + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(n == uv.size()); + eigen_assert(fjac.rows() == m); + eigen_assert(fjac.cols() == n); + int half = n/2; + VectorBlock u(uv, 0, half); + VectorBlock v(uv, half, half); + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + { + fjac.coeffRef(j,i) = -std::exp(-(m_x(j)-i)*(m_x(j)-i)/(v(i)*v(i))); + fjac.coeffRef(j,i+half) = -2.*u(i)*(m_x(j)-i)*(m_x(j)-i)/(std::pow(v(i),3)) * std::exp(-(m_x(j)-i)*(m_x(j)-i)/(v(i)*v(i))); + } + } + return 0; + } + VectorType m_x, m_y; //Data Points +}; + +template +int test_minimizeLM(FunctorType& functor, VectorType& uv) +{ + LevenbergMarquardt lm(functor); + LevenbergMarquardtSpace::Status info; + + info = lm.minimize(uv); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template +int test_lmder(FunctorType& functor, VectorType& uv) +{ + typedef typename VectorType::Scalar Scalar; + LevenbergMarquardtSpace::Status info; + LevenbergMarquardt lm(functor); + info = lm.lmder1(uv); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template +int test_minimizeSteps(FunctorType& functor, VectorType& uv) +{ + LevenbergMarquardtSpace::Status info; + LevenbergMarquardt lm(functor); + info = lm.minimizeInit(uv); + if (info==LevenbergMarquardtSpace::ImproperInputParameters) + return info; + do + { + info = lm.minimizeOneStep(uv); + } while (info==LevenbergMarquardtSpace::Running); + + VERIFY_IS_EQUAL(info, 1); + //FIXME Check other parameters + return info; +} + +template +void test_denseLM_T() +{ + typedef Matrix VectorType; + + int inputs = 10; + int values = 1000; + DenseLM dense_gaussian(inputs, values); + VectorType uv(inputs),uv_ref(inputs); + VectorType x(values); + + // Generate the reference solution + uv_ref << -2, 1, 4 ,8, 6, 1.8, 1.2, 1.1, 1.9 , 3; + + //Generate the reference data points + x.setRandom(); + x = 10*x; + x.array() += 10; + dense_gaussian.initPoints(uv_ref, x); + + // Generate the initial parameters + VectorBlock u(uv, 0, inputs/2); + VectorBlock v(uv, inputs/2, inputs/2); + + // Solve the optimization problem + + //Solve in one go + u.setOnes(); v.setOnes(); + test_minimizeLM(dense_gaussian, uv); + + //Solve until the machine precision + u.setOnes(); v.setOnes(); + test_lmder(dense_gaussian, uv); + + // Solve step by step + v.setOnes(); u.setOnes(); + test_minimizeSteps(dense_gaussian, uv); + +} + +void test_denseLM() +{ + CALL_SUBTEST_2(test_denseLM_T()); + + // CALL_SUBTEST_2(test_sparseLM_T()); +} diff --git a/test/levenberg_marquardt.cpp b/test/levenberg_marquardt.cpp new file mode 100644 index 000000000..805b22c96 --- /dev/null +++ b/test/levenberg_marquardt.cpp @@ -0,0 +1,1447 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2009 Thomas Capricelli +// Copyright (C) 2012 desire Nuentsa + +#include "main.h" +#include + +// This disables some useless Warnings on MSVC. +// It is intended to be done for this test only. +#include + +using std::sqrt; + +struct lmder_functor : DenseFunctor +{ + lmder_functor(void): DenseFunctor(3,15) {} + int operator()(const VectorXd &x, VectorXd &fvec) const + { + double tmp1, tmp2, tmp3; + 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, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34, 2.1, 4.39}; + + for (int i = 0; i < values(); i++) + { + tmp1 = i+1; + tmp2 = 16 - i - 1; + tmp3 = (i>=8)? tmp2 : tmp1; + fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); + } + return 0; + } + + int df(const VectorXd &x, MatrixXd &fjac) const + { + double tmp1, tmp2, tmp3, tmp4; + for (int i = 0; i < values(); i++) + { + tmp1 = i+1; + tmp2 = 16 - i - 1; + tmp3 = (i>=8)? tmp2 : tmp1; + tmp4 = (x[1]*tmp2 + x[2]*tmp3); tmp4 = tmp4*tmp4; + fjac(i,0) = -1; + fjac(i,1) = tmp1*tmp2/tmp4; + fjac(i,2) = tmp1*tmp3/tmp4; + } + return 0; + } +}; + +void testLmder1() +{ + int n=3, info; + + VectorXd x; + + /* the following starting values provide a rough fit. */ + x.setConstant(n, 1.); + + // do the computation + lmder_functor functor; + LevenbergMarquardt lm(functor); + info = lm.lmder1(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 6); + VERIFY_IS_EQUAL(lm.njev(), 5); + + // check norm + VERIFY_IS_APPROX(lm.fvec().blueNorm(), 0.09063596); + + // check x + VectorXd x_ref(n); + x_ref << 0.08241058, 1.133037, 2.343695; + VERIFY_IS_APPROX(x, x_ref); +} + +void testLmder() +{ + const int m=15, n=3; + int info; + double fnorm, covfac; + VectorXd x; + + /* the following starting values provide a rough fit. */ + x.setConstant(n, 1.); + + // do the computation + lmder_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return values + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 6); + VERIFY_IS_EQUAL(lm.njev(), 5); + + // check norm + fnorm = lm.fvec().blueNorm(); + VERIFY_IS_APPROX(fnorm, 0.09063596); + + // check x + VectorXd x_ref(n); + x_ref << 0.08241058, 1.133037, 2.343695; + VERIFY_IS_APPROX(x, x_ref); + + // check covariance + covfac = fnorm*fnorm/(m-n); + internal::covar(lm.fjac(), lm.permutation().indices()); // TODO : move this as a function of lm + + MatrixXd cov_ref(n,n); + cov_ref << + 0.0001531202, 0.002869941, -0.002656662, + 0.002869941, 0.09480935, -0.09098995, + -0.002656662, -0.09098995, 0.08778727; + +// std::cout << fjac*covfac << std::endl; + + MatrixXd cov; + cov = covfac*lm.fjac().topLeftCorner(); + VERIFY_IS_APPROX( cov, cov_ref); + // TODO: why isn't this allowed ? : + // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner() , cov_ref); +} + +struct lmdif_functor : DenseFunctor +{ + lmdif_functor(void) : DenseFunctor(3,15) {} + int operator()(const VectorXd &x, VectorXd &fvec) const + { + int i; + double tmp1,tmp2,tmp3; + 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, + 3.7e-1,5.8e-1,7.3e-1,9.6e-1,1.34e0,2.1e0,4.39e0}; + + assert(x.size()==3); + assert(fvec.size()==15); + for (i=0; i<15; i++) + { + tmp1 = i+1; + tmp2 = 15 - i; + tmp3 = tmp1; + + if (i >= 8) tmp3 = tmp2; + fvec[i] = y[i] - (x[0] + tmp1/(x[1]*tmp2 + x[2]*tmp3)); + } + return 0; + } +}; + +void testLmdif1() +{ + const int n=3; + int info; + + VectorXd x(n), fvec(15); + + /* the following starting values provide a rough fit. */ + x.setConstant(n, 1.); + + // do the computation + lmdif_functor functor; + DenseIndex nfev; + info = LevenbergMarquardt::lmdif1(functor, x, &nfev); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(nfev, 26); + + // check norm + functor(x, fvec); + VERIFY_IS_APPROX(fvec.blueNorm(), 0.09063596); + + // check x + VectorXd x_ref(n); + x_ref << 0.0824106, 1.1330366, 2.3436947; + VERIFY_IS_APPROX(x, x_ref); + +} + +void testLmdif() +{ + const int m=15, n=3; + int info; + double fnorm, covfac; + VectorXd x(n); + + /* the following starting values provide a rough fit. */ + x.setConstant(n, 1.); + + // do the computation + lmdif_functor functor; + NumericalDiff numDiff(functor); + LevenbergMarquardt > lm(numDiff); + info = lm.minimize(x); + + // check return values + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 26); + + // check norm + fnorm = lm.fvec().blueNorm(); + VERIFY_IS_APPROX(fnorm, 0.09063596); + + // check x + VectorXd x_ref(n); + x_ref << 0.08241058, 1.133037, 2.343695; + VERIFY_IS_APPROX(x, x_ref); + + // check covariance + covfac = fnorm*fnorm/(m-n); + internal::covar(lm.fjac(), lm.permutation().indices()); // TODO : move this as a function of lm + + MatrixXd cov_ref(n,n); + cov_ref << + 0.0001531202, 0.002869942, -0.002656662, + 0.002869942, 0.09480937, -0.09098997, + -0.002656662, -0.09098997, 0.08778729; + +// std::cout << fjac*covfac << std::endl; + + MatrixXd cov; + cov = covfac*lm.fjac().topLeftCorner(); + VERIFY_IS_APPROX( cov, cov_ref); + // TODO: why isn't this allowed ? : + // VERIFY_IS_APPROX( covfac*fjac.topLeftCorner() , cov_ref); +} + +struct chwirut2_functor : DenseFunctor +{ + chwirut2_functor(void) : DenseFunctor(3,54) {} + static const double m_x[54]; + static const double m_y[54]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + int i; + + assert(b.size()==3); + assert(fvec.size()==54); + for(i=0; i<54; i++) { + double x = m_x[i]; + fvec[i] = exp(-b[0]*x)/(b[1]+b[2]*x) - m_y[i]; + } + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); + assert(fjac.rows()==54); + assert(fjac.cols()==3); + for(int i=0; i<54; i++) { + double x = m_x[i]; + double factor = 1./(b[1]+b[2]*x); + double e = exp(-b[0]*x); + fjac(i,0) = -x*e*factor; + fjac(i,1) = -e*factor*factor; + fjac(i,2) = -x*e*factor*factor; + } + return 0; + } +}; +const 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}; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/chwirut2.shtml +void testNistChwirut2(void) +{ + const int n=3; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 0.1, 0.01, 0.02; + // do the computation + chwirut2_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 10); + VERIFY_IS_EQUAL(lm.njev(), 8); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); + // check x + VERIFY_IS_APPROX(x[0], 1.6657666537E-01); + VERIFY_IS_APPROX(x[1], 5.1653291286E-03); + VERIFY_IS_APPROX(x[2], 1.2150007096E-02); + + /* + * Second try + */ + x<< 0.15, 0.008, 0.010; + // do the computation + lm.resetParameters(); + lm.setFtol(1.E6*NumTraits::epsilon()); + lm.setXtol(1.E6*NumTraits::epsilon()); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 7); + VERIFY_IS_EQUAL(lm.njev(), 6); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.1304802941E+02); + // check x + VERIFY_IS_APPROX(x[0], 1.6657666537E-01); + VERIFY_IS_APPROX(x[1], 5.1653291286E-03); + VERIFY_IS_APPROX(x[2], 1.2150007096E-02); +} + + +struct misra1a_functor : DenseFunctor +{ + misra1a_functor(void) : DenseFunctor(2,14) {} + static const double m_x[14]; + static const double m_y[14]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==2); + assert(fvec.size()==14); + for(int i=0; i<14; i++) { + fvec[i] = b[0]*(1.-exp(-b[1]*m_x[i])) - m_y[i] ; + } + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==2); + assert(fjac.rows()==14); + assert(fjac.cols()==2); + for(int i=0; i<14; i++) { + fjac(i,0) = (1.-exp(-b[1]*m_x[i])); + fjac(i,1) = (b[0]*m_x[i]*exp(-b[1]*m_x[i])); + } + return 0; + } +}; +const 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}; +const 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}; + +// http://www.itl.nist.gov/div898/strd/nls/data/misra1a.shtml +void testNistMisra1a(void) +{ + const int n=2; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 500., 0.0001; + // do the computation + misra1a_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 19); + VERIFY_IS_EQUAL(lm.njev(), 15); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); + // check x + VERIFY_IS_APPROX(x[0], 2.3894212918E+02); + VERIFY_IS_APPROX(x[1], 5.5015643181E-04); + + /* + * Second try + */ + x<< 250., 0.0005; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 5); + VERIFY_IS_EQUAL(lm.njev(), 4); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.2455138894E-01); + // check x + VERIFY_IS_APPROX(x[0], 2.3894212918E+02); + VERIFY_IS_APPROX(x[1], 5.5015643181E-04); +} + +struct hahn1_functor : DenseFunctor +{ + hahn1_functor(void) : DenseFunctor(7,236) {} + static const double m_x[236]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + static const double m_y[236] = { .591E0 , 1.547E0 , 2.902E0 , 2.894E0 , 4.703E0 , 6.307E0 , 7.03E0 , 7.898E0 , 9.470E0 , 9.484E0 , 10.072E0 , 10.163E0 , 11.615E0 , 12.005E0 , 12.478E0 , 12.982E0 , 12.970E0 , 13.926E0 , 14.452E0 , 14.404E0 , 15.190E0 , 15.550E0 , 15.528E0 , 15.499E0 , 16.131E0 , 16.438E0 , 16.387E0 , 16.549E0 , 16.872E0 , 16.830E0 , 16.926E0 , 16.907E0 , 16.966E0 , 17.060E0 , 17.122E0 , 17.311E0 , 17.355E0 , 17.668E0 , 17.767E0 , 17.803E0 , 17.765E0 , 17.768E0 , 17.736E0 , 17.858E0 , 17.877E0 , 17.912E0 , 18.046E0 , 18.085E0 , 18.291E0 , 18.357E0 , 18.426E0 , 18.584E0 , 18.610E0 , 18.870E0 , 18.795E0 , 19.111E0 , .367E0 , .796E0 , 0.892E0 , 1.903E0 , 2.150E0 , 3.697E0 , 5.870E0 , 6.421E0 , 7.422E0 , 9.944E0 , 11.023E0 , 11.87E0 , 12.786E0 , 14.067E0 , 13.974E0 , 14.462E0 , 14.464E0 , 15.381E0 , 15.483E0 , 15.59E0 , 16.075E0 , 16.347E0 , 16.181E0 , 16.915E0 , 17.003E0 , 16.978E0 , 17.756E0 , 17.808E0 , 17.868E0 , 18.481E0 , 18.486E0 , 19.090E0 , 16.062E0 , 16.337E0 , 16.345E0 , + 16.388E0 , 17.159E0 , 17.116E0 , 17.164E0 , 17.123E0 , 17.979E0 , 17.974E0 , 18.007E0 , 17.993E0 , 18.523E0 , 18.669E0 , 18.617E0 , 19.371E0 , 19.330E0 , 0.080E0 , 0.248E0 , 1.089E0 , 1.418E0 , 2.278E0 , 3.624E0 , 4.574E0 , 5.556E0 , 7.267E0 , 7.695E0 , 9.136E0 , 9.959E0 , 9.957E0 , 11.600E0 , 13.138E0 , 13.564E0 , 13.871E0 , 13.994E0 , 14.947E0 , 15.473E0 , 15.379E0 , 15.455E0 , 15.908E0 , 16.114E0 , 17.071E0 , 17.135E0 , 17.282E0 , 17.368E0 , 17.483E0 , 17.764E0 , 18.185E0 , 18.271E0 , 18.236E0 , 18.237E0 , 18.523E0 , 18.627E0 , 18.665E0 , 19.086E0 , 0.214E0 , 0.943E0 , 1.429E0 , 2.241E0 , 2.951E0 , 3.782E0 , 4.757E0 , 5.602E0 , 7.169E0 , 8.920E0 , 10.055E0 , 12.035E0 , 12.861E0 , 13.436E0 , 14.167E0 , 14.755E0 , 15.168E0 , 15.651E0 , 15.746E0 , 16.216E0 , 16.445E0 , 16.965E0 , 17.121E0 , 17.206E0 , 17.250E0 , 17.339E0 , 17.793E0 , 18.123E0 , 18.49E0 , 18.566E0 , 18.645E0 , 18.706E0 , 18.924E0 , 19.1E0 , 0.375E0 , 0.471E0 , 1.504E0 , 2.204E0 , 2.813E0 , 4.765E0 , 9.835E0 , 10.040E0 , 11.946E0 , 12.596E0 , +13.303E0 , 13.922E0 , 14.440E0 , 14.951E0 , 15.627E0 , 15.639E0 , 15.814E0 , 16.315E0 , 16.334E0 , 16.430E0 , 16.423E0 , 17.024E0 , 17.009E0 , 17.165E0 , 17.134E0 , 17.349E0 , 17.576E0 , 17.848E0 , 18.090E0 , 18.276E0 , 18.404E0 , 18.519E0 , 19.133E0 , 19.074E0 , 19.239E0 , 19.280E0 , 19.101E0 , 19.398E0 , 19.252E0 , 19.89E0 , 20.007E0 , 19.929E0 , 19.268E0 , 19.324E0 , 20.049E0 , 20.107E0 , 20.062E0 , 20.065E0 , 19.286E0 , 19.972E0 , 20.088E0 , 20.743E0 , 20.83E0 , 20.935E0 , 21.035E0 , 20.93E0 , 21.074E0 , 21.085E0 , 20.935E0 }; + + // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; + + assert(b.size()==7); + assert(fvec.size()==236); + for(int i=0; i<236; i++) { + double x=m_x[i], xx=x*x, xxx=xx*x; + 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]; + } + return 0; + } + + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==7); + assert(fjac.rows()==236); + assert(fjac.cols()==7); + for(int i=0; i<236; i++) { + double x=m_x[i], xx=x*x, xxx=xx*x; + double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); + fjac(i,0) = 1.*fact; + fjac(i,1) = x*fact; + fjac(i,2) = xx*fact; + fjac(i,3) = xxx*fact; + fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; + fjac(i,4) = x*fact; + fjac(i,5) = xx*fact; + fjac(i,6) = xxx*fact; + } + return 0; + } +}; +const double hahn1_functor::m_x[236] = { 24.41E0 , 34.82E0 , 44.09E0 , 45.07E0 , 54.98E0 , 65.51E0 , 70.53E0 , 75.70E0 , 89.57E0 , 91.14E0 , 96.40E0 , 97.19E0 , 114.26E0 , 120.25E0 , 127.08E0 , 133.55E0 , 133.61E0 , 158.67E0 , 172.74E0 , 171.31E0 , 202.14E0 , 220.55E0 , 221.05E0 , 221.39E0 , 250.99E0 , 268.99E0 , 271.80E0 , 271.97E0 , 321.31E0 , 321.69E0 , 330.14E0 , 333.03E0 , 333.47E0 , 340.77E0 , 345.65E0 , 373.11E0 , 373.79E0 , 411.82E0 , 419.51E0 , 421.59E0 , 422.02E0 , 422.47E0 , 422.61E0 , 441.75E0 , 447.41E0 , 448.7E0 , 472.89E0 , 476.69E0 , 522.47E0 , 522.62E0 , 524.43E0 , 546.75E0 , 549.53E0 , 575.29E0 , 576.00E0 , 625.55E0 , 20.15E0 , 28.78E0 , 29.57E0 , 37.41E0 , 39.12E0 , 50.24E0 , 61.38E0 , 66.25E0 , 73.42E0 , 95.52E0 , 107.32E0 , 122.04E0 , 134.03E0 , 163.19E0 , 163.48E0 , 175.70E0 , 179.86E0 , 211.27E0 , 217.78E0 , 219.14E0 , 262.52E0 , 268.01E0 , 268.62E0 , 336.25E0 , 337.23E0 , 339.33E0 , 427.38E0 , 428.58E0 , 432.68E0 , 528.99E0 , 531.08E0 , 628.34E0 , 253.24E0 , 273.13E0 , 273.66E0 , +282.10E0 , 346.62E0 , 347.19E0 , 348.78E0 , 351.18E0 , 450.10E0 , 450.35E0 , 451.92E0 , 455.56E0 , 552.22E0 , 553.56E0 , 555.74E0 , 652.59E0 , 656.20E0 , 14.13E0 , 20.41E0 , 31.30E0 , 33.84E0 , 39.70E0 , 48.83E0 , 54.50E0 , 60.41E0 , 72.77E0 , 75.25E0 , 86.84E0 , 94.88E0 , 96.40E0 , 117.37E0 , 139.08E0 , 147.73E0 , 158.63E0 , 161.84E0 , 192.11E0 , 206.76E0 , 209.07E0 , 213.32E0 , 226.44E0 , 237.12E0 , 330.90E0 , 358.72E0 , 370.77E0 , 372.72E0 , 396.24E0 , 416.59E0 , 484.02E0 , 495.47E0 , 514.78E0 , 515.65E0 , 519.47E0 , 544.47E0 , 560.11E0 , 620.77E0 , 18.97E0 , 28.93E0 , 33.91E0 , 40.03E0 , 44.66E0 , 49.87E0 , 55.16E0 , 60.90E0 , 72.08E0 , 85.15E0 , 97.06E0 , 119.63E0 , 133.27E0 , 143.84E0 , 161.91E0 , 180.67E0 , 198.44E0 , 226.86E0 , 229.65E0 , 258.27E0 , 273.77E0 , 339.15E0 , 350.13E0 , 362.75E0 , 371.03E0 , 393.32E0 , 448.53E0 , 473.78E0 , 511.12E0 , 524.70E0 , 548.75E0 , 551.64E0 , 574.02E0 , 623.86E0 , 21.46E0 , 24.33E0 , 33.43E0 , 39.22E0 , 44.18E0 , 55.02E0 , 94.33E0 , 96.44E0 , 118.82E0 , 128.48E0 , +141.94E0 , 156.92E0 , 171.65E0 , 190.00E0 , 223.26E0 , 223.88E0 , 231.50E0 , 265.05E0 , 269.44E0 , 271.78E0 , 273.46E0 , 334.61E0 , 339.79E0 , 349.52E0 , 358.18E0 , 377.98E0 , 394.77E0 , 429.66E0 , 468.22E0 , 487.27E0 , 519.54E0 , 523.03E0 , 612.99E0 , 638.59E0 , 641.36E0 , 622.05E0 , 631.50E0 , 663.97E0 , 646.9E0 , 748.29E0 , 749.21E0 , 750.14E0 , 647.04E0 , 646.89E0 , 746.9E0 , 748.43E0 , 747.35E0 , 749.27E0 , 647.61E0 , 747.78E0 , 750.51E0 , 851.37E0 , 845.97E0 , 847.54E0 , 849.93E0 , 851.61E0 , 849.75E0 , 850.98E0 , 848.23E0}; + +// http://www.itl.nist.gov/div898/strd/nls/data/hahn1.shtml +void testNistHahn1(void) +{ + const int n=7; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 10., -1., .05, -.00001, -.05, .001, -.000001; + // do the computation + hahn1_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 11); + VERIFY_IS_EQUAL(lm.njev(), 10); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); + // check x + VERIFY_IS_APPROX(x[0], 1.0776351733E+00); + VERIFY_IS_APPROX(x[1],-1.2269296921E-01); + VERIFY_IS_APPROX(x[2], 4.0863750610E-03); + VERIFY_IS_APPROX(x[3],-1.426264e-06); // shoulde be : -1.4262662514E-06 + VERIFY_IS_APPROX(x[4],-5.7609940901E-03); + VERIFY_IS_APPROX(x[5], 2.4053735503E-04); + VERIFY_IS_APPROX(x[6],-1.2314450199E-07); + + /* + * Second try + */ + x<< .1, -.1, .005, -.000001, -.005, .0001, -.0000001; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 11); + VERIFY_IS_EQUAL(lm.njev(), 10); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.5324382854E+00); + // check x + VERIFY_IS_APPROX(x[0], 1.077640); // should be : 1.0776351733E+00 + VERIFY_IS_APPROX(x[1], -0.1226933); // should be : -1.2269296921E-01 + VERIFY_IS_APPROX(x[2], 0.004086383); // should be : 4.0863750610E-03 + VERIFY_IS_APPROX(x[3], -1.426277e-06); // shoulde be : -1.4262662514E-06 + VERIFY_IS_APPROX(x[4],-5.7609940901E-03); + VERIFY_IS_APPROX(x[5], 0.00024053772); // should be : 2.4053735503E-04 + VERIFY_IS_APPROX(x[6], -1.231450e-07); // should be : -1.2314450199E-07 + +} + +struct misra1d_functor : DenseFunctor +{ + misra1d_functor(void) : DenseFunctor(2,14) {} + static const double x[14]; + static const double y[14]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==2); + assert(fvec.size()==14); + for(int i=0; i<14; i++) { + fvec[i] = b[0]*b[1]*x[i]/(1.+b[1]*x[i]) - y[i]; + } + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==2); + assert(fjac.rows()==14); + assert(fjac.cols()==2); + for(int i=0; i<14; i++) { + double den = 1.+b[1]*x[i]; + fjac(i,0) = b[1]*x[i] / den; + fjac(i,1) = b[0]*x[i]*(den-b[1]*x[i])/den/den; + } + return 0; + } +}; +const 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}; +const 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}; + +// http://www.itl.nist.gov/div898/strd/nls/data/misra1d.shtml +void testNistMisra1d(void) +{ + const int n=2; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 500., 0.0001; + // do the computation + misra1d_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 9); + VERIFY_IS_EQUAL(lm.njev(), 7); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); + // check x + VERIFY_IS_APPROX(x[0], 4.3736970754E+02); + VERIFY_IS_APPROX(x[1], 3.0227324449E-04); + + /* + * Second try + */ + x<< 450., 0.0003; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 4); + VERIFY_IS_EQUAL(lm.njev(), 3); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6419295283E-02); + // check x + VERIFY_IS_APPROX(x[0], 4.3736970754E+02); + VERIFY_IS_APPROX(x[1], 3.0227324449E-04); +} + + +struct lanczos1_functor : DenseFunctor +{ + lanczos1_functor(void) : DenseFunctor(6,24) {} + static const double x[24]; + static const double y[24]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==6); + assert(fvec.size()==24); + for(int i=0; i<24; i++) + 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]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==6); + assert(fjac.rows()==24); + assert(fjac.cols()==6); + for(int i=0; i<24; i++) { + fjac(i,0) = exp(-b[1]*x[i]); + fjac(i,1) = -b[0]*x[i]*exp(-b[1]*x[i]); + fjac(i,2) = exp(-b[3]*x[i]); + fjac(i,3) = -b[2]*x[i]*exp(-b[3]*x[i]); + fjac(i,4) = exp(-b[5]*x[i]); + fjac(i,5) = -b[4]*x[i]*exp(-b[5]*x[i]); + } + return 0; + } +}; +const 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 }; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/lanczos1.shtml +void testNistLanczos1(void) +{ + const int n=6; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 1.2, 0.3, 5.6, 5.5, 6.5, 7.6; + // do the computation + lanczos1_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 2); + VERIFY_IS_EQUAL(lm.nfev(), 79); + VERIFY_IS_EQUAL(lm.njev(), 72); + // check norm^2 +// VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.430899764097e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats + // check x + VERIFY_IS_APPROX(x[0], 9.5100000027E-02); + VERIFY_IS_APPROX(x[1], 1.0000000001E+00); + VERIFY_IS_APPROX(x[2], 8.6070000013E-01); + VERIFY_IS_APPROX(x[3], 3.0000000002E+00); + VERIFY_IS_APPROX(x[4], 1.5575999998E+00); + VERIFY_IS_APPROX(x[5], 5.0000000001E+00); + + /* + * Second try + */ + x<< 0.5, 0.7, 3.6, 4.2, 4., 6.3; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 2); + VERIFY_IS_EQUAL(lm.nfev(), 9); + VERIFY_IS_EQUAL(lm.njev(), 8); + // check norm^2 +// VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.428595533845e-25); // should be 1.4307867721E-25, but nist results are on 128-bit floats + // check x + VERIFY_IS_APPROX(x[0], 9.5100000027E-02); + VERIFY_IS_APPROX(x[1], 1.0000000001E+00); + VERIFY_IS_APPROX(x[2], 8.6070000013E-01); + VERIFY_IS_APPROX(x[3], 3.0000000002E+00); + VERIFY_IS_APPROX(x[4], 1.5575999998E+00); + VERIFY_IS_APPROX(x[5], 5.0000000001E+00); + +} + +struct rat42_functor : DenseFunctor +{ + rat42_functor(void) : DenseFunctor(3,9) {} + static const double x[9]; + static const double y[9]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==3); + assert(fvec.size()==9); + for(int i=0; i<9; i++) { + fvec[i] = b[0] / (1.+exp(b[1]-b[2]*x[i])) - y[i]; + } + return 0; + } + + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); + assert(fjac.rows()==9); + assert(fjac.cols()==3); + for(int i=0; i<9; i++) { + double e = exp(b[1]-b[2]*x[i]); + fjac(i,0) = 1./(1.+e); + fjac(i,1) = -b[0]*e/(1.+e)/(1.+e); + fjac(i,2) = +b[0]*e*x[i]/(1.+e)/(1.+e); + } + return 0; + } +}; +const double rat42_functor::x[9] = { 9.000E0, 14.000E0, 21.000E0, 28.000E0, 42.000E0, 57.000E0, 63.000E0, 70.000E0, 79.000E0 }; +const double rat42_functor::y[9] = { 8.930E0 ,10.800E0 ,18.590E0 ,22.330E0 ,39.350E0 ,56.110E0 ,61.730E0 ,64.620E0 ,67.080E0 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky2.shtml +void testNistRat42(void) +{ + const int n=3; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 100., 1., 0.1; + // do the computation + rat42_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 10); + VERIFY_IS_EQUAL(lm.njev(), 8); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); + // check x + VERIFY_IS_APPROX(x[0], 7.2462237576E+01); + VERIFY_IS_APPROX(x[1], 2.6180768402E+00); + VERIFY_IS_APPROX(x[2], 6.7359200066E-02); + + /* + * Second try + */ + x<< 75., 2.5, 0.07; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 6); + VERIFY_IS_EQUAL(lm.njev(), 5); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.0565229338E+00); + // check x + VERIFY_IS_APPROX(x[0], 7.2462237576E+01); + VERIFY_IS_APPROX(x[1], 2.6180768402E+00); + VERIFY_IS_APPROX(x[2], 6.7359200066E-02); +} + +struct MGH10_functor : DenseFunctor +{ + MGH10_functor(void) : DenseFunctor(3,16) {} + static const double x[16]; + static const double y[16]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==3); + assert(fvec.size()==16); + for(int i=0; i<16; i++) + fvec[i] = b[0] * exp(b[1]/(x[i]+b[2])) - y[i]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); + assert(fjac.rows()==16); + assert(fjac.cols()==3); + for(int i=0; i<16; i++) { + double factor = 1./(x[i]+b[2]); + double e = exp(b[1]*factor); + fjac(i,0) = e; + fjac(i,1) = b[0]*factor*e; + fjac(i,2) = -b[1]*b[0]*factor*factor*e; + } + return 0; + } +}; +const 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 }; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/mgh10.shtml +void testNistMGH10(void) +{ + const int n=3; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 2., 400000., 25000.; + // do the computation + MGH10_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 284 ); + VERIFY_IS_EQUAL(lm.njev(), 249 ); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); + // check x + VERIFY_IS_APPROX(x[0], 5.6096364710E-03); + VERIFY_IS_APPROX(x[1], 6.1813463463E+03); + VERIFY_IS_APPROX(x[2], 3.4522363462E+02); + + /* + * Second try + */ + x<< 0.02, 4000., 250.; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 126); + VERIFY_IS_EQUAL(lm.njev(), 116); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7945855171E+01); + // check x + VERIFY_IS_APPROX(x[0], 5.6096364710E-03); + VERIFY_IS_APPROX(x[1], 6.1813463463E+03); + VERIFY_IS_APPROX(x[2], 3.4522363462E+02); +} + + +struct BoxBOD_functor : DenseFunctor +{ + BoxBOD_functor(void) : DenseFunctor(2,6) {} + static const double x[6]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + static const double y[6] = { 109., 149., 149., 191., 213., 224. }; + assert(b.size()==2); + assert(fvec.size()==6); + for(int i=0; i<6; i++) + fvec[i] = b[0]*(1.-exp(-b[1]*x[i])) - y[i]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==2); + assert(fjac.rows()==6); + assert(fjac.cols()==2); + for(int i=0; i<6; i++) { + double e = exp(-b[1]*x[i]); + fjac(i,0) = 1.-e; + fjac(i,1) = b[0]*x[i]*e; + } + return 0; + } +}; +const double BoxBOD_functor::x[6] = { 1., 2., 3., 5., 7., 10. }; + +// http://www.itl.nist.gov/div898/strd/nls/data/boxbod.shtml +void testNistBoxBOD(void) +{ + const int n=2; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 1., 1.; + // do the computation + BoxBOD_functor functor; + LevenbergMarquardt lm(functor); + lm.setFtol(1.E6*NumTraits::epsilon()); + lm.setXtol(1.E6*NumTraits::epsilon()); + lm.setFactor(10); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 31); + VERIFY_IS_EQUAL(lm.njev(), 25); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); + // check x + VERIFY_IS_APPROX(x[0], 2.1380940889E+02); + VERIFY_IS_APPROX(x[1], 5.4723748542E-01); + + /* + * Second try + */ + x<< 100., 0.75; + // do the computation + lm.resetParameters(); + lm.setFtol(NumTraits::epsilon()); + lm.setXtol( NumTraits::epsilon()); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 15 ); + VERIFY_IS_EQUAL(lm.njev(), 14 ); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.1680088766E+03); + // check x + VERIFY_IS_APPROX(x[0], 2.1380940889E+02); + VERIFY_IS_APPROX(x[1], 5.4723748542E-01); +} + +struct MGH17_functor : DenseFunctor +{ + MGH17_functor(void) : DenseFunctor(5,33) {} + static const double x[33]; + static const double y[33]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==5); + assert(fvec.size()==33); + for(int i=0; i<33; i++) + fvec[i] = b[0] + b[1]*exp(-b[3]*x[i]) + b[2]*exp(-b[4]*x[i]) - y[i]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==5); + assert(fjac.rows()==33); + assert(fjac.cols()==5); + for(int i=0; i<33; i++) { + fjac(i,0) = 1.; + fjac(i,1) = exp(-b[3]*x[i]); + fjac(i,2) = exp(-b[4]*x[i]); + fjac(i,3) = -x[i]*b[1]*exp(-b[3]*x[i]); + fjac(i,4) = -x[i]*b[2]*exp(-b[4]*x[i]); + } + return 0; + } +}; +const 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 }; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/mgh17.shtml +void testNistMGH17(void) +{ + const int n=5; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 50., 150., -100., 1., 2.; + // do the computation + MGH17_functor functor; + LevenbergMarquardt lm(functor); + lm.setFtol(NumTraits::epsilon()); + lm.setXtol(NumTraits::epsilon()); + lm.setMaxfev(1000); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 2); +// VERIFY_IS_EQUAL(lm.nfev(), 602 ); + VERIFY_IS_EQUAL(lm.njev(), 545 ); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); + // check x + VERIFY_IS_APPROX(x[0], 3.7541005211E-01); + VERIFY_IS_APPROX(x[1], 1.9358469127E+00); + VERIFY_IS_APPROX(x[2], -1.4646871366E+00); + VERIFY_IS_APPROX(x[3], 1.2867534640E-02); + VERIFY_IS_APPROX(x[4], 2.2122699662E-02); + + /* + * Second try + */ + x<< 0.5 ,1.5 ,-1 ,0.01 ,0.02; + // do the computation + lm.resetParameters(); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 18); + VERIFY_IS_EQUAL(lm.njev(), 15); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.4648946975E-05); + // check x + VERIFY_IS_APPROX(x[0], 3.7541005211E-01); + VERIFY_IS_APPROX(x[1], 1.9358469127E+00); + VERIFY_IS_APPROX(x[2], -1.4646871366E+00); + VERIFY_IS_APPROX(x[3], 1.2867534640E-02); + VERIFY_IS_APPROX(x[4], 2.2122699662E-02); +} + +struct MGH09_functor : DenseFunctor +{ + MGH09_functor(void) : DenseFunctor(4,11) {} + static const double _x[11]; + static const double y[11]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==4); + assert(fvec.size()==11); + for(int i=0; i<11; i++) { + double x = _x[i], xx=x*x; + fvec[i] = b[0]*(xx+x*b[1])/(xx+x*b[2]+b[3]) - y[i]; + } + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==4); + assert(fjac.rows()==11); + assert(fjac.cols()==4); + for(int i=0; i<11; i++) { + double x = _x[i], xx=x*x; + double factor = 1./(xx+x*b[2]+b[3]); + fjac(i,0) = (xx+x*b[1]) * factor; + fjac(i,1) = b[0]*x* factor; + fjac(i,2) = - b[0]*(xx+x*b[1]) * x * factor * factor; + fjac(i,3) = - b[0]*(xx+x*b[1]) * factor * factor; + } + return 0; + } +}; +const 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 }; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/mgh09.shtml +void testNistMGH09(void) +{ + const int n=4; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 25., 39, 41.5, 39.; + // do the computation + MGH09_functor functor; + LevenbergMarquardt lm(functor); + lm.setMaxfev(1000); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 490 ); + VERIFY_IS_EQUAL(lm.njev(), 376 ); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); + // check x + VERIFY_IS_APPROX(x[0], 0.1928077089); // should be 1.9280693458E-01 + VERIFY_IS_APPROX(x[1], 0.19126423573); // should be 1.9128232873E-01 + VERIFY_IS_APPROX(x[2], 0.12305309914); // should be 1.2305650693E-01 + VERIFY_IS_APPROX(x[3], 0.13605395375); // should be 1.3606233068E-01 + + /* + * Second try + */ + x<< 0.25, 0.39, 0.415, 0.39; + // do the computation + lm.resetParameters(); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 18); + VERIFY_IS_EQUAL(lm.njev(), 16); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 3.0750560385E-04); + // check x + VERIFY_IS_APPROX(x[0], 0.19280781); // should be 1.9280693458E-01 + VERIFY_IS_APPROX(x[1], 0.19126265); // should be 1.9128232873E-01 + VERIFY_IS_APPROX(x[2], 0.12305280); // should be 1.2305650693E-01 + VERIFY_IS_APPROX(x[3], 0.13605322); // should be 1.3606233068E-01 +} + + + +struct Bennett5_functor : DenseFunctor +{ + Bennett5_functor(void) : DenseFunctor(3,154) {} + static const double x[154]; + static const double y[154]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==3); + assert(fvec.size()==154); + for(int i=0; i<154; i++) + fvec[i] = b[0]* pow(b[1]+x[i],-1./b[2]) - y[i]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); + assert(fjac.rows()==154); + assert(fjac.cols()==3); + for(int i=0; i<154; i++) { + double e = pow(b[1]+x[i],-1./b[2]); + fjac(i,0) = e; + fjac(i,1) = - b[0]*e/b[2]/(b[1]+x[i]); + fjac(i,2) = b[0]*e*log(b[1]+x[i])/b[2]/b[2]; + } + return 0; + } +}; +const double Bennett5_functor::x[154] = { 7.447168E0, 8.102586E0, 8.452547E0, 8.711278E0, 8.916774E0, 9.087155E0, 9.232590E0, 9.359535E0, 9.472166E0, 9.573384E0, 9.665293E0, 9.749461E0, 9.827092E0, 9.899128E0, 9.966321E0, 10.029280E0, 10.088510E0, 10.144430E0, 10.197380E0, 10.247670E0, 10.295560E0, 10.341250E0, 10.384950E0, 10.426820E0, 10.467000E0, 10.505640E0, 10.542830E0, 10.578690E0, 10.613310E0, 10.646780E0, 10.679150E0, 10.710520E0, 10.740920E0, 10.770440E0, 10.799100E0, 10.826970E0, 10.854080E0, 10.880470E0, 10.906190E0, 10.931260E0, 10.955720E0, 10.979590E0, 11.002910E0, 11.025700E0, 11.047980E0, 11.069770E0, 11.091100E0, 11.111980E0, 11.132440E0, 11.152480E0, 11.172130E0, 11.191410E0, 11.210310E0, 11.228870E0, 11.247090E0, 11.264980E0, 11.282560E0, 11.299840E0, 11.316820E0, 11.333520E0, 11.349940E0, 11.366100E0, 11.382000E0, 11.397660E0, 11.413070E0, 11.428240E0, 11.443200E0, 11.457930E0, 11.472440E0, 11.486750E0, 11.500860E0, 11.514770E0, 11.528490E0, 11.542020E0, 11.555380E0, 11.568550E0, +11.581560E0, 11.594420E0, 11.607121E0, 11.619640E0, 11.632000E0, 11.644210E0, 11.656280E0, 11.668200E0, 11.679980E0, 11.691620E0, 11.703130E0, 11.714510E0, 11.725760E0, 11.736880E0, 11.747890E0, 11.758780E0, 11.769550E0, 11.780200E0, 11.790730E0, 11.801160E0, 11.811480E0, 11.821700E0, 11.831810E0, 11.841820E0, 11.851730E0, 11.861550E0, 11.871270E0, 11.880890E0, 11.890420E0, 11.899870E0, 11.909220E0, 11.918490E0, 11.927680E0, 11.936780E0, 11.945790E0, 11.954730E0, 11.963590E0, 11.972370E0, 11.981070E0, 11.989700E0, 11.998260E0, 12.006740E0, 12.015150E0, 12.023490E0, 12.031760E0, 12.039970E0, 12.048100E0, 12.056170E0, 12.064180E0, 12.072120E0, 12.080010E0, 12.087820E0, 12.095580E0, 12.103280E0, 12.110920E0, 12.118500E0, 12.126030E0, 12.133500E0, 12.140910E0, 12.148270E0, 12.155570E0, 12.162830E0, 12.170030E0, 12.177170E0, 12.184270E0, 12.191320E0, 12.198320E0, 12.205270E0, 12.212170E0, 12.219030E0, 12.225840E0, 12.232600E0, 12.239320E0, 12.245990E0, 12.252620E0, 12.259200E0, 12.265750E0, 12.272240E0 }; +const double Bennett5_functor::y[154] = { -34.834702E0 ,-34.393200E0 ,-34.152901E0 ,-33.979099E0 ,-33.845901E0 ,-33.732899E0 ,-33.640301E0 ,-33.559200E0 ,-33.486801E0 ,-33.423100E0 ,-33.365101E0 ,-33.313000E0 ,-33.260899E0 ,-33.217400E0 ,-33.176899E0 ,-33.139198E0 ,-33.101601E0 ,-33.066799E0 ,-33.035000E0 ,-33.003101E0 ,-32.971298E0 ,-32.942299E0 ,-32.916302E0 ,-32.890202E0 ,-32.864101E0 ,-32.841000E0 ,-32.817799E0 ,-32.797501E0 ,-32.774300E0 ,-32.757000E0 ,-32.733799E0 ,-32.716400E0 ,-32.699100E0 ,-32.678799E0 ,-32.661400E0 ,-32.644001E0 ,-32.626701E0 ,-32.612202E0 ,-32.597698E0 ,-32.583199E0 ,-32.568699E0 ,-32.554298E0 ,-32.539799E0 ,-32.525299E0 ,-32.510799E0 ,-32.499199E0 ,-32.487598E0 ,-32.473202E0 ,-32.461601E0 ,-32.435501E0 ,-32.435501E0 ,-32.426800E0 ,-32.412300E0 ,-32.400799E0 ,-32.392101E0 ,-32.380501E0 ,-32.366001E0 ,-32.357300E0 ,-32.348598E0 ,-32.339901E0 ,-32.328400E0 ,-32.319698E0 ,-32.311001E0 ,-32.299400E0 ,-32.290699E0 ,-32.282001E0 ,-32.273300E0 ,-32.264599E0 ,-32.256001E0 ,-32.247299E0 +,-32.238602E0 ,-32.229900E0 ,-32.224098E0 ,-32.215401E0 ,-32.203800E0 ,-32.198002E0 ,-32.189400E0 ,-32.183601E0 ,-32.174900E0 ,-32.169102E0 ,-32.163300E0 ,-32.154598E0 ,-32.145901E0 ,-32.140099E0 ,-32.131401E0 ,-32.125599E0 ,-32.119801E0 ,-32.111198E0 ,-32.105400E0 ,-32.096699E0 ,-32.090900E0 ,-32.088001E0 ,-32.079300E0 ,-32.073502E0 ,-32.067699E0 ,-32.061901E0 ,-32.056099E0 ,-32.050301E0 ,-32.044498E0 ,-32.038799E0 ,-32.033001E0 ,-32.027199E0 ,-32.024300E0 ,-32.018501E0 ,-32.012699E0 ,-32.004002E0 ,-32.001099E0 ,-31.995300E0 ,-31.989500E0 ,-31.983700E0 ,-31.977900E0 ,-31.972099E0 ,-31.969299E0 ,-31.963501E0 ,-31.957701E0 ,-31.951900E0 ,-31.946100E0 ,-31.940300E0 ,-31.937401E0 ,-31.931601E0 ,-31.925800E0 ,-31.922899E0 ,-31.917101E0 ,-31.911301E0 ,-31.908400E0 ,-31.902599E0 ,-31.896900E0 ,-31.893999E0 ,-31.888201E0 ,-31.885300E0 ,-31.882401E0 ,-31.876600E0 ,-31.873699E0 ,-31.867901E0 ,-31.862101E0 ,-31.859200E0 ,-31.856300E0 ,-31.850500E0 ,-31.844700E0 ,-31.841801E0 ,-31.838900E0 ,-31.833099E0 ,-31.830200E0 , +-31.827299E0 ,-31.821600E0 ,-31.818701E0 ,-31.812901E0 ,-31.809999E0 ,-31.807100E0 ,-31.801300E0 ,-31.798401E0 ,-31.795500E0 ,-31.789700E0 ,-31.786800E0 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/bennett5.shtml +void testNistBennett5(void) +{ + const int n=3; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< -2000., 50., 0.8; + // do the computation + Bennett5_functor functor; + LevenbergMarquardt lm(functor); + lm.setMaxfev(1000); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 758); + VERIFY_IS_EQUAL(lm.njev(), 744); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); + // check x + VERIFY_IS_APPROX(x[0], -2.5235058043E+03); + VERIFY_IS_APPROX(x[1], 4.6736564644E+01); + VERIFY_IS_APPROX(x[2], 9.3218483193E-01); + /* + * Second try + */ + x<< -1500., 45., 0.85; + // do the computation + lm.resetParameters(); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 203); + VERIFY_IS_EQUAL(lm.njev(), 192); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.2404744073E-04); + // check x + VERIFY_IS_APPROX(x[0], -2523.3007865); // should be -2.5235058043E+03 + VERIFY_IS_APPROX(x[1], 46.735705771); // should be 4.6736564644E+01); + VERIFY_IS_APPROX(x[2], 0.93219881891); // should be 9.3218483193E-01); +} + +struct thurber_functor : DenseFunctor +{ + thurber_functor(void) : DenseFunctor(7,37) {} + static const double _x[37]; + static const double _y[37]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + // int called=0; printf("call hahn1_functor with iflag=%d, called=%d\n", iflag, called); if (iflag==1) called++; + assert(b.size()==7); + assert(fvec.size()==37); + for(int i=0; i<37; i++) { + double x=_x[i], xx=x*x, xxx=xx*x; + 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]; + } + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==7); + assert(fjac.rows()==37); + assert(fjac.cols()==7); + for(int i=0; i<37; i++) { + double x=_x[i], xx=x*x, xxx=xx*x; + double fact = 1./(1.+b[4]*x+b[5]*xx+b[6]*xxx); + fjac(i,0) = 1.*fact; + fjac(i,1) = x*fact; + fjac(i,2) = xx*fact; + fjac(i,3) = xxx*fact; + fact = - (b[0]+b[1]*x+b[2]*xx+b[3]*xxx) * fact * fact; + fjac(i,4) = x*fact; + fjac(i,5) = xx*fact; + fjac(i,6) = xxx*fact; + } + return 0; + } +}; +const 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 }; +const 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}; + +// http://www.itl.nist.gov/div898/strd/nls/data/thurber.shtml +void testNistThurber(void) +{ + const int n=7; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 1000 ,1000 ,400 ,40 ,0.7,0.3,0.0 ; + // do the computation + thurber_functor functor; + LevenbergMarquardt lm(functor); + lm.setFtol(1.E4*NumTraits::epsilon()); + lm.setXtol(1.E4*NumTraits::epsilon()); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 39); + VERIFY_IS_EQUAL(lm.njev(), 36); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); + // check x + VERIFY_IS_APPROX(x[0], 1.2881396800E+03); + VERIFY_IS_APPROX(x[1], 1.4910792535E+03); + VERIFY_IS_APPROX(x[2], 5.8323836877E+02); + VERIFY_IS_APPROX(x[3], 7.5416644291E+01); + VERIFY_IS_APPROX(x[4], 9.6629502864E-01); + VERIFY_IS_APPROX(x[5], 3.9797285797E-01); + VERIFY_IS_APPROX(x[6], 4.9727297349E-02); + + /* + * Second try + */ + x<< 1300 ,1500 ,500 ,75 ,1 ,0.4 ,0.05 ; + // do the computation + lm.resetParameters(); + lm.setFtol(1.E4*NumTraits::epsilon()); + lm.setXtol(1.E4*NumTraits::epsilon()); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 29); + VERIFY_IS_EQUAL(lm.njev(), 28); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 5.6427082397E+03); + // check x + VERIFY_IS_APPROX(x[0], 1.2881396800E+03); + VERIFY_IS_APPROX(x[1], 1.4910792535E+03); + VERIFY_IS_APPROX(x[2], 5.8323836877E+02); + VERIFY_IS_APPROX(x[3], 7.5416644291E+01); + VERIFY_IS_APPROX(x[4], 9.6629502864E-01); + VERIFY_IS_APPROX(x[5], 3.9797285797E-01); + VERIFY_IS_APPROX(x[6], 4.9727297349E-02); +} + +struct rat43_functor : DenseFunctor +{ + rat43_functor(void) : DenseFunctor(4,15) {} + static const double x[15]; + static const double y[15]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==4); + assert(fvec.size()==15); + for(int i=0; i<15; i++) + fvec[i] = b[0] * pow(1.+exp(b[1]-b[2]*x[i]),-1./b[3]) - y[i]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==4); + assert(fjac.rows()==15); + assert(fjac.cols()==4); + for(int i=0; i<15; i++) { + double e = exp(b[1]-b[2]*x[i]); + double power = -1./b[3]; + fjac(i,0) = pow(1.+e, power); + fjac(i,1) = power*b[0]*e*pow(1.+e, power-1.); + fjac(i,2) = -power*b[0]*e*x[i]*pow(1.+e, power-1.); + fjac(i,3) = b[0]*power*power*log(1.+e)*pow(1.+e, power); + } + return 0; + } +}; +const double rat43_functor::x[15] = { 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15. }; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/ratkowsky3.shtml +void testNistRat43(void) +{ + const int n=4; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 100., 10., 1., 1.; + // do the computation + rat43_functor functor; + LevenbergMarquardt lm(functor); + lm.setFtol(1.E6*NumTraits::epsilon()); + lm.setXtol(1.E6*NumTraits::epsilon()); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 27); + VERIFY_IS_EQUAL(lm.njev(), 20); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); + // check x + VERIFY_IS_APPROX(x[0], 6.9964151270E+02); + VERIFY_IS_APPROX(x[1], 5.2771253025E+00); + VERIFY_IS_APPROX(x[2], 7.5962938329E-01); + VERIFY_IS_APPROX(x[3], 1.2792483859E+00); + + /* + * Second try + */ + x<< 700., 5., 0.75, 1.3; + // do the computation + lm.resetParameters(); + lm.setFtol(1.E5*NumTraits::epsilon()); + lm.setXtol(1.E5*NumTraits::epsilon()); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 9); + VERIFY_IS_EQUAL(lm.njev(), 8); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 8.7864049080E+03); + // check x + VERIFY_IS_APPROX(x[0], 6.9964151270E+02); + VERIFY_IS_APPROX(x[1], 5.2771253025E+00); + VERIFY_IS_APPROX(x[2], 7.5962938329E-01); + VERIFY_IS_APPROX(x[3], 1.2792483859E+00); +} + + + +struct eckerle4_functor : DenseFunctor +{ + eckerle4_functor(void) : DenseFunctor(3,35) {} + static const double x[35]; + static const double y[35]; + int operator()(const VectorXd &b, VectorXd &fvec) + { + assert(b.size()==3); + assert(fvec.size()==35); + for(int i=0; i<35; i++) + fvec[i] = b[0]/b[1] * exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/(b[1]*b[1])) - y[i]; + return 0; + } + int df(const VectorXd &b, MatrixXd &fjac) + { + assert(b.size()==3); + assert(fjac.rows()==35); + assert(fjac.cols()==3); + for(int i=0; i<35; i++) { + double b12 = b[1]*b[1]; + double e = exp(-0.5*(x[i]-b[2])*(x[i]-b[2])/b12); + fjac(i,0) = e / b[1]; + fjac(i,1) = ((x[i]-b[2])*(x[i]-b[2])/b12-1.) * b[0]*e/b12; + fjac(i,2) = (x[i]-b[2])*e*b[0]/b[1]/b12; + } + return 0; + } +}; +const 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}; +const 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 }; + +// http://www.itl.nist.gov/div898/strd/nls/data/eckerle4.shtml +void testNistEckerle4(void) +{ + const int n=3; + int info; + + VectorXd x(n); + + /* + * First try + */ + x<< 1., 10., 500.; + // do the computation + eckerle4_functor functor; + LevenbergMarquardt lm(functor); + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 18); + VERIFY_IS_EQUAL(lm.njev(), 15); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); + // check x + VERIFY_IS_APPROX(x[0], 1.5543827178); + VERIFY_IS_APPROX(x[1], 4.0888321754); + VERIFY_IS_APPROX(x[2], 4.5154121844E+02); + + /* + * Second try + */ + x<< 1.5, 5., 450.; + // do the computation + info = lm.minimize(x); + + // check return value + VERIFY_IS_EQUAL(info, 1); + VERIFY_IS_EQUAL(lm.nfev(), 7); + VERIFY_IS_EQUAL(lm.njev(), 6); + // check norm^2 + VERIFY_IS_APPROX(lm.fvec().squaredNorm(), 1.4635887487E-03); + // check x + VERIFY_IS_APPROX(x[0], 1.5543827178); + VERIFY_IS_APPROX(x[1], 4.0888321754); + VERIFY_IS_APPROX(x[2], 4.5154121844E+02); +} + +void test_levenberg_marquardt() +{ + // Tests using the examples provided by (c)minpack + CALL_SUBTEST(testLmder1()); + CALL_SUBTEST(testLmder()); + CALL_SUBTEST(testLmdif1()); +// CALL_SUBTEST(testLmstr1()); +// CALL_SUBTEST(testLmstr()); + CALL_SUBTEST(testLmdif()); + + // NIST tests, level of difficulty = "Lower" + CALL_SUBTEST(testNistMisra1a()); + CALL_SUBTEST(testNistChwirut2()); + + // NIST tests, level of difficulty = "Average" + CALL_SUBTEST(testNistHahn1()); + CALL_SUBTEST(testNistMisra1d()); + CALL_SUBTEST(testNistMGH17()); + CALL_SUBTEST(testNistLanczos1()); + +// // NIST tests, level of difficulty = "Higher" + CALL_SUBTEST(testNistRat42()); + CALL_SUBTEST(testNistMGH10()); + CALL_SUBTEST(testNistBoxBOD()); +// CALL_SUBTEST(testNistMGH09()); + CALL_SUBTEST(testNistBennett5()); + CALL_SUBTEST(testNistThurber()); + CALL_SUBTEST(testNistRat43()); + CALL_SUBTEST(testNistEckerle4()); +} diff --git a/test/sparseLM.cpp b/test/sparseLM.cpp new file mode 100644 index 000000000..8e148f9bc --- /dev/null +++ b/test/sparseLM.cpp @@ -0,0 +1,176 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// Copyright (C) 2012 Desire Nuentsa +// Copyright (C) 2012 Gael Guennebaud +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include +#include +#include + +#include "main.h" +#include + +using namespace std; +using namespace Eigen; + +template +struct sparseGaussianTest : SparseFunctor +{ + typedef Matrix VectorType; + typedef SparseFunctor Base; + typedef typename Base::JacobianType JacobianType; + sparseGaussianTest(int inputs, int values) : SparseFunctor(inputs,values) + { } + + VectorType model(const VectorType& uv, VectorType& x) + { + VectorType y; //Change this to use expression template + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + eigen_assert(x.size() == m); + y.setZero(m); + int half = n/2; + VectorBlock u(uv, 0, half); + VectorBlock v(uv, half, half); + Scalar coeff; + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + { + coeff = (x(j)-i)/v(i); + coeff *= coeff; + if (coeff < 1. && coeff > 0.) + y(j) += u(i)*std::pow((1-coeff), 2); + } + } + return y; + } + void initPoints(VectorType& uv_ref, VectorType& x) + { + m_x = x; + m_y = this->model(uv_ref,x); + } + int operator()(const VectorType& uv, VectorType& fvec) + { + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(uv.size()%2 == 0); + eigen_assert(uv.size() == n); + int half = n/2; + VectorBlock u(uv, 0, half); + VectorBlock v(uv, half, half); + fvec = m_y; + Scalar coeff; + for (int j = 0; j < m; j++) + { + for (int i = 0; i < half; i++) + { + coeff = (m_x(j)-i)/v(i); + coeff *= coeff; + if (coeff < 1. && coeff > 0.) + fvec(j) -= u(i)*std::pow((1-coeff), 2); + } + } + return 0; + } + + int df(const VectorType& uv, JacobianType& fjac) + { + int m = Base::values(); + int n = Base::inputs(); + eigen_assert(n == uv.size()); + eigen_assert(fjac.rows() == m); + eigen_assert(fjac.cols() == n); + int half = n/2; + VectorBlock u(uv, 0, half); + VectorBlock v(uv, half, half); + Scalar coeff; + + //Derivatives with respect to u + for (int col = 0; col < half; col++) + { + for (int row = 0; row < m; row++) + { + coeff = (m_x(row)-col)/v(col); + coeff = coeff*coeff; + if(coeff < 1. && coeff > 0.) + { + fjac.coeffRef(row,col) = -(1-coeff)*(1-coeff); + } + } + } + //Derivatives with respect to v + for (int col = 0; col < half; col++) + { + for (int row = 0; row < m; row++) + { + coeff = (m_x(row)-col)/v(col); + coeff = coeff*coeff; + if(coeff < 1. && coeff > 0.) + { + fjac.coeffRef(row,col+half) = -4 * (u(col)/v(col))*coeff*(1-coeff); + } + } + } + return 0; + } + + VectorType m_x, m_y; //Data points +}; + + +template +void test_sparseLM_T() +{ + typedef Matrix VectorType; + + int inputs = 10; + int values = 2000; + sparseGaussianTest sparse_gaussian(inputs, values); + VectorType uv(inputs),uv_ref(inputs); + VectorType x(values); + // Generate the reference solution + uv_ref << -2, 1, 4 ,8, 6, 1.8, 1.2, 1.1, 1.9 , 3; + //Generate the reference data points + x.setRandom(); + x = 10*x; + x.array() += 10; + sparse_gaussian.initPoints(uv_ref, x); + + + // Generate the initial parameters + VectorBlock u(uv, 0, inputs/2); + VectorBlock v(uv, inputs/2, inputs/2); + v.setOnes(); + //Generate u or Solve for u from v + u.setOnes(); + + // Solve the optimization problem + LevenbergMarquardt > lm(sparse_gaussian); + int info; +// info = lm.minimize(uv); + + VERIFY_IS_EQUAL(info,1); + // Do a step by step solution and save the residual + int maxiter = 200; + int iter = 0; + MatrixXd Err(values, maxiter); + MatrixXd Mod(values, maxiter); + LevenbergMarquardtSpace::Status status; + status = lm.minimizeInit(uv); + if (status==LevenbergMarquardtSpace::ImproperInputParameters) + return ; + +} +void test_sparseLM() +{ + CALL_SUBTEST_1(test_sparseLM_T()); + + // CALL_SUBTEST_2(test_sparseLM_T()); +}