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

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<Scalar>::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<Scalar>::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<FunctorType,Scalar>::hybrj1(
     parameters.maxfev = 100*(n+1);
     parameters.xtol = tol;
     diag.setConstant(n, 1.);
-    return solve(
-        x,
-        2
-    );
+    useExternalScaling = true;
+    return solve(x);
 }
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveInit(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveInit(FVectorType  &x)
 {
     n = x.size();
 
@@ -179,9 +157,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::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<FunctorType,Scalar>::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<FunctorType,Scalar>::solveInit(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(FVectorType  &x)
 {
     int j;
     std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@@ -231,10 +206,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::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<FunctorType,Scalar>::solveOneStep(
     qtf = fjac.transpose() * fvec;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     while (true) {
@@ -372,14 +347,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveOneStep(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solve(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::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<FunctorType,Scalar>::hybrd1(
     parameters.xtol = tol;
 
     diag.setConstant(n, 1.);
-    return solveNumericalDiff(
-        x,
-        2
-    );
+    useExternalScaling = true;
+    return solveNumericalDiff(x);
 }
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffInit(FVectorType  &x)
 {
     n = x.size();
 
@@ -425,10 +392,9 @@ HybridNonLinearSolver<FunctorType,Scalar>::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<FunctorType,Scalar>::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<FunctorType,Scalar>::solveNumericalDiffInit(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(FVectorType  &x)
 {
     int j;
     std::vector<PlanarRotation<Scalar> > v_givens(n), w_givens(n);
@@ -480,10 +443,10 @@ HybridNonLinearSolver<FunctorType,Scalar>::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<FunctorType,Scalar>::solveNumericalDiffOneStep(
     qtf = fjac.transpose() * fvec;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     while (true) {
@@ -621,14 +584,11 @@ HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiffOneStep(
 
 template<typename FunctorType, typename Scalar>
 HybridNonLinearSolverSpace::Status
-HybridNonLinearSolver<FunctorType,Scalar>::solveNumericalDiff(
-        FVectorType  &x,
-        const int mode
-        )
+HybridNonLinearSolver<FunctorType,Scalar>::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<Scalar>::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<Scalar>::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<FunctorType,Scalar>::lmder1(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimize(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::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<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeInit(FVectorType  &x)
 {
     n = x.size();
     m = functor.values();
@@ -201,9 +178,9 @@ LevenbergMarquardt<FunctorType,Scalar>::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<FunctorType,Scalar>::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<FunctorType,Scalar>::minimizeInit(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOneStep(FVectorType  &x)
 {
     int j;
 
@@ -257,10 +231,10 @@ LevenbergMarquardt<FunctorType,Scalar>::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<FunctorType,Scalar>::minimizeOneStep(
         return LevenbergMarquardtSpace::CosinusTooSmall;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     do {
@@ -406,10 +380,7 @@ LevenbergMarquardt<FunctorType,Scalar>::lmstr1(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageInit(FVectorType  &x)
 {
     n = x.size();
     m = functor.values();
@@ -423,9 +394,9 @@ LevenbergMarquardt<FunctorType,Scalar>::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<FunctorType,Scalar>::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<FunctorType,Scalar>::minimizeOptimumStorageInit(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(FVectorType  &x)
 {
     int i, j;
     bool sing;
@@ -514,10 +482,10 @@ LevenbergMarquardt<FunctorType,Scalar>::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<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
         return LevenbergMarquardtSpace::CosinusTooSmall;
 
     /* rescale if necessary. */
-    if (mode != 2)
+    if (!useExternalScaling)
         diag = diag.cwiseMax(wa2);
 
     do {
@@ -635,14 +603,11 @@ LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorageOneStep(
 
 template<typename FunctorType, typename Scalar>
 LevenbergMarquardtSpace::Status
-LevenbergMarquardt<FunctorType,Scalar>::minimizeOptimumStorage(
-        FVectorType  &x,
-        const int mode
-        )
+LevenbergMarquardt<FunctorType,Scalar>::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<hybrj_functor> 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);