001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.estimation; 019 020 import java.util.ArrayList; 021 import java.util.HashSet; 022 023 import junit.framework.Test; 024 import junit.framework.TestCase; 025 import junit.framework.TestSuite; 026 027 /** 028 * <p>Some of the unit tests are re-implementations of the MINPACK <a 029 * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a 030 * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. 031 * The redistribution policy for MINPACK is available <a 032 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for 033 * convenience, it is reproduced below.</p> 034 035 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> 036 * <tr><td> 037 * Minpack Copyright Notice (1999) University of Chicago. 038 * All rights reserved 039 * </td></tr> 040 * <tr><td> 041 * Redistribution and use in source and binary forms, with or without 042 * modification, are permitted provided that the following conditions 043 * are met: 044 * <ol> 045 * <li>Redistributions of source code must retain the above copyright 046 * notice, this list of conditions and the following disclaimer.</li> 047 * <li>Redistributions in binary form must reproduce the above 048 * copyright notice, this list of conditions and the following 049 * disclaimer in the documentation and/or other materials provided 050 * with the distribution.</li> 051 * <li>The end-user documentation included with the redistribution, if any, 052 * must include the following acknowledgment: 053 * <code>This product includes software developed by the University of 054 * Chicago, as Operator of Argonne National Laboratory.</code> 055 * Alternately, this acknowledgment may appear in the software itself, 056 * if and wherever such third-party acknowledgments normally appear.</li> 057 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" 058 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE 059 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND 060 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR 061 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES 062 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE 063 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY 064 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR 065 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF 066 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) 067 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION 068 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL 069 * BE CORRECTED.</strong></li> 070 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT 071 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF 072 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, 073 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF 074 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF 075 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER 076 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT 077 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, 078 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE 079 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> 080 * <ol></td></tr> 081 * </table> 082 083 * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) 084 * @author Burton S. Garbow (original fortran minpack tests) 085 * @author Kenneth E. Hillstrom (original fortran minpack tests) 086 * @author Jorge J. More (original fortran minpack tests) 087 * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) 088 */ 089 @Deprecated 090 public class LevenbergMarquardtEstimatorTest 091 extends TestCase { 092 093 public LevenbergMarquardtEstimatorTest(String name) { 094 super(name); 095 } 096 097 public void testTrivial() throws EstimationException { 098 LinearProblem problem = 099 new LinearProblem(new LinearMeasurement[] { 100 new LinearMeasurement(new double[] {2}, 101 new EstimatedParameter[] { 102 new EstimatedParameter("p0", 0) 103 }, 3.0) 104 }); 105 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 106 estimator.estimate(problem); 107 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 108 try { 109 estimator.guessParametersErrors(problem); 110 fail("an exception should have been thrown"); 111 } catch (EstimationException ee) { 112 // expected behavior 113 } catch (Exception e) { 114 fail("wrong exception caught"); 115 } 116 assertEquals(1.5, 117 problem.getUnboundParameters()[0].getEstimate(), 118 1.0e-10); 119 } 120 121 public void testQRColumnsPermutation() throws EstimationException { 122 123 EstimatedParameter[] x = { 124 new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0) 125 }; 126 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 127 new LinearMeasurement(new double[] { 1.0, -1.0 }, 128 new EstimatedParameter[] { x[0], x[1] }, 129 4.0), 130 new LinearMeasurement(new double[] { 2.0 }, 131 new EstimatedParameter[] { x[1] }, 132 6.0), 133 new LinearMeasurement(new double[] { 1.0, -2.0 }, 134 new EstimatedParameter[] { x[0], x[1] }, 135 1.0) 136 }); 137 138 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 139 estimator.estimate(problem); 140 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 141 assertEquals(7.0, x[0].getEstimate(), 1.0e-10); 142 assertEquals(3.0, x[1].getEstimate(), 1.0e-10); 143 144 } 145 146 public void testNoDependency() throws EstimationException { 147 EstimatedParameter[] p = new EstimatedParameter[] { 148 new EstimatedParameter("p0", 0), 149 new EstimatedParameter("p1", 0), 150 new EstimatedParameter("p2", 0), 151 new EstimatedParameter("p3", 0), 152 new EstimatedParameter("p4", 0), 153 new EstimatedParameter("p5", 0) 154 }; 155 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 156 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0), 157 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1), 158 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2), 159 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3), 160 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4), 161 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5) 162 }); 163 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 164 estimator.estimate(problem); 165 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 166 for (int i = 0; i < p.length; ++i) { 167 assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10); 168 } 169 } 170 171 public void testOneSet() throws EstimationException { 172 173 EstimatedParameter[] p = { 174 new EstimatedParameter("p0", 0), 175 new EstimatedParameter("p1", 0), 176 new EstimatedParameter("p2", 0) 177 }; 178 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 179 new LinearMeasurement(new double[] { 1.0 }, 180 new EstimatedParameter[] { p[0] }, 181 1.0), 182 new LinearMeasurement(new double[] { -1.0, 1.0 }, 183 new EstimatedParameter[] { p[0], p[1] }, 184 1.0), 185 new LinearMeasurement(new double[] { -1.0, 1.0 }, 186 new EstimatedParameter[] { p[1], p[2] }, 187 1.0) 188 }); 189 190 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 191 estimator.estimate(problem); 192 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 193 assertEquals(1.0, p[0].getEstimate(), 1.0e-10); 194 assertEquals(2.0, p[1].getEstimate(), 1.0e-10); 195 assertEquals(3.0, p[2].getEstimate(), 1.0e-10); 196 197 } 198 199 public void testTwoSets() throws EstimationException { 200 EstimatedParameter[] p = { 201 new EstimatedParameter("p0", 0), 202 new EstimatedParameter("p1", 1), 203 new EstimatedParameter("p2", 2), 204 new EstimatedParameter("p3", 3), 205 new EstimatedParameter("p4", 4), 206 new EstimatedParameter("p5", 5) 207 }; 208 209 double epsilon = 1.0e-7; 210 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 211 212 // 4 elements sub-problem 213 new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 }, 214 new EstimatedParameter[] { p[0], p[1], p[3] }, 215 2.0), 216 new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 }, 217 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 218 -9.0), 219 new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 }, 220 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 221 2.0), 222 new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 }, 223 new EstimatedParameter[] { p[1], p[2], p[3] }, 224 2.0), 225 226 // 2 elements sub-problem 227 new LinearMeasurement(new double[] { epsilon, 1.0 }, 228 new EstimatedParameter[] { p[4], p[5] }, 229 1.0 + epsilon * epsilon), 230 new LinearMeasurement(new double[] { 1.0, 1.0 }, 231 new EstimatedParameter[] { p[4], p[5] }, 232 2.0) 233 234 }); 235 236 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 237 estimator.estimate(problem); 238 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 239 assertEquals( 3.0, p[0].getEstimate(), 1.0e-10); 240 assertEquals( 4.0, p[1].getEstimate(), 1.0e-10); 241 assertEquals(-1.0, p[2].getEstimate(), 1.0e-10); 242 assertEquals(-2.0, p[3].getEstimate(), 1.0e-10); 243 assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10); 244 assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10); 245 246 } 247 248 public void testNonInversible() throws EstimationException { 249 250 EstimatedParameter[] p = { 251 new EstimatedParameter("p0", 0), 252 new EstimatedParameter("p1", 0), 253 new EstimatedParameter("p2", 0) 254 }; 255 LinearMeasurement[] m = new LinearMeasurement[] { 256 new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 }, 257 new EstimatedParameter[] { p[0], p[1], p[2] }, 258 1.0), 259 new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 }, 260 new EstimatedParameter[] { p[0], p[1], p[2] }, 261 1.0), 262 new LinearMeasurement(new double[] { -3.0, -9.0 }, 263 new EstimatedParameter[] { p[0], p[2] }, 264 1.0) 265 }; 266 LinearProblem problem = new LinearProblem(m); 267 268 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 269 double initialCost = estimator.getRMS(problem); 270 estimator.estimate(problem); 271 assertTrue(estimator.getRMS(problem) < initialCost); 272 assertTrue(Math.sqrt(m.length) * estimator.getRMS(problem) > 0.6); 273 try { 274 estimator.getCovariances(problem); 275 fail("an exception should have been thrown"); 276 } catch (EstimationException ee) { 277 // expected behavior 278 } catch (Exception e) { 279 fail("wrong exception caught"); 280 } 281 double dJ0 = 2 * (m[0].getResidual() * m[0].getPartial(p[0]) 282 + m[1].getResidual() * m[1].getPartial(p[0]) 283 + m[2].getResidual() * m[2].getPartial(p[0])); 284 double dJ1 = 2 * (m[0].getResidual() * m[0].getPartial(p[1]) 285 + m[1].getResidual() * m[1].getPartial(p[1])); 286 double dJ2 = 2 * (m[0].getResidual() * m[0].getPartial(p[2]) 287 + m[1].getResidual() * m[1].getPartial(p[2]) 288 + m[2].getResidual() * m[2].getPartial(p[2])); 289 assertEquals(0, dJ0, 1.0e-10); 290 assertEquals(0, dJ1, 1.0e-10); 291 assertEquals(0, dJ2, 1.0e-10); 292 293 } 294 295 public void testIllConditioned() throws EstimationException { 296 EstimatedParameter[] p = { 297 new EstimatedParameter("p0", 0), 298 new EstimatedParameter("p1", 1), 299 new EstimatedParameter("p2", 2), 300 new EstimatedParameter("p3", 3) 301 }; 302 303 LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] { 304 new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 }, 305 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 306 32.0), 307 new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 }, 308 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 309 23.0), 310 new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 }, 311 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 312 33.0), 313 new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 }, 314 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 315 31.0) 316 }); 317 LevenbergMarquardtEstimator estimator1 = new LevenbergMarquardtEstimator(); 318 estimator1.estimate(problem1); 319 assertEquals(0, estimator1.getRMS(problem1), 1.0e-10); 320 assertEquals(1.0, p[0].getEstimate(), 1.0e-10); 321 assertEquals(1.0, p[1].getEstimate(), 1.0e-10); 322 assertEquals(1.0, p[2].getEstimate(), 1.0e-10); 323 assertEquals(1.0, p[3].getEstimate(), 1.0e-10); 324 325 LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] { 326 new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 }, 327 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 328 32.0), 329 new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 }, 330 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 331 23.0), 332 new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 }, 333 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 334 33.0), 335 new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 }, 336 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 337 31.0) 338 }); 339 LevenbergMarquardtEstimator estimator2 = new LevenbergMarquardtEstimator(); 340 estimator2.estimate(problem2); 341 assertEquals(0, estimator2.getRMS(problem2), 1.0e-10); 342 assertEquals(-81.0, p[0].getEstimate(), 1.0e-8); 343 assertEquals(137.0, p[1].getEstimate(), 1.0e-8); 344 assertEquals(-34.0, p[2].getEstimate(), 1.0e-8); 345 assertEquals( 22.0, p[3].getEstimate(), 1.0e-8); 346 347 } 348 349 public void testMoreEstimatedParametersSimple() throws EstimationException { 350 351 EstimatedParameter[] p = { 352 new EstimatedParameter("p0", 7), 353 new EstimatedParameter("p1", 6), 354 new EstimatedParameter("p2", 5), 355 new EstimatedParameter("p3", 4) 356 }; 357 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 358 new LinearMeasurement(new double[] { 3.0, 2.0 }, 359 new EstimatedParameter[] { p[0], p[1] }, 360 7.0), 361 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, 362 new EstimatedParameter[] { p[1], p[2], p[3] }, 363 3.0), 364 new LinearMeasurement(new double[] { 2.0, 1.0 }, 365 new EstimatedParameter[] { p[0], p[2] }, 366 5.0) 367 }); 368 369 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 370 estimator.estimate(problem); 371 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 372 373 } 374 375 public void testMoreEstimatedParametersUnsorted() throws EstimationException { 376 EstimatedParameter[] p = { 377 new EstimatedParameter("p0", 2), 378 new EstimatedParameter("p1", 2), 379 new EstimatedParameter("p2", 2), 380 new EstimatedParameter("p3", 2), 381 new EstimatedParameter("p4", 2), 382 new EstimatedParameter("p5", 2) 383 }; 384 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 385 new LinearMeasurement(new double[] { 1.0, 1.0 }, 386 new EstimatedParameter[] { p[0], p[1] }, 387 3.0), 388 new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 }, 389 new EstimatedParameter[] { p[2], p[3], p[4] }, 390 12.0), 391 new LinearMeasurement(new double[] { 1.0, -1.0 }, 392 new EstimatedParameter[] { p[4], p[5] }, 393 -1.0), 394 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, 395 new EstimatedParameter[] { p[3], p[2], p[5] }, 396 7.0), 397 new LinearMeasurement(new double[] { 1.0, -1.0 }, 398 new EstimatedParameter[] { p[4], p[3] }, 399 1.0) 400 }); 401 402 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 403 estimator.estimate(problem); 404 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 405 assertEquals(3.0, p[2].getEstimate(), 1.0e-10); 406 assertEquals(4.0, p[3].getEstimate(), 1.0e-10); 407 assertEquals(5.0, p[4].getEstimate(), 1.0e-10); 408 assertEquals(6.0, p[5].getEstimate(), 1.0e-10); 409 410 } 411 412 public void testRedundantEquations() throws EstimationException { 413 EstimatedParameter[] p = { 414 new EstimatedParameter("p0", 1), 415 new EstimatedParameter("p1", 1) 416 }; 417 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 418 new LinearMeasurement(new double[] { 1.0, 1.0 }, 419 new EstimatedParameter[] { p[0], p[1] }, 420 3.0), 421 new LinearMeasurement(new double[] { 1.0, -1.0 }, 422 new EstimatedParameter[] { p[0], p[1] }, 423 1.0), 424 new LinearMeasurement(new double[] { 1.0, 3.0 }, 425 new EstimatedParameter[] { p[0], p[1] }, 426 5.0) 427 }); 428 429 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 430 estimator.estimate(problem); 431 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 432 assertEquals(2.0, p[0].getEstimate(), 1.0e-10); 433 assertEquals(1.0, p[1].getEstimate(), 1.0e-10); 434 435 } 436 437 public void testInconsistentEquations() throws EstimationException { 438 EstimatedParameter[] p = { 439 new EstimatedParameter("p0", 1), 440 new EstimatedParameter("p1", 1) 441 }; 442 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 443 new LinearMeasurement(new double[] { 1.0, 1.0 }, 444 new EstimatedParameter[] { p[0], p[1] }, 445 3.0), 446 new LinearMeasurement(new double[] { 1.0, -1.0 }, 447 new EstimatedParameter[] { p[0], p[1] }, 448 1.0), 449 new LinearMeasurement(new double[] { 1.0, 3.0 }, 450 new EstimatedParameter[] { p[0], p[1] }, 451 4.0) 452 }); 453 454 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 455 estimator.estimate(problem); 456 assertTrue(estimator.getRMS(problem) > 0.1); 457 458 } 459 460 public void testControlParameters() { 461 Circle circle = new Circle(98.680, 47.345); 462 circle.addPoint( 30.0, 68.0); 463 circle.addPoint( 50.0, -6.0); 464 circle.addPoint(110.0, -20.0); 465 circle.addPoint( 35.0, 15.0); 466 circle.addPoint( 45.0, 97.0); 467 checkEstimate(circle, 0.1, 10, 1.0e-14, 1.0e-16, 1.0e-10, false); 468 checkEstimate(circle, 0.1, 10, 1.0e-15, 1.0e-17, 1.0e-10, true); 469 checkEstimate(circle, 0.1, 5, 1.0e-15, 1.0e-16, 1.0e-10, true); 470 circle.addPoint(300, -300); 471 checkEstimate(circle, 0.1, 20, 1.0e-18, 1.0e-16, 1.0e-10, true); 472 } 473 474 private void checkEstimate(EstimationProblem problem, 475 double initialStepBoundFactor, int maxCostEval, 476 double costRelativeTolerance, double parRelativeTolerance, 477 double orthoTolerance, boolean shouldFail) { 478 try { 479 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 480 estimator.setInitialStepBoundFactor(initialStepBoundFactor); 481 estimator.setMaxCostEval(maxCostEval); 482 estimator.setCostRelativeTolerance(costRelativeTolerance); 483 estimator.setParRelativeTolerance(parRelativeTolerance); 484 estimator.setOrthoTolerance(orthoTolerance); 485 estimator.estimate(problem); 486 assertTrue(! shouldFail); 487 } catch (EstimationException ee) { 488 assertTrue(shouldFail); 489 } catch (Exception e) { 490 fail("wrong exception type caught"); 491 } 492 } 493 494 public void testCircleFitting() throws EstimationException { 495 Circle circle = new Circle(98.680, 47.345); 496 circle.addPoint( 30.0, 68.0); 497 circle.addPoint( 50.0, -6.0); 498 circle.addPoint(110.0, -20.0); 499 circle.addPoint( 35.0, 15.0); 500 circle.addPoint( 45.0, 97.0); 501 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 502 estimator.estimate(circle); 503 assertTrue(estimator.getCostEvaluations() < 10); 504 assertTrue(estimator.getJacobianEvaluations() < 10); 505 double rms = estimator.getRMS(circle); 506 assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10); 507 assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10); 508 assertEquals(96.07590211815305, circle.getX(), 1.0e-10); 509 assertEquals(48.13516790438953, circle.getY(), 1.0e-10); 510 double[][] cov = estimator.getCovariances(circle); 511 assertEquals(1.839, cov[0][0], 0.001); 512 assertEquals(0.731, cov[0][1], 0.001); 513 assertEquals(cov[0][1], cov[1][0], 1.0e-14); 514 assertEquals(0.786, cov[1][1], 0.001); 515 double[] errors = estimator.guessParametersErrors(circle); 516 assertEquals(1.384, errors[0], 0.001); 517 assertEquals(0.905, errors[1], 0.001); 518 519 // add perfect measurements and check errors are reduced 520 double cx = circle.getX(); 521 double cy = circle.getY(); 522 double r = circle.getRadius(); 523 for (double d= 0; d < 2 * Math.PI; d += 0.01) { 524 circle.addPoint(cx + r * Math.cos(d), cy + r * Math.sin(d)); 525 } 526 estimator = new LevenbergMarquardtEstimator(); 527 estimator.estimate(circle); 528 cov = estimator.getCovariances(circle); 529 assertEquals(0.004, cov[0][0], 0.001); 530 assertEquals(6.40e-7, cov[0][1], 1.0e-9); 531 assertEquals(cov[0][1], cov[1][0], 1.0e-14); 532 assertEquals(0.003, cov[1][1], 0.001); 533 errors = estimator.guessParametersErrors(circle); 534 assertEquals(0.004, errors[0], 0.001); 535 assertEquals(0.004, errors[1], 0.001); 536 537 } 538 539 public void testCircleFittingBadInit() throws EstimationException { 540 Circle circle = new Circle(-12, -12); 541 double[][] points = new double[][] { 542 {-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724}, 543 {-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619}, 544 {-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832}, 545 {-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235}, 546 { 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201}, 547 { 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718}, 548 {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862}, 549 {-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526}, 550 {-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398}, 551 {-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513}, 552 {-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737}, 553 { 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850}, 554 { 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138}, 555 {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578}, 556 {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926}, 557 {-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068}, 558 {-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119}, 559 {-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560}, 560 { 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807}, 561 { 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174}, 562 { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635}, 563 {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251}, 564 {-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597}, 565 {-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428}, 566 {-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380}, 567 {-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077}, 568 { 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681}, 569 { 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022}, 570 {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526} 571 }; 572 for (int i = 0; i < points.length; ++i) { 573 circle.addPoint(points[i][0], points[i][1]); 574 } 575 LevenbergMarquardtEstimator estimator = new LevenbergMarquardtEstimator(); 576 estimator.estimate(circle); 577 assertTrue(estimator.getCostEvaluations() < 15); 578 assertTrue(estimator.getJacobianEvaluations() < 10); 579 assertEquals( 0.030184491196225207, estimator.getRMS(circle), 1.0e-9); 580 assertEquals( 0.2922350065939634, circle.getRadius(), 1.0e-9); 581 assertEquals(-0.15173845023862165, circle.getX(), 1.0e-8); 582 assertEquals( 0.20750021499570379, circle.getY(), 1.0e-8); 583 } 584 585 public void testMath199() { 586 try { 587 QuadraticProblem problem = new QuadraticProblem(); 588 problem.addPoint (0, -3.182591015485607, 0.0); 589 problem.addPoint (1, -2.5581184967730577, 4.4E-323); 590 problem.addPoint (2, -2.1488478161387325, 1.0); 591 problem.addPoint (3, -1.9122489313410047, 4.4E-323); 592 problem.addPoint (4, 1.7785661310051026, 0.0); 593 new LevenbergMarquardtEstimator().estimate(problem); 594 fail("an exception should have been thrown"); 595 } catch (EstimationException ee) { 596 // expected behavior 597 } 598 599 } 600 601 private static class LinearProblem implements EstimationProblem { 602 603 public LinearProblem(LinearMeasurement[] measurements) { 604 this.measurements = measurements; 605 } 606 607 public WeightedMeasurement[] getMeasurements() { 608 return measurements; 609 } 610 611 public EstimatedParameter[] getUnboundParameters() { 612 return getAllParameters(); 613 } 614 615 public EstimatedParameter[] getAllParameters() { 616 HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>(); 617 for (int i = 0; i < measurements.length; ++i) { 618 EstimatedParameter[] parameters = measurements[i].getParameters(); 619 for (int j = 0; j < parameters.length; ++j) { 620 set.add(parameters[j]); 621 } 622 } 623 return set.toArray(new EstimatedParameter[set.size()]); 624 } 625 626 private LinearMeasurement[] measurements; 627 628 } 629 630 private static class LinearMeasurement extends WeightedMeasurement { 631 632 public LinearMeasurement(double[] factors, EstimatedParameter[] parameters, 633 double setPoint) { 634 super(1.0, setPoint); 635 this.factors = factors; 636 this.parameters = parameters; 637 } 638 639 @Override 640 public double getTheoreticalValue() { 641 double v = 0; 642 for (int i = 0; i < factors.length; ++i) { 643 v += factors[i] * parameters[i].getEstimate(); 644 } 645 return v; 646 } 647 648 @Override 649 public double getPartial(EstimatedParameter parameter) { 650 for (int i = 0; i < parameters.length; ++i) { 651 if (parameters[i] == parameter) { 652 return factors[i]; 653 } 654 } 655 return 0; 656 } 657 658 public EstimatedParameter[] getParameters() { 659 return parameters; 660 } 661 662 private double[] factors; 663 private EstimatedParameter[] parameters; 664 private static final long serialVersionUID = -3922448707008868580L; 665 666 } 667 668 private static class Circle implements EstimationProblem { 669 670 public Circle(double cx, double cy) { 671 this.cx = new EstimatedParameter("cx", cx); 672 this.cy = new EstimatedParameter("cy", cy); 673 points = new ArrayList<PointModel>(); 674 } 675 676 public void addPoint(double px, double py) { 677 points.add(new PointModel(this, px, py)); 678 } 679 680 public int getM() { 681 return points.size(); 682 } 683 684 public WeightedMeasurement[] getMeasurements() { 685 return points.toArray(new PointModel[points.size()]); 686 } 687 688 public EstimatedParameter[] getAllParameters() { 689 return new EstimatedParameter[] { cx, cy }; 690 } 691 692 public EstimatedParameter[] getUnboundParameters() { 693 return new EstimatedParameter[] { cx, cy }; 694 } 695 696 public double getPartialRadiusX() { 697 double dRdX = 0; 698 for (PointModel point : points) { 699 dRdX += point.getPartialDiX(); 700 } 701 return dRdX / points.size(); 702 } 703 704 public double getPartialRadiusY() { 705 double dRdY = 0; 706 for (PointModel point : points) { 707 dRdY += point.getPartialDiY(); 708 } 709 return dRdY / points.size(); 710 } 711 712 public double getRadius() { 713 double r = 0; 714 for (PointModel point : points) { 715 r += point.getCenterDistance(); 716 } 717 return r / points.size(); 718 } 719 720 public double getX() { 721 return cx.getEstimate(); 722 } 723 724 public double getY() { 725 return cy.getEstimate(); 726 } 727 728 private static class PointModel extends WeightedMeasurement { 729 730 public PointModel(Circle circle, double px, double py) { 731 super(1.0, 0.0); 732 this.px = px; 733 this.py = py; 734 this.circle = circle; 735 } 736 737 @Override 738 public double getPartial(EstimatedParameter parameter) { 739 if (parameter == circle.cx) { 740 return getPartialDiX() - circle.getPartialRadiusX(); 741 } else if (parameter == circle.cy) { 742 return getPartialDiY() - circle.getPartialRadiusY(); 743 } 744 return 0; 745 } 746 747 public double getCenterDistance() { 748 double dx = px - circle.cx.getEstimate(); 749 double dy = py - circle.cy.getEstimate(); 750 return Math.sqrt(dx * dx + dy * dy); 751 } 752 753 public double getPartialDiX() { 754 return (circle.cx.getEstimate() - px) / getCenterDistance(); 755 } 756 757 public double getPartialDiY() { 758 return (circle.cy.getEstimate() - py) / getCenterDistance(); 759 } 760 761 @Override 762 public double getTheoreticalValue() { 763 return getCenterDistance() - circle.getRadius(); 764 } 765 766 private double px; 767 private double py; 768 private transient final Circle circle; 769 private static final long serialVersionUID = 1L; 770 771 } 772 773 private EstimatedParameter cx; 774 private EstimatedParameter cy; 775 private ArrayList<PointModel> points; 776 777 } 778 779 private static class QuadraticProblem extends SimpleEstimationProblem { 780 781 private EstimatedParameter a; 782 private EstimatedParameter b; 783 private EstimatedParameter c; 784 785 public QuadraticProblem() { 786 a = new EstimatedParameter("a", 0.0); 787 b = new EstimatedParameter("b", 0.0); 788 c = new EstimatedParameter("c", 0.0); 789 addParameter(a); 790 addParameter(b); 791 addParameter(c); 792 } 793 794 public void addPoint(double x, double y, double w) { 795 addMeasurement(new LocalMeasurement(this, x, y, w)); 796 } 797 798 public double theoreticalValue(double x) { 799 return ( (a.getEstimate() * x + b.getEstimate() ) * x + c.getEstimate()); 800 } 801 802 private double partial(double x, EstimatedParameter parameter) { 803 if (parameter == a) { 804 return x * x; 805 } else if (parameter == b) { 806 return x; 807 } else { 808 return 1.0; 809 } 810 } 811 812 private static class LocalMeasurement extends WeightedMeasurement { 813 814 private static final long serialVersionUID = 1555043155023729130L; 815 private final double x; 816 private transient final QuadraticProblem pb; 817 818 // constructor 819 public LocalMeasurement(QuadraticProblem pb, double x, double y, double w) { 820 super(w, y); 821 this.x = x; 822 this.pb = pb; 823 } 824 825 @Override 826 public double getTheoreticalValue() { 827 return pb.theoreticalValue(x); 828 } 829 830 @Override 831 public double getPartial(EstimatedParameter parameter) { 832 return pb.partial(x, parameter); 833 } 834 835 } 836 } 837 838 public static Test suite() { 839 return new TestSuite(LevenbergMarquardtEstimatorTest.class); 840 } 841 842 }