View Javadoc
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  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.fitting;
23  
24  import java.util.Random;
25  
26  import org.hipparchus.UnitTestUtils;
27  import org.hipparchus.analysis.polynomials.PolynomialFunction;
28  import org.hipparchus.exception.MathIllegalStateException;
29  import org.hipparchus.random.RandomDataGenerator;
30  import org.hipparchus.util.FastMath;
31  import org.junit.Assert;
32  import org.junit.Test;
33  
34  /**
35   * Test for class {@link PolynomialCurveFitter}.
36   */
37  public class PolynomialCurveFitterTest {
38      @Test
39      public void testFit() {
40          final RandomDataGenerator randomDataGenerator = new RandomDataGenerator(64925784252L);
41  
42          final double[] coeff = { 12.9, -3.4, 2.1 }; // 12.9 - 3.4 x + 2.1 x^2
43          final PolynomialFunction f = new PolynomialFunction(coeff);
44  
45          // Collect data from a known polynomial.
46          final WeightedObservedPoints obs = new WeightedObservedPoints();
47          for (int i = 0; i < 100; i++) {
48              final double x = randomDataGenerator.nextUniform(-100, 100);
49              obs.add(x, f.value(x));
50          }
51  
52          // Start fit from initial guesses that are far from the optimal values.
53          final PolynomialCurveFitter fitter
54              = PolynomialCurveFitter.create(0).withStartPoint(new double[] { -1e-20, 3e15, -5e25 });
55          final double[] best = fitter.fit(obs.toList());
56  
57          UnitTestUtils.assertEquals("best != coeff", coeff, best, 1e-12);
58      }
59  
60      @Test
61      public void testNoError() {
62          final Random randomizer = new Random(64925784252l);
63          for (int degree = 1; degree < 10; ++degree) {
64              final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
65              final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
66  
67              final WeightedObservedPoints obs = new WeightedObservedPoints();
68              for (int i = 0; i <= degree; ++i) {
69                  obs.add(1.0, i, p.value(i));
70              }
71  
72              final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
73  
74              for (double x = -1.0; x < 1.0; x += 0.01) {
75                  final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
76                      (1.0 + FastMath.abs(p.value(x)));
77                  Assert.assertEquals(0.0, error, 1.0e-6);
78              }
79          }
80      }
81  
82      @Test
83      public void testSmallError() {
84          final Random randomizer = new Random(53882150042l);
85          double maxError = 0;
86          for (int degree = 0; degree < 10; ++degree) {
87              final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
88              final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
89  
90              final WeightedObservedPoints obs = new WeightedObservedPoints();
91              for (double x = -1.0; x < 1.0; x += 0.01) {
92                  obs.add(1.0, x, p.value(x) + 0.1 * randomizer.nextGaussian());
93              }
94  
95              final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
96  
97              for (double x = -1.0; x < 1.0; x += 0.01) {
98                  final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
99                      (1.0 + FastMath.abs(p.value(x)));
100                 maxError = FastMath.max(maxError, error);
101                 Assert.assertTrue(FastMath.abs(error) < 0.1);
102             }
103         }
104         Assert.assertTrue(maxError > 0.01);
105     }
106 
107     @Test
108     public void testRedundantSolvable() {
109         // Levenberg-Marquardt should handle redundant information gracefully
110         checkUnsolvableProblem(true);
111     }
112 
113     @Test
114     public void testLargeSample() {
115         final Random randomizer = new Random(0x5551480dca5b369bl);
116         double maxError = 0;
117         for (int degree = 0; degree < 10; ++degree) {
118             final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
119             final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
120 
121             final WeightedObservedPoints obs = new WeightedObservedPoints();
122             for (int i = 0; i < 40000; ++i) {
123                 final double x = -1.0 + i / 20000.0;
124                 obs.add(1.0, x, p.value(x) + 0.1 * randomizer.nextGaussian());
125             }
126 
127             final PolynomialFunction fitted = new PolynomialFunction(fitter.fit(obs.toList()));
128             for (double x = -1.0; x < 1.0; x += 0.01) {
129                 final double error = FastMath.abs(p.value(x) - fitted.value(x)) /
130                     (1.0 + FastMath.abs(p.value(x)));
131                 maxError = FastMath.max(maxError, error);
132                 Assert.assertTrue(FastMath.abs(error) < 0.01);
133             }
134         }
135         Assert.assertTrue(maxError > 0.001);
136     }
137 
138     private void checkUnsolvableProblem(boolean solvable) {
139         final Random randomizer = new Random(1248788532l);
140 
141         for (int degree = 0; degree < 10; ++degree) {
142             final PolynomialFunction p = buildRandomPolynomial(degree, randomizer);
143             final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(degree);
144             final WeightedObservedPoints obs = new WeightedObservedPoints();
145             // reusing the same point over and over again does not bring
146             // information, the problem cannot be solved in this case for
147             // degrees greater than 1 (but one point is sufficient for
148             // degree 0)
149             for (double x = -1.0; x < 1.0; x += 0.01) {
150                 obs.add(1.0, 0.0, p.value(0.0));
151             }
152 
153             try {
154                 fitter.fit(obs.toList());
155                 Assert.assertTrue(solvable || (degree == 0));
156             } catch(MathIllegalStateException e) {
157                 Assert.assertTrue((! solvable) && (degree > 0));
158             }
159         }
160     }
161 
162     private PolynomialFunction buildRandomPolynomial(int degree, Random randomizer) {
163         final double[] coefficients = new double[degree + 1];
164         for (int i = 0; i <= degree; ++i) {
165             coefficients[i] = randomizer.nextGaussian();
166         }
167         return new PolynomialFunction(coefficients);
168     }
169 }