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 GaussNewtonEstimatorTest 091 extends TestCase { 092 093 public GaussNewtonEstimatorTest(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 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 106 estimator.estimate(problem); 107 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 108 assertEquals(1.5, 109 problem.getUnboundParameters()[0].getEstimate(), 110 1.0e-10); 111 } 112 113 public void testQRColumnsPermutation() throws EstimationException { 114 115 EstimatedParameter[] x = { 116 new EstimatedParameter("p0", 0), new EstimatedParameter("p1", 0) 117 }; 118 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 119 new LinearMeasurement(new double[] { 1.0, -1.0 }, 120 new EstimatedParameter[] { x[0], x[1] }, 121 4.0), 122 new LinearMeasurement(new double[] { 2.0 }, 123 new EstimatedParameter[] { x[1] }, 124 6.0), 125 new LinearMeasurement(new double[] { 1.0, -2.0 }, 126 new EstimatedParameter[] { x[0], x[1] }, 127 1.0) 128 }); 129 130 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 131 estimator.estimate(problem); 132 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 133 assertEquals(7.0, x[0].getEstimate(), 1.0e-10); 134 assertEquals(3.0, x[1].getEstimate(), 1.0e-10); 135 136 } 137 138 public void testNoDependency() throws EstimationException { 139 EstimatedParameter[] p = new EstimatedParameter[] { 140 new EstimatedParameter("p0", 0), 141 new EstimatedParameter("p1", 0), 142 new EstimatedParameter("p2", 0), 143 new EstimatedParameter("p3", 0), 144 new EstimatedParameter("p4", 0), 145 new EstimatedParameter("p5", 0) 146 }; 147 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 148 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[0] }, 0.0), 149 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[1] }, 1.1), 150 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[2] }, 2.2), 151 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[3] }, 3.3), 152 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[4] }, 4.4), 153 new LinearMeasurement(new double[] {2}, new EstimatedParameter[] { p[5] }, 5.5) 154 }); 155 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 156 estimator.estimate(problem); 157 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 158 for (int i = 0; i < p.length; ++i) { 159 assertEquals(0.55 * i, p[i].getEstimate(), 1.0e-10); 160 } 161 } 162 163 public void testOneSet() throws EstimationException { 164 165 EstimatedParameter[] p = { 166 new EstimatedParameter("p0", 0), 167 new EstimatedParameter("p1", 0), 168 new EstimatedParameter("p2", 0) 169 }; 170 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 171 new LinearMeasurement(new double[] { 1.0 }, 172 new EstimatedParameter[] { p[0] }, 173 1.0), 174 new LinearMeasurement(new double[] { -1.0, 1.0 }, 175 new EstimatedParameter[] { p[0], p[1] }, 176 1.0), 177 new LinearMeasurement(new double[] { -1.0, 1.0 }, 178 new EstimatedParameter[] { p[1], p[2] }, 179 1.0) 180 }); 181 182 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 183 estimator.estimate(problem); 184 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 185 assertEquals(1.0, p[0].getEstimate(), 1.0e-10); 186 assertEquals(2.0, p[1].getEstimate(), 1.0e-10); 187 assertEquals(3.0, p[2].getEstimate(), 1.0e-10); 188 189 } 190 191 public void testTwoSets() throws EstimationException { 192 EstimatedParameter[] p = { 193 new EstimatedParameter("p0", 0), 194 new EstimatedParameter("p1", 1), 195 new EstimatedParameter("p2", 2), 196 new EstimatedParameter("p3", 3), 197 new EstimatedParameter("p4", 4), 198 new EstimatedParameter("p5", 5) 199 }; 200 201 double epsilon = 1.0e-7; 202 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 203 204 // 4 elements sub-problem 205 new LinearMeasurement(new double[] { 2.0, 1.0, 4.0 }, 206 new EstimatedParameter[] { p[0], p[1], p[3] }, 207 2.0), 208 new LinearMeasurement(new double[] { -4.0, -2.0, 3.0, -7.0 }, 209 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 210 -9.0), 211 new LinearMeasurement(new double[] { 4.0, 1.0, -2.0, 8.0 }, 212 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 213 2.0), 214 new LinearMeasurement(new double[] { -3.0, -12.0, -1.0 }, 215 new EstimatedParameter[] { p[1], p[2], p[3] }, 216 2.0), 217 218 // 2 elements sub-problem 219 new LinearMeasurement(new double[] { epsilon, 1.0 }, 220 new EstimatedParameter[] { p[4], p[5] }, 221 1.0 + epsilon * epsilon), 222 new LinearMeasurement(new double[] { 1.0, 1.0 }, 223 new EstimatedParameter[] { p[4], p[5] }, 224 2.0) 225 226 }); 227 228 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 229 estimator.estimate(problem); 230 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 231 assertEquals( 3.0, p[0].getEstimate(), 1.0e-10); 232 assertEquals( 4.0, p[1].getEstimate(), 1.0e-10); 233 assertEquals(-1.0, p[2].getEstimate(), 1.0e-10); 234 assertEquals(-2.0, p[3].getEstimate(), 1.0e-10); 235 assertEquals( 1.0 + epsilon, p[4].getEstimate(), 1.0e-10); 236 assertEquals( 1.0 - epsilon, p[5].getEstimate(), 1.0e-10); 237 238 } 239 240 public void testNonInversible() { 241 242 EstimatedParameter[] p = { 243 new EstimatedParameter("p0", 0), 244 new EstimatedParameter("p1", 0), 245 new EstimatedParameter("p2", 0) 246 }; 247 LinearMeasurement[] m = new LinearMeasurement[] { 248 new LinearMeasurement(new double[] { 1.0, 2.0, -3.0 }, 249 new EstimatedParameter[] { p[0], p[1], p[2] }, 250 1.0), 251 new LinearMeasurement(new double[] { 2.0, 1.0, 3.0 }, 252 new EstimatedParameter[] { p[0], p[1], p[2] }, 253 1.0), 254 new LinearMeasurement(new double[] { -3.0, -9.0 }, 255 new EstimatedParameter[] { p[0], p[2] }, 256 1.0) 257 }; 258 LinearProblem problem = new LinearProblem(m); 259 260 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 261 try { 262 estimator.estimate(problem); 263 fail("an exception should have been caught"); 264 } catch (EstimationException ee) { 265 // expected behavior 266 } catch (Exception e) { 267 fail("wrong exception type caught"); 268 } 269 } 270 271 public void testIllConditioned() throws EstimationException { 272 EstimatedParameter[] p = { 273 new EstimatedParameter("p0", 0), 274 new EstimatedParameter("p1", 1), 275 new EstimatedParameter("p2", 2), 276 new EstimatedParameter("p3", 3) 277 }; 278 279 LinearProblem problem1 = new LinearProblem(new LinearMeasurement[] { 280 new LinearMeasurement(new double[] { 10.0, 7.0, 8.0, 7.0 }, 281 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 282 32.0), 283 new LinearMeasurement(new double[] { 7.0, 5.0, 6.0, 5.0 }, 284 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 285 23.0), 286 new LinearMeasurement(new double[] { 8.0, 6.0, 10.0, 9.0 }, 287 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 288 33.0), 289 new LinearMeasurement(new double[] { 7.0, 5.0, 9.0, 10.0 }, 290 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 291 31.0) 292 }); 293 GaussNewtonEstimator estimator1 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 294 estimator1.estimate(problem1); 295 assertEquals(0, estimator1.getRMS(problem1), 1.0e-10); 296 assertEquals(1.0, p[0].getEstimate(), 1.0e-10); 297 assertEquals(1.0, p[1].getEstimate(), 1.0e-10); 298 assertEquals(1.0, p[2].getEstimate(), 1.0e-10); 299 assertEquals(1.0, p[3].getEstimate(), 1.0e-10); 300 301 LinearProblem problem2 = new LinearProblem(new LinearMeasurement[] { 302 new LinearMeasurement(new double[] { 10.0, 7.0, 8.1, 7.2 }, 303 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 304 32.0), 305 new LinearMeasurement(new double[] { 7.08, 5.04, 6.0, 5.0 }, 306 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 307 23.0), 308 new LinearMeasurement(new double[] { 8.0, 5.98, 9.89, 9.0 }, 309 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 310 33.0), 311 new LinearMeasurement(new double[] { 6.99, 4.99, 9.0, 9.98 }, 312 new EstimatedParameter[] { p[0], p[1], p[2], p[3] }, 313 31.0) 314 }); 315 GaussNewtonEstimator estimator2 = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 316 estimator2.estimate(problem2); 317 assertEquals(0, estimator2.getRMS(problem2), 1.0e-10); 318 assertEquals(-81.0, p[0].getEstimate(), 1.0e-8); 319 assertEquals(137.0, p[1].getEstimate(), 1.0e-8); 320 assertEquals(-34.0, p[2].getEstimate(), 1.0e-8); 321 assertEquals( 22.0, p[3].getEstimate(), 1.0e-8); 322 323 } 324 325 public void testMoreEstimatedParametersSimple() { 326 327 EstimatedParameter[] p = { 328 new EstimatedParameter("p0", 7), 329 new EstimatedParameter("p1", 6), 330 new EstimatedParameter("p2", 5), 331 new EstimatedParameter("p3", 4) 332 }; 333 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 334 new LinearMeasurement(new double[] { 3.0, 2.0 }, 335 new EstimatedParameter[] { p[0], p[1] }, 336 7.0), 337 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, 338 new EstimatedParameter[] { p[1], p[2], p[3] }, 339 3.0), 340 new LinearMeasurement(new double[] { 2.0, 1.0 }, 341 new EstimatedParameter[] { p[0], p[2] }, 342 5.0) 343 }); 344 345 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 346 try { 347 estimator.estimate(problem); 348 fail("an exception should have been caught"); 349 } catch (EstimationException ee) { 350 // expected behavior 351 } catch (Exception e) { 352 fail("wrong exception type caught"); 353 } 354 355 } 356 357 public void testMoreEstimatedParametersUnsorted() { 358 EstimatedParameter[] p = { 359 new EstimatedParameter("p0", 2), 360 new EstimatedParameter("p1", 2), 361 new EstimatedParameter("p2", 2), 362 new EstimatedParameter("p3", 2), 363 new EstimatedParameter("p4", 2), 364 new EstimatedParameter("p5", 2) 365 }; 366 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 367 new LinearMeasurement(new double[] { 1.0, 1.0 }, 368 new EstimatedParameter[] { p[0], p[1] }, 369 3.0), 370 new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 }, 371 new EstimatedParameter[] { p[2], p[3], p[4] }, 372 12.0), 373 new LinearMeasurement(new double[] { 1.0, -1.0 }, 374 new EstimatedParameter[] { p[4], p[5] }, 375 -1.0), 376 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, 377 new EstimatedParameter[] { p[3], p[2], p[5] }, 378 7.0), 379 new LinearMeasurement(new double[] { 1.0, -1.0 }, 380 new EstimatedParameter[] { p[4], p[3] }, 381 1.0) 382 }); 383 384 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 385 try { 386 estimator.estimate(problem); 387 fail("an exception should have been caught"); 388 } catch (EstimationException ee) { 389 // expected behavior 390 } catch (Exception e) { 391 fail("wrong exception type caught"); 392 } 393 394 } 395 396 public void testRedundantEquations() throws EstimationException { 397 EstimatedParameter[] p = { 398 new EstimatedParameter("p0", 1), 399 new EstimatedParameter("p1", 1) 400 }; 401 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 402 new LinearMeasurement(new double[] { 1.0, 1.0 }, 403 new EstimatedParameter[] { p[0], p[1] }, 404 3.0), 405 new LinearMeasurement(new double[] { 1.0, -1.0 }, 406 new EstimatedParameter[] { p[0], p[1] }, 407 1.0), 408 new LinearMeasurement(new double[] { 1.0, 3.0 }, 409 new EstimatedParameter[] { p[0], p[1] }, 410 5.0) 411 }); 412 413 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 414 estimator.estimate(problem); 415 assertEquals(0, estimator.getRMS(problem), 1.0e-10); 416 EstimatedParameter[] all = problem.getAllParameters(); 417 for (int i = 0; i < all.length; ++i) { 418 assertEquals(all[i].getName().equals("p0") ? 2.0 : 1.0, 419 all[i].getEstimate(), 1.0e-10); 420 } 421 422 } 423 424 public void testInconsistentEquations() throws EstimationException { 425 EstimatedParameter[] p = { 426 new EstimatedParameter("p0", 1), 427 new EstimatedParameter("p1", 1) 428 }; 429 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 430 new LinearMeasurement(new double[] { 1.0, 1.0 }, 431 new EstimatedParameter[] { p[0], p[1] }, 432 3.0), 433 new LinearMeasurement(new double[] { 1.0, -1.0 }, 434 new EstimatedParameter[] { p[0], p[1] }, 435 1.0), 436 new LinearMeasurement(new double[] { 1.0, 3.0 }, 437 new EstimatedParameter[] { p[0], p[1] }, 438 4.0) 439 }); 440 441 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 442 estimator.estimate(problem); 443 assertTrue(estimator.getRMS(problem) > 0.1); 444 445 } 446 447 public void testBoundParameters() throws EstimationException { 448 EstimatedParameter[] p = { 449 new EstimatedParameter("unbound0", 2, false), 450 new EstimatedParameter("unbound1", 2, false), 451 new EstimatedParameter("bound", 2, true) 452 }; 453 LinearProblem problem = new LinearProblem(new LinearMeasurement[] { 454 new LinearMeasurement(new double[] { 1.0, 1.0, 1.0 }, 455 new EstimatedParameter[] { p[0], p[1], p[2] }, 456 3.0), 457 new LinearMeasurement(new double[] { 1.0, -1.0, 1.0 }, 458 new EstimatedParameter[] { p[0], p[1], p[2] }, 459 1.0), 460 new LinearMeasurement(new double[] { 1.0, 3.0, 2.0 }, 461 new EstimatedParameter[] { p[0], p[1], p[2] }, 462 7.0) 463 }); 464 465 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 466 estimator.estimate(problem); 467 assertTrue(estimator.getRMS(problem) < 1.0e-10); 468 double[][] covariances = estimator.getCovariances(problem); 469 int i0 = 0, i1 = 1; 470 if (problem.getUnboundParameters()[0].getName().endsWith("1")) { 471 i0 = 1; 472 i1 = 0; 473 } 474 assertEquals(11.0 / 24, covariances[i0][i0], 1.0e-10); 475 assertEquals(-3.0 / 24, covariances[i0][i1], 1.0e-10); 476 assertEquals(-3.0 / 24, covariances[i1][i0], 1.0e-10); 477 assertEquals( 3.0 / 24, covariances[i1][i1], 1.0e-10); 478 479 double[] errors = estimator.guessParametersErrors(problem); 480 assertEquals(0, errors[i0], 1.0e-10); 481 assertEquals(0, errors[i1], 1.0e-10); 482 483 } 484 485 public void testMaxIterations() { 486 Circle circle = new Circle(98.680, 47.345); 487 circle.addPoint( 30.0, 68.0); 488 circle.addPoint( 50.0, -6.0); 489 circle.addPoint(110.0, -20.0); 490 circle.addPoint( 35.0, 15.0); 491 circle.addPoint( 45.0, 97.0); 492 try { 493 GaussNewtonEstimator estimator = new GaussNewtonEstimator(4, 1.0e-14, 1.0e-14); 494 estimator.estimate(circle); 495 fail("an exception should have been caught"); 496 } catch (EstimationException ee) { 497 // expected behavior 498 } catch (Exception e) { 499 fail("wrong exception type caught"); 500 } 501 } 502 503 public void testCircleFitting() throws EstimationException { 504 Circle circle = new Circle(98.680, 47.345); 505 circle.addPoint( 30.0, 68.0); 506 circle.addPoint( 50.0, -6.0); 507 circle.addPoint(110.0, -20.0); 508 circle.addPoint( 35.0, 15.0); 509 circle.addPoint( 45.0, 97.0); 510 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-10, 1.0e-10); 511 estimator.estimate(circle); 512 double rms = estimator.getRMS(circle); 513 assertEquals(1.768262623567235, Math.sqrt(circle.getM()) * rms, 1.0e-10); 514 assertEquals(69.96016176931406, circle.getRadius(), 1.0e-10); 515 assertEquals(96.07590211815305, circle.getX(), 1.0e-10); 516 assertEquals(48.13516790438953, circle.getY(), 1.0e-10); 517 } 518 519 public void testCircleFittingBadInit() { 520 Circle circle = new Circle(-12, -12); 521 double[][] points = new double[][] { 522 {-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724}, 523 {-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619}, 524 {-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832}, 525 {-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235}, 526 { 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201}, 527 { 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718}, 528 {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862}, 529 {-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526}, 530 {-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398}, 531 {-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513}, 532 {-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737}, 533 { 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850}, 534 { 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138}, 535 {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578}, 536 {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926}, 537 {-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068}, 538 {-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119}, 539 {-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560}, 540 { 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807}, 541 { 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174}, 542 { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635}, 543 {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251}, 544 {-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597}, 545 {-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428}, 546 {-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380}, 547 {-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077}, 548 { 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681}, 549 { 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022}, 550 {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526} 551 }; 552 for (int i = 0; i < points.length; ++i) { 553 circle.addPoint(points[i][0], points[i][1]); 554 } 555 GaussNewtonEstimator estimator = new GaussNewtonEstimator(100, 1.0e-6, 1.0e-6); 556 try { 557 estimator.estimate(circle); 558 fail("an exception should have been caught"); 559 } catch (EstimationException ee) { 560 // expected behavior 561 } catch (Exception e) { 562 fail("wrong exception type caught"); 563 } 564 } 565 566 private static class LinearProblem extends SimpleEstimationProblem { 567 568 public LinearProblem(LinearMeasurement[] measurements) { 569 HashSet<EstimatedParameter> set = new HashSet<EstimatedParameter>(); 570 for (int i = 0; i < measurements.length; ++i) { 571 addMeasurement(measurements[i]); 572 EstimatedParameter[] parameters = measurements[i].getParameters(); 573 for (int j = 0; j < parameters.length; ++j) { 574 set.add(parameters[j]); 575 } 576 } 577 for (EstimatedParameter p : set) { 578 addParameter(p); 579 } 580 } 581 582 } 583 584 private static class LinearMeasurement extends WeightedMeasurement { 585 586 public LinearMeasurement(double[] factors, EstimatedParameter[] parameters, 587 double setPoint) { 588 super(1.0, setPoint, true); 589 this.factors = factors; 590 this.parameters = parameters; 591 setIgnored(false); 592 } 593 594 @Override 595 public double getTheoreticalValue() { 596 double v = 0; 597 for (int i = 0; i < factors.length; ++i) { 598 v += factors[i] * parameters[i].getEstimate(); 599 } 600 return v; 601 } 602 603 @Override 604 public double getPartial(EstimatedParameter parameter) { 605 for (int i = 0; i < parameters.length; ++i) { 606 if (parameters[i] == parameter) { 607 return factors[i]; 608 } 609 } 610 return 0; 611 } 612 613 public EstimatedParameter[] getParameters() { 614 return parameters; 615 } 616 617 private double[] factors; 618 private EstimatedParameter[] parameters; 619 private static final long serialVersionUID = -3922448707008868580L; 620 621 } 622 623 private static class Circle implements EstimationProblem { 624 625 public Circle(double cx, double cy) { 626 this.cx = new EstimatedParameter("cx", cx); 627 this.cy = new EstimatedParameter(new EstimatedParameter("cy", cy)); 628 points = new ArrayList<PointModel>(); 629 } 630 631 public void addPoint(double px, double py) { 632 points.add(new PointModel(this, px, py)); 633 } 634 635 public int getM() { 636 return points.size(); 637 } 638 639 public WeightedMeasurement[] getMeasurements() { 640 return points.toArray(new PointModel[points.size()]); 641 } 642 643 public EstimatedParameter[] getAllParameters() { 644 return new EstimatedParameter[] { cx, cy }; 645 } 646 647 public EstimatedParameter[] getUnboundParameters() { 648 return new EstimatedParameter[] { cx, cy }; 649 } 650 651 public double getPartialRadiusX() { 652 double dRdX = 0; 653 for (PointModel point : points) { 654 dRdX += point.getPartialDiX(); 655 } 656 return dRdX / points.size(); 657 } 658 659 public double getPartialRadiusY() { 660 double dRdY = 0; 661 for (PointModel point : points) { 662 dRdY += point.getPartialDiY(); 663 } 664 return dRdY / points.size(); 665 } 666 667 public double getRadius() { 668 double r = 0; 669 for (PointModel point : points) { 670 r += point.getCenterDistance(); 671 } 672 return r / points.size(); 673 } 674 675 public double getX() { 676 return cx.getEstimate(); 677 } 678 679 public double getY() { 680 return cy.getEstimate(); 681 } 682 683 private static class PointModel extends WeightedMeasurement { 684 685 public PointModel(Circle circle, double px, double py) { 686 super(1.0, 0.0); 687 this.px = px; 688 this.py = py; 689 this.circle = circle; 690 } 691 692 @Override 693 public double getPartial(EstimatedParameter parameter) { 694 if (parameter == circle.cx) { 695 return getPartialDiX() - circle.getPartialRadiusX(); 696 } else if (parameter == circle.cy) { 697 return getPartialDiY() - circle.getPartialRadiusY(); 698 } 699 return 0; 700 } 701 702 public double getCenterDistance() { 703 double dx = px - circle.cx.getEstimate(); 704 double dy = py - circle.cy.getEstimate(); 705 return Math.sqrt(dx * dx + dy * dy); 706 } 707 708 public double getPartialDiX() { 709 return (circle.cx.getEstimate() - px) / getCenterDistance(); 710 } 711 712 public double getPartialDiY() { 713 return (circle.cy.getEstimate() - py) / getCenterDistance(); 714 } 715 716 @Override 717 public double getTheoreticalValue() { 718 return getCenterDistance() - circle.getRadius(); 719 } 720 721 private double px; 722 private double py; 723 private transient final Circle circle; 724 private static final long serialVersionUID = 1L; 725 726 } 727 728 private EstimatedParameter cx; 729 private EstimatedParameter cy; 730 private ArrayList<PointModel> points; 731 732 } 733 734 public static Test suite() { 735 return new TestSuite(GaussNewtonEstimatorTest.class); 736 } 737 738 }