1 // Licensed to the Apache Software Foundation (ASF) under one 2 // or more contributor license agreements. See the NOTICE file 3 // distributed with this work for additional information 4 // regarding copyright ownership. The ASF licenses this file 5 // to you under the Apache License, Version 2.0 (the 6 // "License"); you may not use this file except in compliance 7 // with 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, 12 // software distributed under the License is distributed on an 13 // "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY 14 // KIND, either express or implied. See the License for the 15 // specific language governing permissions and limitations 16 // under the License. 17 18 package org.apache.commons.math3.optimization.fitting; 19 20 import java.util.Random; 21 22 import org.apache.commons.math3.analysis.polynomials.PolynomialFunction; 23 import org.apache.commons.math3.analysis.polynomials.PolynomialFunction.Parametric; 24 import org.apache.commons.math3.exception.ConvergenceException; 25 import org.apache.commons.math3.exception.TooManyEvaluationsException; 26 import org.apache.commons.math3.optimization.DifferentiableMultivariateVectorOptimizer; 27 import org.apache.commons.math3.optimization.general.GaussNewtonOptimizer; 28 import org.apache.commons.math3.optimization.general.LevenbergMarquardtOptimizer; 29 import org.apache.commons.math3.optimization.SimpleVectorValueChecker; 30 import org.apache.commons.math3.util.FastMath; 31 import org.apache.commons.math3.distribution.RealDistribution; 32 import org.apache.commons.math3.distribution.UniformRealDistribution; 33 import org.apache.commons.math3.TestUtils; 34 import org.junit.Test; 35 import org.junit.Assert; 36 37 /** 38 * Test for class {@link CurveFitter} where the function to fit is a 39 * polynomial. 40 */ 41 @Deprecated 42 public class PolynomialFitterTest { 43 @Test testFit()44 public void testFit() { 45 final RealDistribution rng = new UniformRealDistribution(-100, 100); 46 rng.reseedRandomGenerator(64925784252L); 47 48 final LevenbergMarquardtOptimizer optim = new LevenbergMarquardtOptimizer(); 49 final PolynomialFitter fitter = new PolynomialFitter(optim); 50 final double[] coeff = { 12.9, -3.4, 2.1 }; // 12.9 - 3.4 x + 2.1 x^2 51 final PolynomialFunction f = new PolynomialFunction(coeff); 52 53 // Collect data from a known polynomial. 54 for (int i = 0; i < 100; i++) { 55 final double x = rng.sample(); 56 fitter.addObservedPoint(x, f.value(x)); 57 } 58 59 // Start fit from initial guesses that are far from the optimal values. 60 final double[] best = fitter.fit(new double[] { -1e-20, 3e15, -5e25 }); 61 62 TestUtils.assertEquals("best != coeff", coeff, best, 1e-12); 63 } 64 65 @Test testNoError()66 public void testNoError() { 67 Random randomizer = new Random(64925784252l); 68 for (int degree = 1; degree < 10; ++degree) { 69 PolynomialFunction p = buildRandomPolynomial(degree, randomizer); 70 71 PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer()); 72 for (int i = 0; i <= degree; ++i) { 73 fitter.addObservedPoint(1.0, i, p.value(i)); 74 } 75 76 final double[] init = new double[degree + 1]; 77 PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init)); 78 79 for (double x = -1.0; x < 1.0; x += 0.01) { 80 double error = FastMath.abs(p.value(x) - fitted.value(x)) / 81 (1.0 + FastMath.abs(p.value(x))); 82 Assert.assertEquals(0.0, error, 1.0e-6); 83 } 84 } 85 } 86 87 @Test testSmallError()88 public void testSmallError() { 89 Random randomizer = new Random(53882150042l); 90 double maxError = 0; 91 for (int degree = 0; degree < 10; ++degree) { 92 PolynomialFunction p = buildRandomPolynomial(degree, randomizer); 93 94 PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer()); 95 for (double x = -1.0; x < 1.0; x += 0.01) { 96 fitter.addObservedPoint(1.0, x, 97 p.value(x) + 0.1 * randomizer.nextGaussian()); 98 } 99 100 final double[] init = new double[degree + 1]; 101 PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init)); 102 103 for (double x = -1.0; x < 1.0; x += 0.01) { 104 double error = FastMath.abs(p.value(x) - fitted.value(x)) / 105 (1.0 + FastMath.abs(p.value(x))); 106 maxError = FastMath.max(maxError, error); 107 Assert.assertTrue(FastMath.abs(error) < 0.1); 108 } 109 } 110 Assert.assertTrue(maxError > 0.01); 111 } 112 113 @Test testMath798()114 public void testMath798() { 115 final double tol = 1e-14; 116 final SimpleVectorValueChecker checker = new SimpleVectorValueChecker(tol, tol); 117 final double[] init = new double[] { 0, 0 }; 118 final int maxEval = 3; 119 120 final double[] lm = doMath798(new LevenbergMarquardtOptimizer(checker), maxEval, init); 121 final double[] gn = doMath798(new GaussNewtonOptimizer(checker), maxEval, init); 122 123 for (int i = 0; i <= 1; i++) { 124 Assert.assertEquals(lm[i], gn[i], tol); 125 } 126 } 127 128 /** 129 * This test shows that the user can set the maximum number of iterations 130 * to avoid running for too long. 131 * But in the test case, the real problem is that the tolerance is way too 132 * stringent. 133 */ 134 @Test(expected=TooManyEvaluationsException.class) testMath798WithToleranceTooLow()135 public void testMath798WithToleranceTooLow() { 136 final double tol = 1e-100; 137 final SimpleVectorValueChecker checker = new SimpleVectorValueChecker(tol, tol); 138 final double[] init = new double[] { 0, 0 }; 139 final int maxEval = 10000; // Trying hard to fit. 140 141 doMath798(new GaussNewtonOptimizer(checker), maxEval, init); 142 } 143 144 /** 145 * This test shows that the user can set the maximum number of iterations 146 * to avoid running for too long. 147 * Even if the real problem is that the tolerance is way too stringent, it 148 * is possible to get the best solution so far, i.e. a checker will return 149 * the point when the maximum iteration count has been reached. 150 */ 151 @Test testMath798WithToleranceTooLowButNoException()152 public void testMath798WithToleranceTooLowButNoException() { 153 final double tol = 1e-100; 154 final double[] init = new double[] { 0, 0 }; 155 final int maxEval = 10000; // Trying hard to fit. 156 final SimpleVectorValueChecker checker = new SimpleVectorValueChecker(tol, tol, maxEval); 157 158 final double[] lm = doMath798(new LevenbergMarquardtOptimizer(checker), maxEval, init); 159 final double[] gn = doMath798(new GaussNewtonOptimizer(checker), maxEval, init); 160 161 for (int i = 0; i <= 1; i++) { 162 Assert.assertEquals(lm[i], gn[i], 1e-15); 163 } 164 } 165 166 /** 167 * @param optimizer Optimizer. 168 * @param maxEval Maximum number of function evaluations. 169 * @param init First guess. 170 * @return the solution found by the given optimizer. 171 */ doMath798(DifferentiableMultivariateVectorOptimizer optimizer, int maxEval, double[] init)172 private double[] doMath798(DifferentiableMultivariateVectorOptimizer optimizer, 173 int maxEval, 174 double[] init) { 175 final CurveFitter<Parametric> fitter = new CurveFitter<Parametric>(optimizer); 176 177 fitter.addObservedPoint(-0.2, -7.12442E-13); 178 fitter.addObservedPoint(-0.199, -4.33397E-13); 179 fitter.addObservedPoint(-0.198, -2.823E-13); 180 fitter.addObservedPoint(-0.197, -1.40405E-13); 181 fitter.addObservedPoint(-0.196, -7.80821E-15); 182 fitter.addObservedPoint(-0.195, 6.20484E-14); 183 fitter.addObservedPoint(-0.194, 7.24673E-14); 184 fitter.addObservedPoint(-0.193, 1.47152E-13); 185 fitter.addObservedPoint(-0.192, 1.9629E-13); 186 fitter.addObservedPoint(-0.191, 2.12038E-13); 187 fitter.addObservedPoint(-0.19, 2.46906E-13); 188 fitter.addObservedPoint(-0.189, 2.77495E-13); 189 fitter.addObservedPoint(-0.188, 2.51281E-13); 190 fitter.addObservedPoint(-0.187, 2.64001E-13); 191 fitter.addObservedPoint(-0.186, 2.8882E-13); 192 fitter.addObservedPoint(-0.185, 3.13604E-13); 193 fitter.addObservedPoint(-0.184, 3.14248E-13); 194 fitter.addObservedPoint(-0.183, 3.1172E-13); 195 fitter.addObservedPoint(-0.182, 3.12912E-13); 196 fitter.addObservedPoint(-0.181, 3.06761E-13); 197 fitter.addObservedPoint(-0.18, 2.8559E-13); 198 fitter.addObservedPoint(-0.179, 2.86806E-13); 199 fitter.addObservedPoint(-0.178, 2.985E-13); 200 fitter.addObservedPoint(-0.177, 2.67148E-13); 201 fitter.addObservedPoint(-0.176, 2.94173E-13); 202 fitter.addObservedPoint(-0.175, 3.27528E-13); 203 fitter.addObservedPoint(-0.174, 3.33858E-13); 204 fitter.addObservedPoint(-0.173, 2.97511E-13); 205 fitter.addObservedPoint(-0.172, 2.8615E-13); 206 fitter.addObservedPoint(-0.171, 2.84624E-13); 207 208 final double[] coeff = fitter.fit(maxEval, 209 new PolynomialFunction.Parametric(), 210 init); 211 return coeff; 212 } 213 214 @Test testRedundantSolvable()215 public void testRedundantSolvable() { 216 // Levenberg-Marquardt should handle redundant information gracefully 217 checkUnsolvableProblem(new LevenbergMarquardtOptimizer(), true); 218 } 219 220 @Test testRedundantUnsolvable()221 public void testRedundantUnsolvable() { 222 // Gauss-Newton should not be able to solve redundant information 223 checkUnsolvableProblem(new GaussNewtonOptimizer(true, new SimpleVectorValueChecker(1e-15, 1e-15)), false); 224 } 225 226 @Test testLargeSample()227 public void testLargeSample() { 228 Random randomizer = new Random(0x5551480dca5b369bl); 229 double maxError = 0; 230 for (int degree = 0; degree < 10; ++degree) { 231 PolynomialFunction p = buildRandomPolynomial(degree, randomizer); 232 233 PolynomialFitter fitter = new PolynomialFitter(new LevenbergMarquardtOptimizer()); 234 for (int i = 0; i < 40000; ++i) { 235 double x = -1.0 + i / 20000.0; 236 fitter.addObservedPoint(1.0, x, 237 p.value(x) + 0.1 * randomizer.nextGaussian()); 238 } 239 240 final double[] init = new double[degree + 1]; 241 PolynomialFunction fitted = new PolynomialFunction(fitter.fit(init)); 242 243 for (double x = -1.0; x < 1.0; x += 0.01) { 244 double error = FastMath.abs(p.value(x) - fitted.value(x)) / 245 (1.0 + FastMath.abs(p.value(x))); 246 maxError = FastMath.max(maxError, error); 247 Assert.assertTrue(FastMath.abs(error) < 0.01); 248 } 249 } 250 Assert.assertTrue(maxError > 0.001); 251 } 252 checkUnsolvableProblem(DifferentiableMultivariateVectorOptimizer optimizer, boolean solvable)253 private void checkUnsolvableProblem(DifferentiableMultivariateVectorOptimizer optimizer, 254 boolean solvable) { 255 Random randomizer = new Random(1248788532l); 256 for (int degree = 0; degree < 10; ++degree) { 257 PolynomialFunction p = buildRandomPolynomial(degree, randomizer); 258 259 PolynomialFitter fitter = new PolynomialFitter(optimizer); 260 261 // reusing the same point over and over again does not bring 262 // information, the problem cannot be solved in this case for 263 // degrees greater than 1 (but one point is sufficient for 264 // degree 0) 265 for (double x = -1.0; x < 1.0; x += 0.01) { 266 fitter.addObservedPoint(1.0, 0.0, p.value(0.0)); 267 } 268 269 try { 270 final double[] init = new double[degree + 1]; 271 fitter.fit(init); 272 Assert.assertTrue(solvable || (degree == 0)); 273 } catch(ConvergenceException e) { 274 Assert.assertTrue((! solvable) && (degree > 0)); 275 } 276 } 277 } 278 buildRandomPolynomial(int degree, Random randomizer)279 private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) { 280 final double[] coefficients = new double[degree + 1]; 281 for (int i = 0; i <= degree; ++i) { 282 coefficients[i] = randomizer.nextGaussian(); 283 } 284 return new PolynomialFunction(coefficients); 285 } 286 } 287