From a7d085eb4ecedd45f091eeadb93277c7f3878b27 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Sun, 21 Feb 2010 12:41:37 +0100 Subject: NonLinearOptimization : clean 'mode' handling from the old minpack code : * this is actually a boolean, not an int * use a better name * can be set at initialization time instead of bloating all methods signatures --- .../NonLinearOptimization/HybridNonLinearSolver.h | 108 +++++++-------------- .../src/NonLinearOptimization/LevenbergMarquardt.h | 95 ++++++------------ unsupported/test/NonLinearOptimization.cpp | 6 +- 3 files changed, 68 insertions(+), 141 deletions(-) diff --git a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h index 0754a426b..35dc332e0 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h +++ b/unsupported/Eigen/src/NonLinearOptimization/HybridNonLinearSolver.h @@ -57,7 +57,7 @@ class HybridNonLinearSolver { public: HybridNonLinearSolver(FunctorType &_functor) - : functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; } + : functor(_functor) { nfev=njev=iter = 0; fnorm= 0.; useExternalScaling=false;} struct Parameters { Parameters() @@ -84,36 +84,18 @@ public: const Scalar tol = ei_sqrt(NumTraits::epsilon()) ); - HybridNonLinearSolverSpace::Status solveInit( - FVectorType &x, - const int mode=1 - ); - HybridNonLinearSolverSpace::Status solveOneStep( - FVectorType &x, - const int mode=1 - ); - HybridNonLinearSolverSpace::Status solve( - FVectorType &x, - const int mode=1 - ); + HybridNonLinearSolverSpace::Status solveInit(FVectorType &x); + HybridNonLinearSolverSpace::Status solveOneStep(FVectorType &x); + HybridNonLinearSolverSpace::Status solve(FVectorType &x); HybridNonLinearSolverSpace::Status hybrd1( FVectorType &x, const Scalar tol = ei_sqrt(NumTraits::epsilon()) ); - HybridNonLinearSolverSpace::Status solveNumericalDiffInit( - FVectorType &x, - const int mode=1 - ); - HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep( - FVectorType &x, - const int mode=1 - ); - HybridNonLinearSolverSpace::Status solveNumericalDiff( - FVectorType &x, - const int mode=1 - ); + HybridNonLinearSolverSpace::Status solveNumericalDiffInit(FVectorType &x); + HybridNonLinearSolverSpace::Status solveNumericalDiffOneStep(FVectorType &x); + HybridNonLinearSolverSpace::Status solveNumericalDiff(FVectorType &x); void resetParameters(void) { parameters = Parameters(); } Parameters parameters; @@ -124,6 +106,7 @@ public: int njev; int iter; Scalar fnorm; + bool useExternalScaling; private: FunctorType &functor; int n; @@ -160,18 +143,13 @@ HybridNonLinearSolver::hybrj1( parameters.maxfev = 100*(n+1); parameters.xtol = tol; diag.setConstant(n, 1.); - return solve( - x, - 2 - ); + useExternalScaling = true; + return solve(x); } template HybridNonLinearSolverSpace::Status -HybridNonLinearSolver::solveInit( - FVectorType &x, - const int mode - ) +HybridNonLinearSolver::solveInit(FVectorType &x) { n = x.size(); @@ -179,9 +157,9 @@ HybridNonLinearSolver::solveInit( fvec.resize(n); qtf.resize(n); fjac.resize(n, n); - if (mode != 2) + if (!useExternalScaling) diag.resize(n); - assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'"); + assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'"); /* Function Body */ nfev = 0; @@ -190,7 +168,7 @@ HybridNonLinearSolver::solveInit( /* check the input parameters for errors. */ if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0. ) return HybridNonLinearSolverSpace::ImproperInputParameters; - if (mode == 2) + if (useExternalScaling) for (int j = 0; j < n; ++j) if (diag[j] <= 0.) return HybridNonLinearSolverSpace::ImproperInputParameters; @@ -214,10 +192,7 @@ HybridNonLinearSolver::solveInit( template HybridNonLinearSolverSpace::Status -HybridNonLinearSolver::solveOneStep( - FVectorType &x, - const int mode - ) +HybridNonLinearSolver::solveOneStep(FVectorType &x) { int j; std::vector > v_givens(n), w_givens(n); @@ -231,10 +206,10 @@ HybridNonLinearSolver::solveOneStep( wa2 = fjac.colwise().blueNorm(); - /* on the first iteration and if mode is 1, scale according */ + /* on the first iteration and if external scaling is not used, scale according */ /* to the norms of the columns of the initial jacobian. */ if (iter == 1) { - if (mode != 2) + if (!useExternalScaling) for (j = 0; j < n; ++j) diag[j] = (wa2[j]==0.) ? 1. : wa2[j]; @@ -260,7 +235,7 @@ HybridNonLinearSolver::solveOneStep( qtf = fjac.transpose() * fvec; /* rescale if necessary. */ - if (mode != 2) + if (!useExternalScaling) diag = diag.cwiseMax(wa2); while (true) { @@ -372,14 +347,11 @@ HybridNonLinearSolver::solveOneStep( template HybridNonLinearSolverSpace::Status -HybridNonLinearSolver::solve( - FVectorType &x, - const int mode - ) +HybridNonLinearSolver::solve(FVectorType &x) { - HybridNonLinearSolverSpace::Status status = solveInit(x, mode); + HybridNonLinearSolverSpace::Status status = solveInit(x); while (status==HybridNonLinearSolverSpace::Running) - status = solveOneStep(x, mode); + status = solveOneStep(x); return status; } @@ -403,18 +375,13 @@ HybridNonLinearSolver::hybrd1( parameters.xtol = tol; diag.setConstant(n, 1.); - return solveNumericalDiff( - x, - 2 - ); + useExternalScaling = true; + return solveNumericalDiff(x); } template HybridNonLinearSolverSpace::Status -HybridNonLinearSolver::solveNumericalDiffInit( - FVectorType &x, - const int mode - ) +HybridNonLinearSolver::solveNumericalDiffInit(FVectorType &x) { n = x.size(); @@ -425,10 +392,9 @@ HybridNonLinearSolver::solveNumericalDiffInit( qtf.resize(n); fjac.resize(n, n); fvec.resize(n); - if (mode != 2) + if (!useExternalScaling) diag.resize(n); - assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'"); - + assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'"); /* Function Body */ nfev = 0; @@ -437,7 +403,7 @@ HybridNonLinearSolver::solveNumericalDiffInit( /* check the input parameters for errors. */ if (n <= 0 || parameters.xtol < 0. || parameters.maxfev <= 0 || parameters.nb_of_subdiagonals< 0 || parameters.nb_of_superdiagonals< 0 || parameters.factor <= 0. ) return HybridNonLinearSolverSpace::ImproperInputParameters; - if (mode == 2) + if (useExternalScaling) for (int j = 0; j < n; ++j) if (diag[j] <= 0.) return HybridNonLinearSolverSpace::ImproperInputParameters; @@ -461,10 +427,7 @@ HybridNonLinearSolver::solveNumericalDiffInit( template HybridNonLinearSolverSpace::Status -HybridNonLinearSolver::solveNumericalDiffOneStep( - FVectorType &x, - const int mode - ) +HybridNonLinearSolver::solveNumericalDiffOneStep(FVectorType &x) { int j; std::vector > v_givens(n), w_givens(n); @@ -480,10 +443,10 @@ HybridNonLinearSolver::solveNumericalDiffOneStep( wa2 = fjac.colwise().blueNorm(); - /* on the first iteration and if mode is 1, scale according */ + /* on the first iteration and if external scaling is not used, scale according */ /* to the norms of the columns of the initial jacobian. */ if (iter == 1) { - if (mode != 2) + if (!useExternalScaling) for (j = 0; j < n; ++j) diag[j] = (wa2[j]==0.) ? 1. : wa2[j]; @@ -509,7 +472,7 @@ HybridNonLinearSolver::solveNumericalDiffOneStep( qtf = fjac.transpose() * fvec; /* rescale if necessary. */ - if (mode != 2) + if (!useExternalScaling) diag = diag.cwiseMax(wa2); while (true) { @@ -621,14 +584,11 @@ HybridNonLinearSolver::solveNumericalDiffOneStep( template HybridNonLinearSolverSpace::Status -HybridNonLinearSolver::solveNumericalDiff( - FVectorType &x, - const int mode - ) +HybridNonLinearSolver::solveNumericalDiff(FVectorType &x) { - HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x, mode); + HybridNonLinearSolverSpace::Status status = solveNumericalDiffInit(x); while (status==HybridNonLinearSolverSpace::Running) - status = solveNumericalDiffOneStep(x, mode); + status = solveNumericalDiffOneStep(x); return status; } diff --git a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h index f21331e3e..8bae1e131 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h +++ b/unsupported/Eigen/src/NonLinearOptimization/LevenbergMarquardt.h @@ -61,7 +61,7 @@ class LevenbergMarquardt { public: LevenbergMarquardt(FunctorType &_functor) - : functor(_functor) { nfev = njev = iter = 0; fnorm=gnorm = 0.; } + : functor(_functor) { nfev = njev = iter = 0; fnorm = gnorm = 0.; useExternalScaling=false; } struct Parameters { Parameters() @@ -87,18 +87,9 @@ public: const Scalar tol = ei_sqrt(NumTraits::epsilon()) ); - LevenbergMarquardtSpace::Status minimize( - FVectorType &x, - const int mode=1 - ); - LevenbergMarquardtSpace::Status minimizeInit( - FVectorType &x, - const int mode=1 - ); - LevenbergMarquardtSpace::Status minimizeOneStep( - FVectorType &x, - const int mode=1 - ); + LevenbergMarquardtSpace::Status minimize(FVectorType &x); + LevenbergMarquardtSpace::Status minimizeInit(FVectorType &x); + LevenbergMarquardtSpace::Status minimizeOneStep(FVectorType &x); static LevenbergMarquardtSpace::Status lmdif1( FunctorType &functor, @@ -112,18 +103,9 @@ public: const Scalar tol = ei_sqrt(NumTraits::epsilon()) ); - LevenbergMarquardtSpace::Status minimizeOptimumStorage( - FVectorType &x, - const int mode=1 - ); - LevenbergMarquardtSpace::Status minimizeOptimumStorageInit( - FVectorType &x, - const int mode=1 - ); - LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep( - FVectorType &x, - const int mode=1 - ); + LevenbergMarquardtSpace::Status minimizeOptimumStorage(FVectorType &x); + LevenbergMarquardtSpace::Status minimizeOptimumStorageInit(FVectorType &x); + LevenbergMarquardtSpace::Status minimizeOptimumStorageOneStep(FVectorType &x); void resetParameters(void) { parameters = Parameters(); } @@ -135,6 +117,7 @@ public: int njev; int iter; Scalar fnorm, gnorm; + bool useExternalScaling; Scalar lm_param(void) { return par; } private: @@ -175,24 +158,18 @@ LevenbergMarquardt::lmder1( template LevenbergMarquardtSpace::Status -LevenbergMarquardt::minimize( - FVectorType &x, - const int mode - ) +LevenbergMarquardt::minimize(FVectorType &x) { - LevenbergMarquardtSpace::Status status = minimizeInit(x, mode); + LevenbergMarquardtSpace::Status status = minimizeInit(x); do { - status = minimizeOneStep(x, mode); + status = minimizeOneStep(x); } while (status==LevenbergMarquardtSpace::Running); return status; } template LevenbergMarquardtSpace::Status -LevenbergMarquardt::minimizeInit( - FVectorType &x, - const int mode - ) +LevenbergMarquardt::minimizeInit(FVectorType &x) { n = x.size(); m = functor.values(); @@ -201,9 +178,9 @@ LevenbergMarquardt::minimizeInit( wa4.resize(m); fvec.resize(m); fjac.resize(m, n); - if (mode != 2) + if (!useExternalScaling) diag.resize(n); - assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'"); + assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'"); qtf.resize(n); /* Function Body */ @@ -214,7 +191,7 @@ LevenbergMarquardt::minimizeInit( if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.) return LevenbergMarquardtSpace::ImproperInputParameters; - if (mode == 2) + if (useExternalScaling) for (int j = 0; j < n; ++j) if (diag[j] <= 0.) return LevenbergMarquardtSpace::ImproperInputParameters; @@ -235,10 +212,7 @@ LevenbergMarquardt::minimizeInit( template LevenbergMarquardtSpace::Status -LevenbergMarquardt::minimizeOneStep( - FVectorType &x, - const int mode - ) +LevenbergMarquardt::minimizeOneStep(FVectorType &x) { int j; @@ -257,10 +231,10 @@ LevenbergMarquardt::minimizeOneStep( fjac = qrfac.matrixQR(); permutation = qrfac.colsPermutation(); - /* on the first iteration and if mode is 1, scale according */ + /* on the first iteration and if external scaling is not used, scale according */ /* to the norms of the columns of the initial jacobian. */ if (iter == 1) { - if (mode != 2) + if (!useExternalScaling) for (j = 0; j < n; ++j) diag[j] = (wa2[j]==0.)? 1. : wa2[j]; @@ -290,7 +264,7 @@ LevenbergMarquardt::minimizeOneStep( return LevenbergMarquardtSpace::CosinusTooSmall; /* rescale if necessary. */ - if (mode != 2) + if (!useExternalScaling) diag = diag.cwiseMax(wa2); do { @@ -406,10 +380,7 @@ LevenbergMarquardt::lmstr1( template LevenbergMarquardtSpace::Status -LevenbergMarquardt::minimizeOptimumStorageInit( - FVectorType &x, - const int mode - ) +LevenbergMarquardt::minimizeOptimumStorageInit(FVectorType &x) { n = x.size(); m = functor.values(); @@ -423,9 +394,9 @@ LevenbergMarquardt::minimizeOptimumStorageInit( // The purpose it to only use a nxn matrix, instead of mxn here, so // that we can handle cases where m>>n : fjac.resize(n, n); - if (mode != 2) + if (!useExternalScaling) diag.resize(n); - assert( (mode!=2 || diag.size()==n) || "When using mode==2, the caller must provide a valid 'diag'"); + assert( (!useExternalScaling || diag.size()==n) || "When useExternalScaling is set, the caller must provide a valid 'diag'"); qtf.resize(n); /* Function Body */ @@ -436,7 +407,7 @@ LevenbergMarquardt::minimizeOptimumStorageInit( if (n <= 0 || m < n || parameters.ftol < 0. || parameters.xtol < 0. || parameters.gtol < 0. || parameters.maxfev <= 0 || parameters.factor <= 0.) return LevenbergMarquardtSpace::ImproperInputParameters; - if (mode == 2) + if (useExternalScaling) for (int j = 0; j < n; ++j) if (diag[j] <= 0.) return LevenbergMarquardtSpace::ImproperInputParameters; @@ -458,10 +429,7 @@ LevenbergMarquardt::minimizeOptimumStorageInit( template LevenbergMarquardtSpace::Status -LevenbergMarquardt::minimizeOptimumStorageOneStep( - FVectorType &x, - const int mode - ) +LevenbergMarquardt::minimizeOptimumStorageOneStep(FVectorType &x) { int i, j; bool sing; @@ -514,10 +482,10 @@ LevenbergMarquardt::minimizeOptimumStorageOneStep( } } - /* on the first iteration and if mode is 1, scale according */ + /* on the first iteration and if external scaling is not used, scale according */ /* to the norms of the columns of the initial jacobian. */ if (iter == 1) { - if (mode != 2) + if (!useExternalScaling) for (j = 0; j < n; ++j) diag[j] = (wa2[j]==0.)? 1. : wa2[j]; @@ -541,7 +509,7 @@ LevenbergMarquardt::minimizeOptimumStorageOneStep( return LevenbergMarquardtSpace::CosinusTooSmall; /* rescale if necessary. */ - if (mode != 2) + if (!useExternalScaling) diag = diag.cwiseMax(wa2); do { @@ -635,14 +603,11 @@ LevenbergMarquardt::minimizeOptimumStorageOneStep( template LevenbergMarquardtSpace::Status -LevenbergMarquardt::minimizeOptimumStorage( - FVectorType &x, - const int mode - ) +LevenbergMarquardt::minimizeOptimumStorage(FVectorType &x) { - LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x, mode); + LevenbergMarquardtSpace::Status status = minimizeOptimumStorageInit(x); do { - status = minimizeOptimumStorageOneStep(x, mode); + status = minimizeOptimumStorageOneStep(x); } while (status==LevenbergMarquardtSpace::Running); return status; } diff --git a/unsupported/test/NonLinearOptimization.cpp b/unsupported/test/NonLinearOptimization.cpp index 1313726a1..7aea7b361 100644 --- a/unsupported/test/NonLinearOptimization.cpp +++ b/unsupported/test/NonLinearOptimization.cpp @@ -317,7 +317,8 @@ void testHybrj() hybrj_functor functor; HybridNonLinearSolver solver(functor); solver.diag.setConstant(n, 1.); - info = solver.solve(x, 2); + solver.useExternalScaling = true; + info = solver.solve(x); // check return value VERIFY( 1 == info); @@ -401,7 +402,8 @@ void testHybrd() solver.parameters.nb_of_subdiagonals = 1; solver.parameters.nb_of_superdiagonals = 1; solver.diag.setConstant(n, 1.); - info = solver.solveNumericalDiff(x, 2); + solver.useExternalScaling = true; + info = solver.solveNumericalDiff(x); // check return value VERIFY( 1 == info); -- cgit v1.2.3