Return-Path: X-Original-To: apmail-commons-commits-archive@minotaur.apache.org Delivered-To: apmail-commons-commits-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id 65BB4F1D7 for ; Sat, 6 Apr 2013 23:43:44 +0000 (UTC) Received: (qmail 9657 invoked by uid 500); 6 Apr 2013 23:43:44 -0000 Delivered-To: apmail-commons-commits-archive@commons.apache.org Received: (qmail 9592 invoked by uid 500); 6 Apr 2013 23:43:44 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 9583 invoked by uid 99); 6 Apr 2013 23:43:44 -0000 Received: from nike.apache.org (HELO nike.apache.org) (192.87.106.230) by apache.org (qpsmtpd/0.29) with ESMTP; Sat, 06 Apr 2013 23:43:44 +0000 X-ASF-Spam-Status: No, hits=-2000.0 required=5.0 tests=ALL_TRUSTED X-Spam-Check-By: apache.org Received: from [140.211.11.4] (HELO eris.apache.org) (140.211.11.4) by apache.org (qpsmtpd/0.29) with ESMTP; Sat, 06 Apr 2013 23:43:21 +0000 Received: from eris.apache.org (localhost [127.0.0.1]) by eris.apache.org (Postfix) with ESMTP id 845E92388C48 for ; Sat, 6 Apr 2013 23:42:06 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Subject: svn commit: r857558 [27/39] - in /websites/production/commons/content/proper/commons-math/testapidocs/src-html/org/apache/commons/math3: ./ analysis/ analysis/differentiation/ analysis/interpolation/ complex/ dfp/ distribution/ distribution/fitting/ ex... Date: Sat, 06 Apr 2013 23:42:02 -0000 To: commits@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.8-patched Message-Id: <20130406234206.845E92388C48@eris.apache.org> X-Virus-Checked: Checked by ClamAV on apache.org Modified: websites/production/commons/content/proper/commons-math/testapidocs/src-html/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.html ============================================================================== --- websites/production/commons/content/proper/commons-math/testapidocs/src-html/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.html (original) +++ websites/production/commons/content/proper/commons-math/testapidocs/src-html/org/apache/commons/math3/optim/nonlinear/vector/jacobian/LevenbergMarquardtOptimizerTest.html Sat Apr 6 23:42:01 2013 @@ -26,384 +26,396 @@ 023 import org.apache.commons.math3.optim.PointVectorValuePair; 024 import org.apache.commons.math3.optim.InitialGuess; 025 import org.apache.commons.math3.optim.MaxEval; -026 import org.apache.commons.math3.optim.nonlinear.vector.Target; -027 import org.apache.commons.math3.optim.nonlinear.vector.Weight; -028 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction; -029 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian; -030 import org.apache.commons.math3.analysis.MultivariateVectorFunction; -031 import org.apache.commons.math3.analysis.MultivariateMatrixFunction; -032 import org.apache.commons.math3.exception.ConvergenceException; -033 import org.apache.commons.math3.exception.DimensionMismatchException; -034 import org.apache.commons.math3.exception.TooManyEvaluationsException; -035 import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; -036 import org.apache.commons.math3.linear.SingularMatrixException; -037 import org.apache.commons.math3.util.FastMath; -038 import org.apache.commons.math3.util.Precision; -039 import org.junit.Assert; -040 import org.junit.Test; -041 import org.junit.Ignore; -042 -043 /** -044 * <p>Some of the unit tests are re-implementations of the MINPACK <a -045 * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a -046 * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. -047 * The redistribution policy for MINPACK is available <a -048 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for -049 * convenience, it is reproduced below.</p> -050 -051 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> -052 * <tr><td> -053 * Minpack Copyright Notice (1999) University of Chicago. -054 * All rights reserved -055 * </td></tr> -056 * <tr><td> -057 * Redistribution and use in source and binary forms, with or without -058 * modification, are permitted provided that the following conditions -059 * are met: -060 * <ol> -061 * <li>Redistributions of source code must retain the above copyright -062 * notice, this list of conditions and the following disclaimer.</li> -063 * <li>Redistributions in binary form must reproduce the above -064 * copyright notice, this list of conditions and the following -065 * disclaimer in the documentation and/or other materials provided -066 * with the distribution.</li> -067 * <li>The end-user documentation included with the redistribution, if any, -068 * must include the following acknowledgment: -069 * <code>This product includes software developed by the University of -070 * Chicago, as Operator of Argonne National Laboratory.</code> -071 * Alternately, this acknowledgment may appear in the software itself, -072 * if and wherever such third-party acknowledgments normally appear.</li> -073 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" -074 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE -075 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND -076 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR -077 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES -078 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE -079 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY -080 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR -081 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF -082 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) -083 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION -084 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL -085 * BE CORRECTED.</strong></li> -086 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT -087 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF -088 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, -089 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF -090 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF -091 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER -092 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT -093 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, -094 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE -095 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> -096 * <ol></td></tr> -097 * </table> -098 -099 * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) -100 * @author Burton S. Garbow (original fortran minpack tests) -101 * @author Kenneth E. Hillstrom (original fortran minpack tests) -102 * @author Jorge J. More (original fortran minpack tests) -103 * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) -104 */ -105 public class LevenbergMarquardtOptimizerTest -106 extends AbstractLeastSquaresOptimizerAbstractTest { -107 @Override -108 public AbstractLeastSquaresOptimizer createOptimizer() { -109 return new LevenbergMarquardtOptimizer(); -110 } -111 -112 @Override -113 @Test(expected=SingularMatrixException.class) -114 public void testNonInvertible() { -115 /* -116 * Overrides the method from parent class, since the default singularity -117 * threshold (1e-14) does not trigger the expected exception. -118 */ -119 LinearProblem problem = new LinearProblem(new double[][] { -120 { 1, 2, -3 }, -121 { 2, 1, 3 }, -122 { -3, 0, -9 } -123 }, new double[] { 1, 1, 1 }); -124 -125 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); -126 PointVectorValuePair optimum -127 = optimizer.optimize(new MaxEval(100), -128 problem.getModelFunction(), -129 problem.getModelFunctionJacobian(), -130 problem.getTarget(), -131 new Weight(new double[] { 1, 1, 1 }), -132 new InitialGuess(new double[] { 0, 0, 0 })); -133 Assert.assertTrue(FastMath.sqrt(optimizer.getTargetSize()) * optimizer.getRMS() > 0.6); -134 -135 optimizer.computeCovariances(optimum.getPoint(), 1.5e-14); -136 } -137 -138 @Test -139 public void testControlParameters() { -140 CircleVectorial circle = new CircleVectorial(); -141 circle.addPoint( 30.0, 68.0); -142 circle.addPoint( 50.0, -6.0); -143 circle.addPoint(110.0, -20.0); -144 circle.addPoint( 35.0, 15.0); -145 circle.addPoint( 45.0, 97.0); -146 checkEstimate(circle.getModelFunction(), -147 circle.getModelFunctionJacobian(), -148 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false); -149 checkEstimate(circle.getModelFunction(), -150 circle.getModelFunctionJacobian(), -151 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true); -152 checkEstimate(circle.getModelFunction(), -153 circle.getModelFunctionJacobian(), -154 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true); -155 circle.addPoint(300, -300); -156 checkEstimate(circle.getModelFunction(), -157 circle.getModelFunctionJacobian(), -158 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true); -159 } -160 -161 private void checkEstimate(ModelFunction problem, -162 ModelFunctionJacobian problemJacobian, -163 double initialStepBoundFactor, int maxCostEval, -164 double costRelativeTolerance, double parRelativeTolerance, -165 double orthoTolerance, boolean shouldFail) { -166 try { -167 LevenbergMarquardtOptimizer optimizer -168 = new LevenbergMarquardtOptimizer(initialStepBoundFactor, -169 costRelativeTolerance, -170 parRelativeTolerance, -171 orthoTolerance, -172 Precision.SAFE_MIN); -173 optimizer.optimize(new MaxEval(maxCostEval), -174 problem, -175 problemJacobian, -176 new Target(new double[] { 0, 0, 0, 0, 0 }), -177 new Weight(new double[] { 1, 1, 1, 1, 1 }), -178 new InitialGuess(new double[] { 98.680, 47.345 })); -179 Assert.assertTrue(!shouldFail); -180 } catch (DimensionMismatchException ee) { -181 Assert.assertTrue(shouldFail); -182 } catch (TooManyEvaluationsException ee) { -183 Assert.assertTrue(shouldFail); -184 } -185 } -186 -187 /** -188 * Non-linear test case: fitting of decay curve (from Chapter 8 of -189 * Bevington's textbook, "Data reduction and analysis for the physical sciences"). -190 * XXX The expected ("reference") values may not be accurate and the tolerance too -191 * relaxed for this test to be currently really useful (the issue is under -192 * investigation). -193 */ -194 @Test -195 public void testBevington() { -196 final double[][] dataPoints = { -197 // column 1 = times -198 { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, -199 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, -200 315, 330, 345, 360, 375, 390, 405, 420, 435, 450, -201 465, 480, 495, 510, 525, 540, 555, 570, 585, 600, -202 615, 630, 645, 660, 675, 690, 705, 720, 735, 750, -203 765, 780, 795, 810, 825, 840, 855, 870, 885, }, -204 // column 2 = measured counts -205 { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89, -206 74, 61, 66, 68, 48, 54, 51, 46, 55, 29, -207 28, 37, 49, 26, 35, 29, 31, 24, 25, 35, -208 24, 30, 26, 28, 21, 18, 20, 27, 17, 17, -209 14, 17, 24, 11, 22, 17, 12, 10, 13, 16, -210 9, 9, 14, 21, 17, 13, 12, 18, 10, }, -211 }; -212 -213 final BevingtonProblem problem = new BevingtonProblem(); -214 -215 final int len = dataPoints[0].length; -216 final double[] weights = new double[len]; -217 for (int i = 0; i < len; i++) { -218 problem.addPoint(dataPoints[0][i], -219 dataPoints[1][i]); -220 -221 weights[i] = 1 / dataPoints[1][i]; -222 } -223 -224 final LevenbergMarquardtOptimizer optimizer -225 = new LevenbergMarquardtOptimizer(); +026 import org.apache.commons.math3.optim.SimpleBounds; +027 import org.apache.commons.math3.optim.nonlinear.vector.Target; +028 import org.apache.commons.math3.optim.nonlinear.vector.Weight; +029 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction; +030 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian; +031 import org.apache.commons.math3.analysis.MultivariateVectorFunction; +032 import org.apache.commons.math3.analysis.MultivariateMatrixFunction; +033 import org.apache.commons.math3.exception.ConvergenceException; +034 import org.apache.commons.math3.exception.DimensionMismatchException; +035 import org.apache.commons.math3.exception.TooManyEvaluationsException; +036 import org.apache.commons.math3.exception.MathUnsupportedOperationException; +037 import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; +038 import org.apache.commons.math3.linear.SingularMatrixException; +039 import org.apache.commons.math3.util.FastMath; +040 import org.apache.commons.math3.util.Precision; +041 import org.junit.Assert; +042 import org.junit.Test; +043 import org.junit.Ignore; +044 +045 /** +046 * <p>Some of the unit tests are re-implementations of the MINPACK <a +047 * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a +048 * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. +049 * The redistribution policy for MINPACK is available <a +050 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for +051 * convenience, it is reproduced below.</p> +052 +053 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> +054 * <tr><td> +055 * Minpack Copyright Notice (1999) University of Chicago. +056 * All rights reserved +057 * </td></tr> +058 * <tr><td> +059 * Redistribution and use in source and binary forms, with or without +060 * modification, are permitted provided that the following conditions +061 * are met: +062 * <ol> +063 * <li>Redistributions of source code must retain the above copyright +064 * notice, this list of conditions and the following disclaimer.</li> +065 * <li>Redistributions in binary form must reproduce the above +066 * copyright notice, this list of conditions and the following +067 * disclaimer in the documentation and/or other materials provided +068 * with the distribution.</li> +069 * <li>The end-user documentation included with the redistribution, if any, +070 * must include the following acknowledgment: +071 * <code>This product includes software developed by the University of +072 * Chicago, as Operator of Argonne National Laboratory.</code> +073 * Alternately, this acknowledgment may appear in the software itself, +074 * if and wherever such third-party acknowledgments normally appear.</li> +075 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" +076 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE +077 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND +078 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR +079 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES +080 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE +081 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY +082 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR +083 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF +084 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) +085 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION +086 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL +087 * BE CORRECTED.</strong></li> +088 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT +089 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF +090 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, +091 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF +092 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF +093 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER +094 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT +095 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, +096 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE +097 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> +098 * <ol></td></tr> +099 * </table> +100 +101 * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) +102 * @author Burton S. Garbow (original fortran minpack tests) +103 * @author Kenneth E. Hillstrom (original fortran minpack tests) +104 * @author Jorge J. More (original fortran minpack tests) +105 * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) +106 */ +107 public class LevenbergMarquardtOptimizerTest +108 extends AbstractLeastSquaresOptimizerAbstractTest { +109 @Override +110 public AbstractLeastSquaresOptimizer createOptimizer() { +111 return new LevenbergMarquardtOptimizer(); +112 } +113 +114 @Test(expected=MathUnsupportedOperationException.class) +115 public void testConstraintsUnsupported() { +116 createOptimizer().optimize(new MaxEval(100), +117 new Target(new double[] { 2 }), +118 new Weight(new double[] { 1 }), +119 new InitialGuess(new double[] { 1, 2 }), +120 new SimpleBounds(new double[] { -10, 0 }, +121 new double[] { 20, 30 })); +122 } +123 +124 @Override +125 @Test(expected=SingularMatrixException.class) +126 public void testNonInvertible() { +127 /* +128 * Overrides the method from parent class, since the default singularity +129 * threshold (1e-14) does not trigger the expected exception. +130 */ +131 LinearProblem problem = new LinearProblem(new double[][] { +132 { 1, 2, -3 }, +133 { 2, 1, 3 }, +134 { -3, 0, -9 } +135 }, new double[] { 1, 1, 1 }); +136 +137 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); +138 PointVectorValuePair optimum +139 = optimizer.optimize(new MaxEval(100), +140 problem.getModelFunction(), +141 problem.getModelFunctionJacobian(), +142 problem.getTarget(), +143 new Weight(new double[] { 1, 1, 1 }), +144 new InitialGuess(new double[] { 0, 0, 0 })); +145 Assert.assertTrue(FastMath.sqrt(optimizer.getTargetSize()) * optimizer.getRMS() > 0.6); +146 +147 optimizer.computeCovariances(optimum.getPoint(), 1.5e-14); +148 } +149 +150 @Test +151 public void testControlParameters() { +152 CircleVectorial circle = new CircleVectorial(); +153 circle.addPoint( 30.0, 68.0); +154 circle.addPoint( 50.0, -6.0); +155 circle.addPoint(110.0, -20.0); +156 circle.addPoint( 35.0, 15.0); +157 circle.addPoint( 45.0, 97.0); +158 checkEstimate(circle.getModelFunction(), +159 circle.getModelFunctionJacobian(), +160 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false); +161 checkEstimate(circle.getModelFunction(), +162 circle.getModelFunctionJacobian(), +163 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true); +164 checkEstimate(circle.getModelFunction(), +165 circle.getModelFunctionJacobian(), +166 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true); +167 circle.addPoint(300, -300); +168 checkEstimate(circle.getModelFunction(), +169 circle.getModelFunctionJacobian(), +170 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true); +171 } +172 +173 private void checkEstimate(ModelFunction problem, +174 ModelFunctionJacobian problemJacobian, +175 double initialStepBoundFactor, int maxCostEval, +176 double costRelativeTolerance, double parRelativeTolerance, +177 double orthoTolerance, boolean shouldFail) { +178 try { +179 LevenbergMarquardtOptimizer optimizer +180 = new LevenbergMarquardtOptimizer(initialStepBoundFactor, +181 costRelativeTolerance, +182 parRelativeTolerance, +183 orthoTolerance, +184 Precision.SAFE_MIN); +185 optimizer.optimize(new MaxEval(maxCostEval), +186 problem, +187 problemJacobian, +188 new Target(new double[] { 0, 0, 0, 0, 0 }), +189 new Weight(new double[] { 1, 1, 1, 1, 1 }), +190 new InitialGuess(new double[] { 98.680, 47.345 })); +191 Assert.assertTrue(!shouldFail); +192 } catch (DimensionMismatchException ee) { +193 Assert.assertTrue(shouldFail); +194 } catch (TooManyEvaluationsException ee) { +195 Assert.assertTrue(shouldFail); +196 } +197 } +198 +199 /** +200 * Non-linear test case: fitting of decay curve (from Chapter 8 of +201 * Bevington's textbook, "Data reduction and analysis for the physical sciences"). +202 * XXX The expected ("reference") values may not be accurate and the tolerance too +203 * relaxed for this test to be currently really useful (the issue is under +204 * investigation). +205 */ +206 @Test +207 public void testBevington() { +208 final double[][] dataPoints = { +209 // column 1 = times +210 { 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, +211 165, 180, 195, 210, 225, 240, 255, 270, 285, 300, +212 315, 330, 345, 360, 375, 390, 405, 420, 435, 450, +213 465, 480, 495, 510, 525, 540, 555, 570, 585, 600, +214 615, 630, 645, 660, 675, 690, 705, 720, 735, 750, +215 765, 780, 795, 810, 825, 840, 855, 870, 885, }, +216 // column 2 = measured counts +217 { 775, 479, 380, 302, 185, 157, 137, 119, 110, 89, +218 74, 61, 66, 68, 48, 54, 51, 46, 55, 29, +219 28, 37, 49, 26, 35, 29, 31, 24, 25, 35, +220 24, 30, 26, 28, 21, 18, 20, 27, 17, 17, +221 14, 17, 24, 11, 22, 17, 12, 10, 13, 16, +222 9, 9, 14, 21, 17, 13, 12, 18, 10, }, +223 }; +224 +225 final BevingtonProblem problem = new BevingtonProblem(); 226 -227 final PointVectorValuePair optimum -228 = optimizer.optimize(new MaxEval(100), -229 problem.getModelFunction(), -230 problem.getModelFunctionJacobian(), -231 new Target(dataPoints[1]), -232 new Weight(weights), -233 new InitialGuess(new double[] { 10, 900, 80, 27, 225 })); -234 -235 final double[] solution = optimum.getPoint(); -236 final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 }; -237 -238 final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14); -239 final double[][] expectedCovarMatrix = { -240 { 3.38, -3.69, 27.98, -2.34, -49.24 }, -241 { -3.69, 2492.26, 81.89, -69.21, -8.9 }, -242 { 27.98, 81.89, 468.99, -44.22, -615.44 }, -243 { -2.34, -69.21, -44.22, 6.39, 53.80 }, -244 { -49.24, -8.9, -615.44, 53.8, 929.45 } -245 }; +227 final int len = dataPoints[0].length; +228 final double[] weights = new double[len]; +229 for (int i = 0; i < len; i++) { +230 problem.addPoint(dataPoints[0][i], +231 dataPoints[1][i]); +232 +233 weights[i] = 1 / dataPoints[1][i]; +234 } +235 +236 final LevenbergMarquardtOptimizer optimizer +237 = new LevenbergMarquardtOptimizer(); +238 +239 final PointVectorValuePair optimum +240 = optimizer.optimize(new MaxEval(100), +241 problem.getModelFunction(), +242 problem.getModelFunctionJacobian(), +243 new Target(dataPoints[1]), +244 new Weight(weights), +245 new InitialGuess(new double[] { 10, 900, 80, 27, 225 })); 246 -247 final int numParams = expectedSolution.length; -248 -249 // Check that the computed solution is within the reference error range. -250 for (int i = 0; i < numParams; i++) { -251 final double error = FastMath.sqrt(expectedCovarMatrix[i][i]); -252 Assert.assertEquals("Parameter " + i, expectedSolution[i], solution[i], error); -253 } -254 -255 // Check that each entry of the computed covariance matrix is within 10% -256 // of the reference matrix entry. -257 for (int i = 0; i < numParams; i++) { -258 for (int j = 0; j < numParams; j++) { -259 Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]", -260 expectedCovarMatrix[i][j], -261 covarMatrix[i][j], -262 FastMath.abs(0.1 * expectedCovarMatrix[i][j])); -263 } -264 } -265 } +247 final double[] solution = optimum.getPoint(); +248 final double[] expectedSolution = { 10.4, 958.3, 131.4, 33.9, 205.0 }; +249 +250 final double[][] covarMatrix = optimizer.computeCovariances(solution, 1e-14); +251 final double[][] expectedCovarMatrix = { +252 { 3.38, -3.69, 27.98, -2.34, -49.24 }, +253 { -3.69, 2492.26, 81.89, -69.21, -8.9 }, +254 { 27.98, 81.89, 468.99, -44.22, -615.44 }, +255 { -2.34, -69.21, -44.22, 6.39, 53.80 }, +256 { -49.24, -8.9, -615.44, 53.8, 929.45 } +257 }; +258 +259 final int numParams = expectedSolution.length; +260 +261 // Check that the computed solution is within the reference error range. +262 for (int i = 0; i < numParams; i++) { +263 final double error = FastMath.sqrt(expectedCovarMatrix[i][i]); +264 Assert.assertEquals("Parameter " + i, expectedSolution[i], solution[i], error); +265 } 266 -267 @Test -268 public void testCircleFitting2() { -269 final double xCenter = 123.456; -270 final double yCenter = 654.321; -271 final double xSigma = 10; -272 final double ySigma = 15; -273 final double radius = 111.111; -274 // The test is extremely sensitive to the seed. -275 final long seed = 59421061L; -276 final RandomCirclePointGenerator factory -277 = new RandomCirclePointGenerator(xCenter, yCenter, radius, -278 xSigma, ySigma, -279 seed); -280 final CircleProblem circle = new CircleProblem(xSigma, ySigma); -281 -282 final int numPoints = 10; -283 for (Vector2D p : factory.generate(numPoints)) { -284 circle.addPoint(p.getX(), p.getY()); -285 } -286 -287 // First guess for the center's coordinates and radius. -288 final double[] init = { 90, 659, 115 }; -289 -290 final LevenbergMarquardtOptimizer optimizer -291 = new LevenbergMarquardtOptimizer(); -292 final PointVectorValuePair optimum = optimizer.optimize(new MaxEval(100), -293 circle.getModelFunction(), -294 circle.getModelFunctionJacobian(), -295 new Target(circle.target()), -296 new Weight(circle.weight()), -297 new InitialGuess(init)); +267 // Check that each entry of the computed covariance matrix is within 10% +268 // of the reference matrix entry. +269 for (int i = 0; i < numParams; i++) { +270 for (int j = 0; j < numParams; j++) { +271 Assert.assertEquals("Covariance matrix [" + i + "][" + j + "]", +272 expectedCovarMatrix[i][j], +273 covarMatrix[i][j], +274 FastMath.abs(0.1 * expectedCovarMatrix[i][j])); +275 } +276 } +277 } +278 +279 @Test +280 public void testCircleFitting2() { +281 final double xCenter = 123.456; +282 final double yCenter = 654.321; +283 final double xSigma = 10; +284 final double ySigma = 15; +285 final double radius = 111.111; +286 // The test is extremely sensitive to the seed. +287 final long seed = 59421061L; +288 final RandomCirclePointGenerator factory +289 = new RandomCirclePointGenerator(xCenter, yCenter, radius, +290 xSigma, ySigma, +291 seed); +292 final CircleProblem circle = new CircleProblem(xSigma, ySigma); +293 +294 final int numPoints = 10; +295 for (Vector2D p : factory.generate(numPoints)) { +296 circle.addPoint(p.getX(), p.getY()); +297 } 298 -299 final double[] paramFound = optimum.getPoint(); -300 -301 // Retrieve errors estimation. -302 final double[] asymptoticStandardErrorFound = optimizer.computeSigma(paramFound, 1e-14); -303 -304 // Check that the parameters are found within the assumed error bars. -305 Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]); -306 Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]); -307 Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]); -308 } -309 -310 private static class QuadraticProblem { -311 private List<Double> x; -312 private List<Double> y; -313 -314 public QuadraticProblem() { -315 x = new ArrayList<Double>(); -316 y = new ArrayList<Double>(); -317 } -318 -319 public void addPoint(double x, double y) { -320 this.x.add(x); -321 this.y.add(y); -322 } -323 -324 public ModelFunction getModelFunction() { -325 return new ModelFunction(new MultivariateVectorFunction() { -326 public double[] value(double[] variables) { -327 double[] values = new double[x.size()]; -328 for (int i = 0; i < values.length; ++i) { -329 values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2]; -330 } -331 return values; -332 } -333 }); +299 // First guess for the center's coordinates and radius. +300 final double[] init = { 90, 659, 115 }; +301 +302 final LevenbergMarquardtOptimizer optimizer +303 = new LevenbergMarquardtOptimizer(); +304 final PointVectorValuePair optimum = optimizer.optimize(new MaxEval(100), +305 circle.getModelFunction(), +306 circle.getModelFunctionJacobian(), +307 new Target(circle.target()), +308 new Weight(circle.weight()), +309 new InitialGuess(init)); +310 +311 final double[] paramFound = optimum.getPoint(); +312 +313 // Retrieve errors estimation. +314 final double[] asymptoticStandardErrorFound = optimizer.computeSigma(paramFound, 1e-14); +315 +316 // Check that the parameters are found within the assumed error bars. +317 Assert.assertEquals(xCenter, paramFound[0], asymptoticStandardErrorFound[0]); +318 Assert.assertEquals(yCenter, paramFound[1], asymptoticStandardErrorFound[1]); +319 Assert.assertEquals(radius, paramFound[2], asymptoticStandardErrorFound[2]); +320 } +321 +322 private static class QuadraticProblem { +323 private List<Double> x; +324 private List<Double> y; +325 +326 public QuadraticProblem() { +327 x = new ArrayList<Double>(); +328 y = new ArrayList<Double>(); +329 } +330 +331 public void addPoint(double x, double y) { +332 this.x.add(x); +333 this.y.add(y); 334 } 335 -336 public ModelFunctionJacobian getModelFunctionJacobian() { -337 return new ModelFunctionJacobian(new MultivariateMatrixFunction() { -338 public double[][] value(double[] params) { -339 double[][] jacobian = new double[x.size()][3]; -340 for (int i = 0; i < jacobian.length; ++i) { -341 jacobian[i][0] = x.get(i) * x.get(i); -342 jacobian[i][1] = x.get(i); -343 jacobian[i][2] = 1.0; -344 } -345 return jacobian; -346 } -347 }); -348 } -349 } -350 -351 private static class BevingtonProblem { -352 private List<Double> time; -353 private List<Double> count; -354 -355 public BevingtonProblem() { -356 time = new ArrayList<Double>(); -357 count = new ArrayList<Double>(); -358 } -359 -360 public void addPoint(double t, double c) { -361 time.add(t); -362 count.add(c); -363 } -364 -365 public ModelFunction getModelFunction() { -366 return new ModelFunction(new MultivariateVectorFunction() { -367 public double[] value(double[] params) { -368 double[] values = new double[time.size()]; -369 for (int i = 0; i < values.length; ++i) { -370 final double t = time.get(i); -371 values[i] = params[0] + -372 params[1] * Math.exp(-t / params[3]) + -373 params[2] * Math.exp(-t / params[4]); -374 } -375 return values; -376 } -377 }); -378 } -379 -380 public ModelFunctionJacobian getModelFunctionJacobian() { -381 return new ModelFunctionJacobian(new MultivariateMatrixFunction() { -382 public double[][] value(double[] params) { -383 double[][] jacobian = new double[time.size()][5]; -384 -385 for (int i = 0; i < jacobian.length; ++i) { -386 final double t = time.get(i); -387 jacobian[i][0] = 1; -388 -389 final double p3 = params[3]; -390 final double p4 = params[4]; -391 final double tOp3 = t / p3; -392 final double tOp4 = t / p4; -393 jacobian[i][1] = Math.exp(-tOp3); -394 jacobian[i][2] = Math.exp(-tOp4); -395 jacobian[i][3] = params[1] * Math.exp(-tOp3) * tOp3 / p3; -396 jacobian[i][4] = params[2] * Math.exp(-tOp4) * tOp4 / p4; -397 } -398 return jacobian; -399 } -400 }); -401 } -402 } -403 } +336 public ModelFunction getModelFunction() { +337 return new ModelFunction(new MultivariateVectorFunction() { +338 public double[] value(double[] variables) { +339 double[] values = new double[x.size()]; +340 for (int i = 0; i < values.length; ++i) { +341 values[i] = (variables[0] * x.get(i) + variables[1]) * x.get(i) + variables[2]; +342 } +343 return values; +344 } +345 }); +346 } +347 +348 public ModelFunctionJacobian getModelFunctionJacobian() { +349 return new ModelFunctionJacobian(new MultivariateMatrixFunction() { +350 public double[][] value(double[] params) { +351 double[][] jacobian = new double[x.size()][3]; +352 for (int i = 0; i < jacobian.length; ++i) { +353 jacobian[i][0] = x.get(i) * x.get(i); +354 jacobian[i][1] = x.get(i); +355 jacobian[i][2] = 1.0; +356 } +357 return jacobian; +358 } +359 }); +360 } +361 } +362 +363 private static class BevingtonProblem { +364 private List<Double> time; +365 private List<Double> count; +366 +367 public BevingtonProblem() { +368 time = new ArrayList<Double>(); +369 count = new ArrayList<Double>(); +370 } +371 +372 public void addPoint(double t, double c) { +373 time.add(t); +374 count.add(c); +375 } +376 +377 public ModelFunction getModelFunction() { +378 return new ModelFunction(new MultivariateVectorFunction() { +379 public double[] value(double[] params) { +380 double[] values = new double[time.size()]; +381 for (int i = 0; i < values.length; ++i) { +382 final double t = time.get(i); +383 values[i] = params[0] + +384 params[1] * Math.exp(-t / params[3]) + +385 params[2] * Math.exp(-t / params[4]); +386 } +387 return values; +388 } +389 }); +390 } +391 +392 public ModelFunctionJacobian getModelFunctionJacobian() { +393 return new ModelFunctionJacobian(new MultivariateMatrixFunction() { +394 public double[][] value(double[] params) { +395 double[][] jacobian = new double[time.size()][5]; +396 +397 for (int i = 0; i < jacobian.length; ++i) { +398 final double t = time.get(i); +399 jacobian[i][0] = 1; +400 +401 final double p3 = params[3]; +402 final double p4 = params[4]; +403 final double tOp3 = t / p3; +404 final double tOp4 = t / p4; +405 jacobian[i][1] = Math.exp(-tOp3); +406 jacobian[i][2] = Math.exp(-tOp4); +407 jacobian[i][3] = params[1] * Math.exp(-tOp3) * tOp3 / p3; +408 jacobian[i][4] = params[2] * Math.exp(-tOp4) * tOp4 / p4; +409 } +410 return jacobian; +411 } +412 }); +413 } +414 } +415 }