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