From 784e4f69ec67d511a0620a2b030688f186ebb959 Mon Sep 17 00:00:00 2001
From: Dimitri Pourbaix <dimpbx@apache.org>
Date: Wed, 11 Aug 2010 13:40:39 +0000
Subject: [PATCH] MATH-405 corrected

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@984404 13f79535-47bb-0310-9956-ffa450edef68
---
 .../AbstractLeastSquaresOptimizer.java        |  7 +--
 .../general/LevenbergMarquardtOptimizer.java  | 55 +++++++++++--------
 .../optimization/general/MinpackTest.java     | 11 ++--
 3 files changed, 39 insertions(+), 34 deletions(-)

diff --git a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
index 152c30c3a..c4b19855a 100644
--- a/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
+++ b/src/main/java/org/apache/commons/math/optimization/general/AbstractLeastSquaresOptimizer.java
@@ -247,12 +247,7 @@ public abstract class AbstractLeastSquaresOptimizer implements DifferentiableMul
      * @return chi-square value
      */
     public double getChiSquare() {
-        double chiSquare = 0;
-        for (int i = 0; i < rows; ++i) {
-            final double residual = residuals[i];
-            chiSquare += residual * residual * residualsWeights[i];
-        }
-        return chiSquare;
+        return cost*cost;
     }
 
     /**
diff --git a/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java b/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
index 4cf4ee006..b41456eff 100644
--- a/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
+++ b/src/main/java/org/apache/commons/math/optimization/general/LevenbergMarquardtOptimizer.java
@@ -255,6 +255,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
         double[] diag    = new double[cols];
         double[] oldX    = new double[cols];
         double[] oldRes  = new double[rows];
+        double[] oldObj  = new double[rows];
+        double[] qtf     = new double[rows];
         double[] work1   = new double[cols];
         double[] work2   = new double[cols];
         double[] work3   = new double[cols];
@@ -267,7 +269,9 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
         boolean firstIteration = true;
         VectorialPointValuePair current = new VectorialPointValuePair(point, objective);
         while (true) {
-
+            for (int i=0;i<rows;i++) {
+                qtf[i]=residuals[i];
+            }
             incrementIterationsCounter();
 
             // compute the Q.R. decomposition of the jacobian matrix
@@ -276,8 +280,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             qrDecomposition();
 
             // compute Qt.res
-            qTy(residuals);
-
+            qTy(qtf);
             // now we don't need Q anymore,
             // so let jacobian contain the R matrix with its diagonal elements
             for (int k = 0; k < solvedCols; ++k) {
@@ -315,7 +318,7 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                     if (s != 0) {
                         double sum = 0;
                         for (int i = 0; i <= j; ++i) {
-                            sum += jacobian[i][pj] * residuals[i];
+                            sum += jacobian[i][pj] * qtf[i];
                         }
                         maxCosine = Math.max(maxCosine, Math.abs(sum) / (s * cost));
                     }
@@ -323,6 +326,8 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
             }
             if (maxCosine <= orthoTolerance) {
                 // convergence has been reached
+            	updateResidualsAndCost();
+            	current = new VectorialPointValuePair(point, objective);
                 return current;
             }
 
@@ -343,9 +348,12 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                 double[] tmpVec = residuals;
                 residuals = oldRes;
                 oldRes    = tmpVec;
+                tmpVec    = objective;
+                objective = oldObj;
+                oldObj    = tmpVec;
 
                 // determine the Levenberg-Marquardt parameter
-                determineLMParameter(oldRes, delta, diag, work1, work2, work3);
+                determineLMParameter(qtf, delta, diag, work1, work2, work3);
 
                 // compute the new point and the norm of the evolution direction
                 double lmNorm = 0;
@@ -357,7 +365,6 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                     lmNorm  += s * s;
                 }
                 lmNorm = Math.sqrt(lmNorm);
-
                 // on the first iteration, adjust the initial step bound.
                 if (firstIteration) {
                     delta = Math.min(delta, lmNorm);
@@ -365,7 +372,6 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
 
                 // evaluate the function at x + p and calculate its norm
                 updateResidualsAndCost();
-                current = new VectorialPointValuePair(point, objective);
 
                 // compute the scaled actual reduction
                 double actRed = -1.0;
@@ -421,6 +427,15 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                         xNorm    += xK * xK;
                     }
                     xNorm = Math.sqrt(xNorm);
+                    current = new VectorialPointValuePair(point, objective);
+
+                    // tests for convergence.
+                    if (checker != null) {
+                    // we use the vectorial convergence checker
+                    	if (checker.converged(getIterations(), previous, current)) {
+                    		return current;
+                    	}
+                    }
                 } else {
                     // failed iteration, reset the previous values
                     cost = previousCost;
@@ -431,24 +446,18 @@ public class LevenbergMarquardtOptimizer extends AbstractLeastSquaresOptimizer {
                     tmpVec    = residuals;
                     residuals = oldRes;
                     oldRes    = tmpVec;
+                    tmpVec    = objective;
+                    objective = oldObj;
+                    oldObj    = tmpVec;
                 }
-
-                // tests for convergence.
-                if (checker != null) {
-                    // we use the vectorial convergence checker
-                    if (checker.converged(getIterations(), previous, current)) {
-                        return current;
-                    }
-                } else {
-                    // we use the Levenberg-Marquardt specific convergence parameters
-                    if (((Math.abs(actRed) <= costRelativeTolerance) &&
-                         (preRed <= costRelativeTolerance) &&
-                         (ratio <= 2.0)) ||
-                        (delta <= parRelativeTolerance * xNorm)) {
-                        return current;
-                    }
+                if (checker==null) {
+                	if (((Math.abs(actRed) <= costRelativeTolerance) &&
+                        (preRed <= costRelativeTolerance) &&
+                        (ratio <= 2.0)) ||
+                       (delta <= parRelativeTolerance * xNorm)) {
+                       return current;
+                   }
                 }
-
                 // tests for termination and stringent tolerances
                 // (2.2204e-16 is the machine epsilon for IEEE754)
                 if ((Math.abs(actRed) <= 2.2204e-16) && (preRed <= 2.2204e-16) && (ratio <= 2.0)) {
diff --git a/src/test/java/org/apache/commons/math/optimization/general/MinpackTest.java b/src/test/java/org/apache/commons/math/optimization/general/MinpackTest.java
index 59467b224..620df2d4b 100644
--- a/src/test/java/org/apache/commons/math/optimization/general/MinpackTest.java
+++ b/src/test/java/org/apache/commons/math/optimization/general/MinpackTest.java
@@ -152,14 +152,14 @@ public class MinpackTest extends TestCase {
     minpackTest(new FreudensteinRothFunction(new double[] { 5.0, -20.0 },
                                              12432.833948863, 6.9988751744895,
                                              new double[] {
-                                               11.4121122022341,
-                                               -0.8968550851268697
+                                                11.41300466147456,
+                                                -0.896796038685959
                                              }), false);
     minpackTest(new FreudensteinRothFunction(new double[] { 50.0, -200.0 },
                                              11426454.595762, 6.99887517242903,
                                              new double[] {
-                                               11.412069435091231,
-                                               -0.8968582807605691
+                                                 11.412781785788564,
+                                                 -0.8968051074920405
                                              }), false);
   }
 
@@ -325,7 +325,8 @@ public class MinpackTest extends TestCase {
     minpackTest(new JennrichSampsonFunction(10, new double[] { 0.3, 0.4 },
                                             64.5856498144943, 11.1517793413499,
                                             new double[] {
-                                             0.2578330049, 0.257829976764542
+ //                                            0.2578330049, 0.257829976764542
+                                               0.2578199266368004, 0.25782997676455244
                                             }), false);
   }