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.fitting;
18 
19 import java.util.ArrayList;
20 import java.util.List;
21 import java.util.Random;
22 
23 import org.apache.commons.math3.analysis.function.HarmonicOscillator;
24 import org.apache.commons.math3.exception.MathIllegalStateException;
25 import org.apache.commons.math3.exception.NumberIsTooSmallException;
26 import org.apache.commons.math3.util.FastMath;
27 import org.apache.commons.math3.util.MathUtils;
28 import org.junit.Assert;
29 import org.junit.Test;
30 
31 public class HarmonicCurveFitterTest {
32     /**
33      * Zero points is not enough observed points.
34      */
35     @Test(expected=NumberIsTooSmallException.class)
testPreconditions1()36     public void testPreconditions1() {
37         HarmonicCurveFitter.create().fit(new WeightedObservedPoints().toList());
38     }
39 
40     @Test
testNoError()41     public void testNoError() {
42         final double a = 0.2;
43         final double w = 3.4;
44         final double p = 4.1;
45         final HarmonicOscillator f = new HarmonicOscillator(a, w, p);
46 
47         final WeightedObservedPoints points = new WeightedObservedPoints();
48         for (double x = 0.0; x < 1.3; x += 0.01) {
49             points.add(1, x, f.value(x));
50         }
51 
52         final HarmonicCurveFitter fitter = HarmonicCurveFitter.create();
53         final double[] fitted = fitter.fit(points.toList());
54         Assert.assertEquals(a, fitted[0], 1.0e-13);
55         Assert.assertEquals(w, fitted[1], 1.0e-13);
56         Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1e-13);
57 
58         final HarmonicOscillator ff = new HarmonicOscillator(fitted[0], fitted[1], fitted[2]);
59         for (double x = -1.0; x < 1.0; x += 0.01) {
60             Assert.assertTrue(FastMath.abs(f.value(x) - ff.value(x)) < 1e-13);
61         }
62     }
63 
64     @Test
65     public void test1PercentError() {
66         final Random randomizer = new Random(64925784252L);
67         final double a = 0.2;
68         final double w = 3.4;
69         final double p = 4.1;
70         final HarmonicOscillator f = new HarmonicOscillator(a, w, p);
71 
72         final WeightedObservedPoints points = new WeightedObservedPoints();
73         for (double x = 0.0; x < 10.0; x += 0.1) {
74             points.add(1, x, f.value(x) + 0.01 * randomizer.nextGaussian());
75         }
76 
77         final HarmonicCurveFitter fitter = HarmonicCurveFitter.create();
78         final double[] fitted = fitter.fit(points.toList());
79         Assert.assertEquals(a, fitted[0], 7.6e-4);
80         Assert.assertEquals(w, fitted[1], 2.7e-3);
81         Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.3e-2);
82     }
83 
84     @Test
85     public void testTinyVariationsData() {
86         final Random randomizer = new Random(64925784252L);
87 
88         final WeightedObservedPoints points = new WeightedObservedPoints();
89         for (double x = 0.0; x < 10.0; x += 0.1) {
90             points.add(1, x, 1e-7 * randomizer.nextGaussian());
91         }
92 
93         final HarmonicCurveFitter fitter = HarmonicCurveFitter.create();
94         fitter.fit(points.toList());
95 
96         // This test serves to cover the part of the code of "guessAOmega"
97         // when the algorithm using integrals fails.
98     }
99 
100     @Test
101     public void testInitialGuess() {
102         final Random randomizer = new Random(45314242L);
103         final double a = 0.2;
104         final double w = 3.4;
105         final double p = 4.1;
106         final HarmonicOscillator f = new HarmonicOscillator(a, w, p);
107 
108         final WeightedObservedPoints points = new WeightedObservedPoints();
109         for (double x = 0.0; x < 10.0; x += 0.1) {
110             points.add(1, x, f.value(x) + 0.01 * randomizer.nextGaussian());
111         }
112 
113         final HarmonicCurveFitter fitter = HarmonicCurveFitter.create()
114             .withStartPoint(new double[] { 0.15, 3.6, 4.5 });
115         final double[] fitted = fitter.fit(points.toList());
116         Assert.assertEquals(a, fitted[0], 1.2e-3);
117         Assert.assertEquals(w, fitted[1], 3.3e-3);
118         Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.7e-2);
119     }
120 
121     @Test
122     public void testUnsorted() {
123         Random randomizer = new Random(64925784252L);
124         final double a = 0.2;
125         final double w = 3.4;
126         final double p = 4.1;
127         final HarmonicOscillator f = new HarmonicOscillator(a, w, p);
128 
129         // Build a regularly spaced array of measurements.
130         final int size = 100;
131         final double[] xTab = new double[size];
132         final double[] yTab = new double[size];
133         for (int i = 0; i < size; i++) {
134             xTab[i] = 0.1 * i;
135             yTab[i] = f.value(xTab[i]) + 0.01 * randomizer.nextGaussian();
136         }
137 
138         // shake it
139         for (int i = 0; i < size; i++) {
140             int i1 = randomizer.nextInt(size);
141             int i2 = randomizer.nextInt(size);
142             double xTmp = xTab[i1];
143             double yTmp = yTab[i1];
144             xTab[i1] = xTab[i2];
145             yTab[i1] = yTab[i2];
146             xTab[i2] = xTmp;
147             yTab[i2] = yTmp;
148         }
149 
150         // Pass it to the fitter.
151         final WeightedObservedPoints points = new WeightedObservedPoints();
152         for (int i = 0; i < size; ++i) {
153             points.add(1, xTab[i], yTab[i]);
154         }
155 
156         final HarmonicCurveFitter fitter = HarmonicCurveFitter.create();
157         final double[] fitted = fitter.fit(points.toList());
158         Assert.assertEquals(a, fitted[0], 7.6e-4);
159         Assert.assertEquals(w, fitted[1], 3.5e-3);
160         Assert.assertEquals(p, MathUtils.normalizeAngle(fitted[2], p), 1.5e-2);
161     }
162 
163     @Test(expected=MathIllegalStateException.class)
164     public void testMath844() {
165         final double[] y = { 0, 1, 2, 3, 2, 1,
166                              0, -1, -2, -3, -2, -1,
167                              0, 1, 2, 3, 2, 1,
168                              0, -1, -2, -3, -2, -1,
169                              0, 1, 2, 3, 2, 1, 0 };
170         final List<WeightedObservedPoint> points = new ArrayList<WeightedObservedPoint>();
171         for (int i = 0; i < y.length; i++) {
172             points.add(new WeightedObservedPoint(1, i, y[i]));
173         }
174 
175         // The guesser fails because the function is far from an harmonic
176         // function: It is a triangular periodic function with amplitude 3
177         // and period 12, and all sample points are taken at integer abscissae
178         // so function values all belong to the integer subset {-3, -2, -1, 0,
179         // 1, 2, 3}.
180         new HarmonicCurveFitter.ParameterGuesser(points);
181     }
182 }
183