From faead6c3fbe343ba97e70f8eda74a488c916b29a Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@apache.org>
Date: Sun, 23 Mar 2008 13:36:03 +0000
Subject: [PATCH] detect numerical problems in Q.R decomposition for
 Levenberg-Marquardt estimator and report them appropriately JIRA: MATH-199

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@640204 13f79535-47bb-0310-9956-ffa450edef68
---
 .../LevenbergMarquardtEstimator.java          |  7 +-
 src/site/xdoc/changes.xml                     |  4 +
 .../LevenbergMarquardtEstimatorTest.java      | 82 +++++++++++++++++++
 3 files changed, 92 insertions(+), 1 deletion(-)

diff --git a/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java b/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java
index 248ebf3c9..07bffbad1 100644
--- a/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java
+++ b/src/java/org/apache/commons/math/estimation/LevenbergMarquardtEstimator.java
@@ -733,8 +733,9 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
    * are performed in non-increasing columns norms order thanks to columns
    * pivoting. The diagonal elements of the R matrix are therefore also in
    * non-increasing absolute values order.</p>
+   * @exception EstimationException if the decomposition cannot be performed
    */
-  private void qrDecomposition() {
+  private void qrDecomposition() throws EstimationException {
 
     // initializations
     for (int k = 0; k < cols; ++k) {
@@ -760,6 +761,10 @@ public class LevenbergMarquardtEstimator extends AbstractEstimator implements Se
           double aki = jacobian[index];
           norm2 += aki * aki;
         }
+        if (Double.isInfinite(norm2) || Double.isNaN(norm2)) {
+            throw new EstimationException("unable to perform Q.R decomposition on the {0}x{1} jacobian matrix",
+                                          new Object[] { new Integer(rows), new Integer(cols) });
+        }
         if (norm2 > ak2) {
           nextColumn = i;
           ak2        = norm2;
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index f5aabb0d1..52c9b8a0d 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -47,6 +47,10 @@ Commons Math Release Notes</title>
       <action dev="luc" type="fix" issue="MATH-198" due-to="Frederick Salardi">
         added an error detection for missing imaginary character while parsing complex string
       </action>
+      <action dev="luc" type="fix" issue="MATH-199" due-to="Mick">
+        detect numerical problems in Q.R decomposition for Levenberg-Marquardt estimator
+        and report them appropriately
+      </action>
     </release>
     <release version="1.2" date="2008-02-24"
     description="This release combines bug fixes and new features. Most notable
diff --git a/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java b/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java
index dc64695a1..b0e6fe43d 100644
--- a/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java
+++ b/src/test/org/apache/commons/math/estimation/LevenbergMarquardtEstimatorTest.java
@@ -20,6 +20,7 @@ package org.apache.commons.math.estimation;
 import java.util.ArrayList;
 import java.util.HashMap;
 import java.util.Iterator;
+import java.util.Locale;
 import java.util.Set;
 
 import org.apache.commons.math.estimation.EstimatedParameter;
@@ -587,6 +588,22 @@ public class LevenbergMarquardtEstimatorTest
     assertEquals( 0.20750021499570379,  circle.getY(),      1.0e-8);
   }
 
+  public void testMath199() {
+      try {
+          QuadraticProblem problem = new QuadraticProblem();
+          problem.addPoint (0, -3.182591015485607, 0.0);
+          problem.addPoint (1, -2.5581184967730577, 4.4E-323);
+          problem.addPoint (2, -2.1488478161387325, 1.0);
+          problem.addPoint (3, -1.9122489313410047, 4.4E-323);
+          problem.addPoint (4, 1.7785661310051026, 0.0);
+          new LevenbergMarquardtEstimator().estimate(problem);
+          fail("an exception should have been thrown");
+      } catch (EstimationException ee) {
+          // expected behavior
+      }
+
+  }
+
   private static class LinearProblem implements EstimationProblem {
 
     public LinearProblem(LinearMeasurement[] measurements) {
@@ -759,6 +776,71 @@ public class LevenbergMarquardtEstimatorTest
     private ArrayList points;
 
   }
+  public class QuadraticProblem extends SimpleEstimationProblem {
+
+      private EstimatedParameter a;
+      private EstimatedParameter b;
+      private EstimatedParameter c;
+
+      public QuadraticProblem() {
+          a = new EstimatedParameter("a", 0.0);
+          b = new EstimatedParameter("b", 0.0);
+          c = new EstimatedParameter("c", 0.0);
+          addParameter(a);
+          addParameter(b);
+          addParameter(c);
+      }
+
+      public void addPoint(double x, double y, double w) {
+          addMeasurement(new LocalMeasurement(x, y, w));
+      }
+
+      public double getA() {
+          return a.getEstimate();
+      }
+
+      public double getB() {
+          return b.getEstimate();
+      }
+
+      public double getC() {
+          return c.getEstimate();
+      }
+
+      public double theoreticalValue(double x) {
+          return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate());
+      }
+
+      private double partial(double x, EstimatedParameter parameter) {
+          if (parameter == a) {
+              return x * x;
+          } else if (parameter == b) {
+              return x;
+          } else {
+              return 1.0;
+          }
+      }
+
+      private class LocalMeasurement extends WeightedMeasurement {
+
+         private final double x;
+
+          // constructor
+          public LocalMeasurement(double x, double y, double w) {
+              super(w, y);
+              this.x = x;
+          }
+
+          public double getTheoreticalValue() {
+              return theoreticalValue(x);
+          }
+
+          public double getPartial(EstimatedParameter parameter) {
+              return partial(x, parameter);
+          }
+
+      }
+  }
 
   public static Test suite() {
     return new TestSuite(LevenbergMarquardtEstimatorTest.class);