1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math3.optim.nonlinear.vector.jacobian; 18 19 import java.io.IOException; 20 import java.util.Arrays; 21 22 import org.apache.commons.math3.analysis.MultivariateVectorFunction; 23 import org.apache.commons.math3.analysis.MultivariateMatrixFunction; 24 import org.apache.commons.math3.exception.ConvergenceException; 25 import org.apache.commons.math3.exception.DimensionMismatchException; 26 import org.apache.commons.math3.geometry.euclidean.twod.Vector2D; 27 import org.apache.commons.math3.linear.BlockRealMatrix; 28 import org.apache.commons.math3.linear.RealMatrix; 29 import org.apache.commons.math3.optim.PointVectorValuePair; 30 import org.apache.commons.math3.optim.InitialGuess; 31 import org.apache.commons.math3.optim.MaxEval; 32 import org.apache.commons.math3.optim.nonlinear.vector.Target; 33 import org.apache.commons.math3.optim.nonlinear.vector.Weight; 34 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunction; 35 import org.apache.commons.math3.optim.nonlinear.vector.ModelFunctionJacobian; 36 import org.apache.commons.math3.util.FastMath; 37 import org.junit.Assert; 38 import org.junit.Test; 39 40 /** 41 * <p>Some of the unit tests are re-implementations of the MINPACK <a 42 * href="http://www.netlib.org/minpack/ex/file17">file17</a> and <a 43 * href="http://www.netlib.org/minpack/ex/file22">file22</a> test files. 44 * The redistribution policy for MINPACK is available <a 45 * href="http://www.netlib.org/minpack/disclaimer">here</a>, for 46 * convenience, it is reproduced below.</p> 47 48 * <table border="0" width="80%" cellpadding="10" align="center" bgcolor="#E0E0E0"> 49 * <tr><td> 50 * Minpack Copyright Notice (1999) University of Chicago. 51 * All rights reserved 52 * </td></tr> 53 * <tr><td> 54 * Redistribution and use in source and binary forms, with or without 55 * modification, are permitted provided that the following conditions 56 * are met: 57 * <ol> 58 * <li>Redistributions of source code must retain the above copyright 59 * notice, this list of conditions and the following disclaimer.</li> 60 * <li>Redistributions in binary form must reproduce the above 61 * copyright notice, this list of conditions and the following 62 * disclaimer in the documentation and/or other materials provided 63 * with the distribution.</li> 64 * <li>The end-user documentation included with the redistribution, if any, 65 * must include the following acknowledgment: 66 * <code>This product includes software developed by the University of 67 * Chicago, as Operator of Argonne National Laboratory.</code> 68 * Alternately, this acknowledgment may appear in the software itself, 69 * if and wherever such third-party acknowledgments normally appear.</li> 70 * <li><strong>WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS" 71 * WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE 72 * UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND 73 * THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR 74 * IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES 75 * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE 76 * OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY 77 * OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR 78 * USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF 79 * THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4) 80 * DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION 81 * UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL 82 * BE CORRECTED.</strong></li> 83 * <li><strong>LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT 84 * HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF 85 * ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT, 86 * INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF 87 * ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF 88 * PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER 89 * SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT 90 * (INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE, 91 * EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE 92 * POSSIBILITY OF SUCH LOSS OR DAMAGES.</strong></li> 93 * <ol></td></tr> 94 * </table> 95 96 * @author Argonne National Laboratory. MINPACK project. March 1980 (original fortran minpack tests) 97 * @author Burton S. Garbow (original fortran minpack tests) 98 * @author Kenneth E. Hillstrom (original fortran minpack tests) 99 * @author Jorge J. More (original fortran minpack tests) 100 * @author Luc Maisonobe (non-minpack tests and minpack tests Java translation) 101 */ 102 @Deprecated 103 public abstract class AbstractLeastSquaresOptimizerAbstractTest { 104 createOptimizer()105 public abstract AbstractLeastSquaresOptimizer createOptimizer(); 106 107 @Test testGetIterations()108 public void testGetIterations() { 109 AbstractLeastSquaresOptimizer optim = createOptimizer(); 110 optim.optimize(new MaxEval(100), new Target(new double[] { 1 }), 111 new Weight(new double[] { 1 }), 112 new InitialGuess(new double[] { 3 }), 113 new ModelFunction(new MultivariateVectorFunction() { 114 public double[] value(double[] point) { 115 return new double[] { 116 FastMath.pow(point[0], 4) 117 }; 118 } 119 }), 120 new ModelFunctionJacobian(new MultivariateMatrixFunction() { 121 public double[][] value(double[] point) { 122 return new double[][] { 123 { 0.25 * FastMath.pow(point[0], 3) } 124 }; 125 } 126 })); 127 128 Assert.assertTrue(optim.getIterations() > 0); 129 } 130 131 @Test testTrivial()132 public void testTrivial() { 133 LinearProblem problem 134 = new LinearProblem(new double[][] { { 2 } }, new double[] { 3 }); 135 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 136 PointVectorValuePair optimum = 137 optimizer.optimize(new MaxEval(100), 138 problem.getModelFunction(), 139 problem.getModelFunctionJacobian(), 140 problem.getTarget(), 141 new Weight(new double[] { 1 }), 142 new InitialGuess(new double[] { 0 })); 143 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 144 Assert.assertEquals(1.5, optimum.getPoint()[0], 1e-10); 145 Assert.assertEquals(3.0, optimum.getValue()[0], 1e-10); 146 } 147 148 @Test testQRColumnsPermutation()149 public void testQRColumnsPermutation() { 150 151 LinearProblem problem 152 = new LinearProblem(new double[][] { { 1, -1 }, { 0, 2 }, { 1, -2 } }, 153 new double[] { 4, 6, 1 }); 154 155 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 156 PointVectorValuePair optimum = 157 optimizer.optimize(new MaxEval(100), 158 problem.getModelFunction(), 159 problem.getModelFunctionJacobian(), 160 problem.getTarget(), 161 new Weight(new double[] { 1, 1, 1 }), 162 new InitialGuess(new double[] { 0, 0 })); 163 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 164 Assert.assertEquals(7, optimum.getPoint()[0], 1e-10); 165 Assert.assertEquals(3, optimum.getPoint()[1], 1e-10); 166 Assert.assertEquals(4, optimum.getValue()[0], 1e-10); 167 Assert.assertEquals(6, optimum.getValue()[1], 1e-10); 168 Assert.assertEquals(1, optimum.getValue()[2], 1e-10); 169 } 170 171 @Test testNoDependency()172 public void testNoDependency() { 173 LinearProblem problem = new LinearProblem(new double[][] { 174 { 2, 0, 0, 0, 0, 0 }, 175 { 0, 2, 0, 0, 0, 0 }, 176 { 0, 0, 2, 0, 0, 0 }, 177 { 0, 0, 0, 2, 0, 0 }, 178 { 0, 0, 0, 0, 2, 0 }, 179 { 0, 0, 0, 0, 0, 2 } 180 }, new double[] { 0, 1.1, 2.2, 3.3, 4.4, 5.5 }); 181 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 182 PointVectorValuePair optimum = 183 optimizer.optimize(new MaxEval(100), 184 problem.getModelFunction(), 185 problem.getModelFunctionJacobian(), 186 problem.getTarget(), 187 new Weight(new double[] { 1, 1, 1, 1, 1, 1 }), 188 new InitialGuess(new double[] { 0, 0, 0, 0, 0, 0 })); 189 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 190 for (int i = 0; i < problem.target.length; ++i) { 191 Assert.assertEquals(0.55 * i, optimum.getPoint()[i], 1e-10); 192 } 193 } 194 195 @Test testOneSet()196 public void testOneSet() { 197 198 LinearProblem problem = new LinearProblem(new double[][] { 199 { 1, 0, 0 }, 200 { -1, 1, 0 }, 201 { 0, -1, 1 } 202 }, new double[] { 1, 1, 1}); 203 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 204 PointVectorValuePair optimum = 205 optimizer.optimize(new MaxEval(100), 206 problem.getModelFunction(), 207 problem.getModelFunctionJacobian(), 208 problem.getTarget(), 209 new Weight(new double[] { 1, 1, 1 }), 210 new InitialGuess(new double[] { 0, 0, 0 })); 211 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 212 Assert.assertEquals(1, optimum.getPoint()[0], 1e-10); 213 Assert.assertEquals(2, optimum.getPoint()[1], 1e-10); 214 Assert.assertEquals(3, optimum.getPoint()[2], 1e-10); 215 } 216 217 @Test testTwoSets()218 public void testTwoSets() { 219 double epsilon = 1e-7; 220 LinearProblem problem = new LinearProblem(new double[][] { 221 { 2, 1, 0, 4, 0, 0 }, 222 { -4, -2, 3, -7, 0, 0 }, 223 { 4, 1, -2, 8, 0, 0 }, 224 { 0, -3, -12, -1, 0, 0 }, 225 { 0, 0, 0, 0, epsilon, 1 }, 226 { 0, 0, 0, 0, 1, 1 } 227 }, new double[] { 2, -9, 2, 2, 1 + epsilon * epsilon, 2}); 228 229 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 230 PointVectorValuePair optimum = 231 optimizer.optimize(new MaxEval(100), 232 problem.getModelFunction(), 233 problem.getModelFunctionJacobian(), 234 problem.getTarget(), 235 new Weight(new double[] { 1, 1, 1, 1, 1, 1 }), 236 new InitialGuess(new double[] { 0, 0, 0, 0, 0, 0 })); 237 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 238 Assert.assertEquals(3, optimum.getPoint()[0], 1e-10); 239 Assert.assertEquals(4, optimum.getPoint()[1], 1e-10); 240 Assert.assertEquals(-1, optimum.getPoint()[2], 1e-10); 241 Assert.assertEquals(-2, optimum.getPoint()[3], 1e-10); 242 Assert.assertEquals(1 + epsilon, optimum.getPoint()[4], 1e-10); 243 Assert.assertEquals(1 - epsilon, optimum.getPoint()[5], 1e-10); 244 } 245 246 @Test(expected=ConvergenceException.class) testNonInvertible()247 public void testNonInvertible() throws Exception { 248 249 LinearProblem problem = new LinearProblem(new double[][] { 250 { 1, 2, -3 }, 251 { 2, 1, 3 }, 252 { -3, 0, -9 } 253 }, new double[] { 1, 1, 1 }); 254 255 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 256 257 optimizer.optimize(new MaxEval(100), 258 problem.getModelFunction(), 259 problem.getModelFunctionJacobian(), 260 problem.getTarget(), 261 new Weight(new double[] { 1, 1, 1 }), 262 new InitialGuess(new double[] { 0, 0, 0 })); 263 } 264 265 @Test testIllConditioned()266 public void testIllConditioned() { 267 LinearProblem problem1 = new LinearProblem(new double[][] { 268 { 10, 7, 8, 7 }, 269 { 7, 5, 6, 5 }, 270 { 8, 6, 10, 9 }, 271 { 7, 5, 9, 10 } 272 }, new double[] { 32, 23, 33, 31 }); 273 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 274 PointVectorValuePair optimum1 = 275 optimizer.optimize(new MaxEval(100), 276 problem1.getModelFunction(), 277 problem1.getModelFunctionJacobian(), 278 problem1.getTarget(), 279 new Weight(new double[] { 1, 1, 1, 1 }), 280 new InitialGuess(new double[] { 0, 1, 2, 3 })); 281 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 282 Assert.assertEquals(1, optimum1.getPoint()[0], 1e-10); 283 Assert.assertEquals(1, optimum1.getPoint()[1], 1e-10); 284 Assert.assertEquals(1, optimum1.getPoint()[2], 1e-10); 285 Assert.assertEquals(1, optimum1.getPoint()[3], 1e-10); 286 287 LinearProblem problem2 = new LinearProblem(new double[][] { 288 { 10.00, 7.00, 8.10, 7.20 }, 289 { 7.08, 5.04, 6.00, 5.00 }, 290 { 8.00, 5.98, 9.89, 9.00 }, 291 { 6.99, 4.99, 9.00, 9.98 } 292 }, new double[] { 32, 23, 33, 31 }); 293 PointVectorValuePair optimum2 = 294 optimizer.optimize(new MaxEval(100), 295 problem2.getModelFunction(), 296 problem2.getModelFunctionJacobian(), 297 problem2.getTarget(), 298 new Weight(new double[] { 1, 1, 1, 1 }), 299 new InitialGuess(new double[] { 0, 1, 2, 3 })); 300 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 301 Assert.assertEquals(-81, optimum2.getPoint()[0], 1e-8); 302 Assert.assertEquals(137, optimum2.getPoint()[1], 1e-8); 303 Assert.assertEquals(-34, optimum2.getPoint()[2], 1e-8); 304 Assert.assertEquals( 22, optimum2.getPoint()[3], 1e-8); 305 } 306 307 @Test testMoreEstimatedParametersSimple()308 public void testMoreEstimatedParametersSimple() { 309 310 LinearProblem problem = new LinearProblem(new double[][] { 311 { 3, 2, 0, 0 }, 312 { 0, 1, -1, 1 }, 313 { 2, 0, 1, 0 } 314 }, new double[] { 7, 3, 5 }); 315 316 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 317 optimizer.optimize(new MaxEval(100), 318 problem.getModelFunction(), 319 problem.getModelFunctionJacobian(), 320 problem.getTarget(), 321 new Weight(new double[] { 1, 1, 1 }), 322 new InitialGuess(new double[] { 7, 6, 5, 4 })); 323 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 324 } 325 326 @Test testMoreEstimatedParametersUnsorted()327 public void testMoreEstimatedParametersUnsorted() { 328 LinearProblem problem = new LinearProblem(new double[][] { 329 { 1, 1, 0, 0, 0, 0 }, 330 { 0, 0, 1, 1, 1, 0 }, 331 { 0, 0, 0, 0, 1, -1 }, 332 { 0, 0, -1, 1, 0, 1 }, 333 { 0, 0, 0, -1, 1, 0 } 334 }, new double[] { 3, 12, -1, 7, 1 }); 335 336 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 337 PointVectorValuePair optimum = 338 optimizer.optimize(new MaxEval(100), 339 problem.getModelFunction(), 340 problem.getModelFunctionJacobian(), 341 problem.getTarget(), 342 new Weight(new double[] { 1, 1, 1, 1, 1 }), 343 new InitialGuess(new double[] { 2, 2, 2, 2, 2, 2 })); 344 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 345 Assert.assertEquals(3, optimum.getPointRef()[2], 1e-10); 346 Assert.assertEquals(4, optimum.getPointRef()[3], 1e-10); 347 Assert.assertEquals(5, optimum.getPointRef()[4], 1e-10); 348 Assert.assertEquals(6, optimum.getPointRef()[5], 1e-10); 349 } 350 351 @Test testRedundantEquations()352 public void testRedundantEquations() { 353 LinearProblem problem = new LinearProblem(new double[][] { 354 { 1, 1 }, 355 { 1, -1 }, 356 { 1, 3 } 357 }, new double[] { 3, 1, 5 }); 358 359 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 360 PointVectorValuePair optimum = 361 optimizer.optimize(new MaxEval(100), 362 problem.getModelFunction(), 363 problem.getModelFunctionJacobian(), 364 problem.getTarget(), 365 new Weight(new double[] { 1, 1, 1 }), 366 new InitialGuess(new double[] { 1, 1 })); 367 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 368 Assert.assertEquals(2, optimum.getPointRef()[0], 1e-10); 369 Assert.assertEquals(1, optimum.getPointRef()[1], 1e-10); 370 } 371 372 @Test testInconsistentEquations()373 public void testInconsistentEquations() { 374 LinearProblem problem = new LinearProblem(new double[][] { 375 { 1, 1 }, 376 { 1, -1 }, 377 { 1, 3 } 378 }, new double[] { 3, 1, 4 }); 379 380 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 381 optimizer.optimize(new MaxEval(100), 382 problem.getModelFunction(), 383 problem.getModelFunctionJacobian(), 384 problem.getTarget(), 385 new Weight(new double[] { 1, 1, 1 }), 386 new InitialGuess(new double[] { 1, 1 })); 387 Assert.assertTrue(optimizer.getRMS() > 0.1); 388 } 389 390 @Test(expected=DimensionMismatchException.class) testInconsistentSizes1()391 public void testInconsistentSizes1() { 392 LinearProblem problem 393 = new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, 394 new double[] { -1, 1 }); 395 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 396 PointVectorValuePair optimum = 397 optimizer.optimize(new MaxEval(100), 398 problem.getModelFunction(), 399 problem.getModelFunctionJacobian(), 400 problem.getTarget(), 401 new Weight(new double[] { 1, 1 }), 402 new InitialGuess(new double[] { 0, 0 })); 403 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 404 Assert.assertEquals(-1, optimum.getPoint()[0], 1e-10); 405 Assert.assertEquals(1, optimum.getPoint()[1], 1e-10); 406 407 optimizer.optimize(new MaxEval(100), 408 problem.getModelFunction(), 409 problem.getModelFunctionJacobian(), 410 problem.getTarget(), 411 new Weight(new double[] { 1 }), 412 new InitialGuess(new double[] { 0, 0 })); 413 } 414 415 @Test(expected=DimensionMismatchException.class) testInconsistentSizes2()416 public void testInconsistentSizes2() { 417 LinearProblem problem 418 = new LinearProblem(new double[][] { { 1, 0 }, { 0, 1 } }, 419 new double[] { -1, 1 }); 420 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 421 PointVectorValuePair optimum 422 = optimizer.optimize(new MaxEval(100), 423 problem.getModelFunction(), 424 problem.getModelFunctionJacobian(), 425 problem.getTarget(), 426 new Weight(new double[] { 1, 1 }), 427 new InitialGuess(new double[] { 0, 0 })); 428 Assert.assertEquals(0, optimizer.getRMS(), 1e-10); 429 Assert.assertEquals(-1, optimum.getPoint()[0], 1e-10); 430 Assert.assertEquals(1, optimum.getPoint()[1], 1e-10); 431 432 optimizer.optimize(new MaxEval(100), 433 problem.getModelFunction(), 434 problem.getModelFunctionJacobian(), 435 new Target(new double[] { 1 }), 436 new Weight(new double[] { 1 }), 437 new InitialGuess(new double[] { 0, 0 })); 438 } 439 440 @Test testCircleFitting()441 public void testCircleFitting() { 442 CircleVectorial circle = new CircleVectorial(); 443 circle.addPoint( 30, 68); 444 circle.addPoint( 50, -6); 445 circle.addPoint(110, -20); 446 circle.addPoint( 35, 15); 447 circle.addPoint( 45, 97); 448 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 449 PointVectorValuePair optimum 450 = optimizer.optimize(new MaxEval(100), 451 circle.getModelFunction(), 452 circle.getModelFunctionJacobian(), 453 new Target(new double[] { 0, 0, 0, 0, 0 }), 454 new Weight(new double[] { 1, 1, 1, 1, 1 }), 455 new InitialGuess(new double[] { 98.680, 47.345 })); 456 Assert.assertTrue(optimizer.getEvaluations() < 10); 457 double rms = optimizer.getRMS(); 458 Assert.assertEquals(1.768262623567235, FastMath.sqrt(circle.getN()) * rms, 1e-10); 459 Vector2D center = new Vector2D(optimum.getPointRef()[0], optimum.getPointRef()[1]); 460 Assert.assertEquals(69.96016176931406, circle.getRadius(center), 1e-6); 461 Assert.assertEquals(96.07590211815305, center.getX(), 1e-6); 462 Assert.assertEquals(48.13516790438953, center.getY(), 1e-6); 463 double[][] cov = optimizer.computeCovariances(optimum.getPoint(), 1e-14); 464 Assert.assertEquals(1.839, cov[0][0], 0.001); 465 Assert.assertEquals(0.731, cov[0][1], 0.001); 466 Assert.assertEquals(cov[0][1], cov[1][0], 1e-14); 467 Assert.assertEquals(0.786, cov[1][1], 0.001); 468 469 // add perfect measurements and check errors are reduced 470 double r = circle.getRadius(center); 471 for (double d= 0; d < 2 * FastMath.PI; d += 0.01) { 472 circle.addPoint(center.getX() + r * FastMath.cos(d), center.getY() + r * FastMath.sin(d)); 473 } 474 double[] target = new double[circle.getN()]; 475 Arrays.fill(target, 0); 476 double[] weights = new double[circle.getN()]; 477 Arrays.fill(weights, 2); 478 optimum = optimizer.optimize(new MaxEval(100), 479 circle.getModelFunction(), 480 circle.getModelFunctionJacobian(), 481 new Target(target), 482 new Weight(weights), 483 new InitialGuess(new double[] { 98.680, 47.345 })); 484 cov = optimizer.computeCovariances(optimum.getPoint(), 1e-14); 485 Assert.assertEquals(0.0016, cov[0][0], 0.001); 486 Assert.assertEquals(3.2e-7, cov[0][1], 1e-9); 487 Assert.assertEquals(cov[0][1], cov[1][0], 1e-14); 488 Assert.assertEquals(0.0016, cov[1][1], 0.001); 489 } 490 491 @Test 492 public void testCircleFittingBadInit() { 493 CircleVectorial circle = new CircleVectorial(); 494 double[][] points = circlePoints; 495 double[] target = new double[points.length]; 496 Arrays.fill(target, 0); 497 double[] weights = new double[points.length]; 498 Arrays.fill(weights, 2); 499 for (int i = 0; i < points.length; ++i) { 500 circle.addPoint(points[i][0], points[i][1]); 501 } 502 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 503 PointVectorValuePair optimum 504 = optimizer.optimize(new MaxEval(100), 505 circle.getModelFunction(), 506 circle.getModelFunctionJacobian(), 507 new Target(target), 508 new Weight(weights), 509 new InitialGuess(new double[] { -12, -12 })); 510 Vector2D center = new Vector2D(optimum.getPointRef()[0], optimum.getPointRef()[1]); 511 Assert.assertTrue(optimizer.getEvaluations() < 25); 512 Assert.assertEquals( 0.043, optimizer.getRMS(), 1e-3); 513 Assert.assertEquals( 0.292235, circle.getRadius(center), 1e-6); 514 Assert.assertEquals(-0.151738, center.getX(), 1e-6); 515 Assert.assertEquals( 0.2075001, center.getY(), 1e-6); 516 } 517 518 @Test 519 public void testCircleFittingGoodInit() { 520 CircleVectorial circle = new CircleVectorial(); 521 double[][] points = circlePoints; 522 double[] target = new double[points.length]; 523 Arrays.fill(target, 0); 524 double[] weights = new double[points.length]; 525 Arrays.fill(weights, 2); 526 for (int i = 0; i < points.length; ++i) { 527 circle.addPoint(points[i][0], points[i][1]); 528 } 529 AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 530 PointVectorValuePair optimum = 531 optimizer.optimize(new MaxEval(100), 532 circle.getModelFunction(), 533 circle.getModelFunctionJacobian(), 534 new Target(target), 535 new Weight(weights), 536 new InitialGuess(new double[] { 0, 0 })); 537 Assert.assertEquals(-0.1517383071957963, optimum.getPointRef()[0], 1e-6); 538 Assert.assertEquals(0.2074999736353867, optimum.getPointRef()[1], 1e-6); 539 Assert.assertEquals(0.04268731682389561, optimizer.getRMS(), 1e-8); 540 } 541 542 private final double[][] circlePoints = new double[][] { 543 {-0.312967, 0.072366}, {-0.339248, 0.132965}, {-0.379780, 0.202724}, 544 {-0.390426, 0.260487}, {-0.361212, 0.328325}, {-0.346039, 0.392619}, 545 {-0.280579, 0.444306}, {-0.216035, 0.470009}, {-0.149127, 0.493832}, 546 {-0.075133, 0.483271}, {-0.007759, 0.452680}, { 0.060071, 0.410235}, 547 { 0.103037, 0.341076}, { 0.118438, 0.273884}, { 0.131293, 0.192201}, 548 { 0.115869, 0.129797}, { 0.072223, 0.058396}, { 0.022884, 0.000718}, 549 {-0.053355, -0.020405}, {-0.123584, -0.032451}, {-0.216248, -0.032862}, 550 {-0.278592, -0.005008}, {-0.337655, 0.056658}, {-0.385899, 0.112526}, 551 {-0.405517, 0.186957}, {-0.415374, 0.262071}, {-0.387482, 0.343398}, 552 {-0.347322, 0.397943}, {-0.287623, 0.458425}, {-0.223502, 0.475513}, 553 {-0.135352, 0.478186}, {-0.061221, 0.483371}, { 0.003711, 0.422737}, 554 { 0.065054, 0.375830}, { 0.108108, 0.297099}, { 0.123882, 0.222850}, 555 { 0.117729, 0.134382}, { 0.085195, 0.056820}, { 0.029800, -0.019138}, 556 {-0.027520, -0.072374}, {-0.102268, -0.091555}, {-0.200299, -0.106578}, 557 {-0.292731, -0.091473}, {-0.356288, -0.051108}, {-0.420561, 0.014926}, 558 {-0.471036, 0.074716}, {-0.488638, 0.182508}, {-0.485990, 0.254068}, 559 {-0.463943, 0.338438}, {-0.406453, 0.404704}, {-0.334287, 0.466119}, 560 {-0.254244, 0.503188}, {-0.161548, 0.495769}, {-0.075733, 0.495560}, 561 { 0.001375, 0.434937}, { 0.082787, 0.385806}, { 0.115490, 0.323807}, 562 { 0.141089, 0.223450}, { 0.138693, 0.131703}, { 0.126415, 0.049174}, 563 { 0.066518, -0.010217}, {-0.005184, -0.070647}, {-0.080985, -0.103635}, 564 {-0.177377, -0.116887}, {-0.260628, -0.100258}, {-0.335756, -0.056251}, 565 {-0.405195, -0.000895}, {-0.444937, 0.085456}, {-0.484357, 0.175597}, 566 {-0.472453, 0.248681}, {-0.438580, 0.347463}, {-0.402304, 0.422428}, 567 {-0.326777, 0.479438}, {-0.247797, 0.505581}, {-0.152676, 0.519380}, 568 {-0.071754, 0.516264}, { 0.015942, 0.472802}, { 0.076608, 0.419077}, 569 { 0.127673, 0.330264}, { 0.159951, 0.262150}, { 0.153530, 0.172681}, 570 { 0.140653, 0.089229}, { 0.078666, 0.024981}, { 0.023807, -0.037022}, 571 {-0.048837, -0.077056}, {-0.127729, -0.075338}, {-0.221271, -0.067526} 572 }; 573 574 public void doTestStRD(final StatisticalReferenceDataset dataset, 575 final double errParams, 576 final double errParamsSd) { 577 final AbstractLeastSquaresOptimizer optimizer = createOptimizer(); 578 final double[] w = new double[dataset.getNumObservations()]; 579 Arrays.fill(w, 1); 580 581 final double[][] data = dataset.getData(); 582 final double[] initial = dataset.getStartingPoint(0); 583 final StatisticalReferenceDataset.LeastSquaresProblem problem = dataset.getLeastSquaresProblem(); 584 final PointVectorValuePair optimum 585 = optimizer.optimize(new MaxEval(100), 586 problem.getModelFunction(), 587 problem.getModelFunctionJacobian(), 588 new Target(data[1]), 589 new Weight(w), 590 new InitialGuess(initial)); 591 592 final double[] actual = optimum.getPoint(); 593 for (int i = 0; i < actual.length; i++) { 594 double expected = dataset.getParameter(i); 595 double delta = FastMath.abs(errParams * expected); 596 Assert.assertEquals(dataset.getName() + ", param #" + i, 597 expected, actual[i], delta); 598 } 599 } 600 601 @Test 602 public void testKirby2() throws IOException { 603 doTestStRD(StatisticalReferenceDatasetFactory.createKirby2(), 1E-7, 1E-7); 604 } 605 606 @Test 607 public void testHahn1() throws IOException { 608 doTestStRD(StatisticalReferenceDatasetFactory.createHahn1(), 1E-7, 1E-4); 609 } 610 611 static class LinearProblem { 612 private final RealMatrix factors; 613 private final double[] target; 614 615 public LinearProblem(double[][] factors, double[] target) { 616 this.factors = new BlockRealMatrix(factors); 617 this.target = target; 618 } 619 620 public Target getTarget() { 621 return new Target(target); 622 } 623 624 public ModelFunction getModelFunction() { 625 return new ModelFunction(new MultivariateVectorFunction() { 626 public double[] value(double[] params) { 627 return factors.operate(params); 628 } 629 }); 630 } 631 632 public ModelFunctionJacobian getModelFunctionJacobian() { 633 return new ModelFunctionJacobian(new MultivariateMatrixFunction() { 634 public double[][] value(double[] params) { 635 return factors.getData(); 636 } 637 }); 638 } 639 } 640 } 641