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    *      https://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.analysis.interpolation;
23  
24  import java.util.Random;
25  
26  import org.hipparchus.analysis.differentiation.DSFactory;
27  import org.hipparchus.analysis.differentiation.DerivativeStructure;
28  import org.hipparchus.analysis.differentiation.Gradient;
29  import org.hipparchus.analysis.polynomials.PolynomialFunction;
30  import org.hipparchus.exception.MathIllegalArgumentException;
31  import org.hipparchus.util.FastMath;
32  import org.junit.Assert;
33  import org.junit.Test;
34  
35  public class HermiteInterpolatorTest {
36  
37      @Test
38      public void testZero() {
39          HermiteInterpolator interpolator = new HermiteInterpolator();
40          interpolator.addSamplePoint(0.0, new double[] { 0.0 });
41          DSFactory factory = new DSFactory(1, 1);
42          for (double x = -10; x < 10; x += 1.0) {
43              DerivativeStructure y = interpolator.value(factory.variable(0, x))[0];
44              Assert.assertEquals(0.0, y.getValue(), 1.0e-15);
45              Assert.assertEquals(0.0, y.getPartialDerivative(1), 1.0e-15);
46          }
47          checkPolynomial(new PolynomialFunction(new double[] { 0.0 }),
48                          interpolator.getPolynomials()[0]);
49      }
50  
51      @Test
52      public void testQuadratic() {
53          HermiteInterpolator interpolator = new HermiteInterpolator();
54          interpolator.addSamplePoint(0.0, new double[] { 2.0 });
55          interpolator.addSamplePoint(1.0, new double[] { 0.0 });
56          interpolator.addSamplePoint(2.0, new double[] { 0.0 });
57          DSFactory factory = new DSFactory(1, 1);
58          for (double x = -10; x < 10; x += 1.0) {
59              DerivativeStructure y = interpolator.value(factory.variable(0, x))[0];
60              Assert.assertEquals((x - 1.0) * (x - 2.0), y.getValue(), 1.0e-15);
61              Assert.assertEquals(2 * x - 3.0, y.getPartialDerivative(1), 1.0e-15);
62          }
63          checkPolynomial(new PolynomialFunction(new double[] { 2.0, -3.0, 1.0 }),
64                          interpolator.getPolynomials()[0]);
65      }
66  
67      @Test
68      public void testMixedDerivatives() {
69          HermiteInterpolator interpolator = new HermiteInterpolator();
70          interpolator.addSamplePoint(0.0, new double[] { 1.0 }, new double[] { 2.0 });
71          interpolator.addSamplePoint(1.0, new double[] { 4.0 });
72          interpolator.addSamplePoint(2.0, new double[] { 5.0 }, new double[] { 2.0 });
73          DSFactory factory = new DSFactory(1, 1);
74          Assert.assertEquals(4, interpolator.getPolynomials()[0].degree());
75          DerivativeStructure y0 = interpolator.value(factory.variable(0, 0.0))[0];
76          Assert.assertEquals(1.0, y0.getValue(), 1.0e-15);
77          Assert.assertEquals(2.0, y0.getPartialDerivative(1), 1.0e-15);
78          double[][] d0 = interpolator.derivatives(0.0, 1);
79          Assert.assertEquals(1.0, d0[0][0], 1.0e-15);
80          Assert.assertEquals(2.0, d0[1][0], 1.0e-15);
81          Assert.assertEquals(4.0, interpolator.value(1.0)[0], 1.0e-15);
82          DerivativeStructure y2 = interpolator.value(factory.variable(0, 2.0))[0];
83          Assert.assertEquals(5.0, y2.getValue(), 1.0e-15);
84          Assert.assertEquals(2.0, y2.getPartialDerivative(1), 1.0e-15);
85          double[][] d2 = interpolator.derivatives(2.0, 1);
86          Assert.assertEquals(5.0, d2[0][0], 1.0e-15);
87          Assert.assertEquals(2.0, d2[1][0], 1.0e-15);
88          checkPolynomial(new PolynomialFunction(new double[] { 1.0, 2.0, 4.0, -4.0, 1.0 }),
89                          interpolator.getPolynomials()[0]);
90      }
91  
92      @Test
93      public void testRandomPolynomialsValuesOnly() {
94  
95          Random random = new Random(0x42b1e7dbd361a932l);
96  
97          for (int i = 0; i < 100; ++i) {
98  
99              int maxDegree = 0;
100             PolynomialFunction[] p = new PolynomialFunction[5];
101             for (int k = 0; k < p.length; ++k) {
102                 int degree = random.nextInt(7);
103                 p[k] = randomPolynomial(degree, random);
104                 maxDegree = FastMath.max(maxDegree, degree);
105             }
106 
107             HermiteInterpolator interpolator = new HermiteInterpolator();
108             for (int j = 0; j < 1 + maxDegree; ++j) {
109                 double x = 0.1 * j;
110                 double[] values = new double[p.length];
111                 for (int k = 0; k < p.length; ++k) {
112                     values[k] = p[k].value(x);
113                 }
114                 interpolator.addSamplePoint(x, values);
115             }
116 
117             for (double x = 0; x < 2; x += 0.1) {
118                 double[] values = interpolator.value(x);
119                 Assert.assertEquals(p.length, values.length);
120                 for (int k = 0; k < p.length; ++k) {
121                     Assert.assertEquals(p[k].value(x), values[k], 1.0e-8 * FastMath.abs(p[k].value(x)));
122                 }
123             }
124 
125             PolynomialFunction[] result = interpolator.getPolynomials();
126             for (int k = 0; k < p.length; ++k) {
127                 checkPolynomial(p[k], result[k]);
128             }
129 
130         }
131     }
132 
133     @Test
134     public void testRandomPolynomialsFirstDerivative() {
135 
136         Random random = new Random(0x570803c982ca5d3bl);
137 
138         for (int i = 0; i < 100; ++i) {
139 
140             int maxDegree = 0;
141             PolynomialFunction[] p      = new PolynomialFunction[5];
142             PolynomialFunction[] pPrime = new PolynomialFunction[5];
143             for (int k = 0; k < p.length; ++k) {
144                 int degree = random.nextInt(7);
145                 p[k]      = randomPolynomial(degree, random);
146                 pPrime[k] = p[k].polynomialDerivative();
147                 maxDegree = FastMath.max(maxDegree, degree);
148             }
149 
150             HermiteInterpolator interpolator = new HermiteInterpolator();
151             for (int j = 0; j < 1 + maxDegree / 2; ++j) {
152                 double x = 0.1 * j;
153                 double[] values      = new double[p.length];
154                 double[] derivatives = new double[p.length];
155                 for (int k = 0; k < p.length; ++k) {
156                     values[k]      = p[k].value(x);
157                     derivatives[k] = pPrime[k].value(x);
158                 }
159                 interpolator.addSamplePoint(x, values, derivatives);
160             }
161 
162             DSFactory factory = new DSFactory(1, 1);
163             for (double x = 0; x < 2; x += 0.1) {
164                 DerivativeStructure[] y = interpolator.value(factory.variable(0, x));
165                 Assert.assertEquals(p.length, y.length);
166                 for (int k = 0; k < p.length; ++k) {
167                     Assert.assertEquals(p[k].value(x), y[k].getValue(), 1.0e-8 * FastMath.abs(p[k].value(x)));
168                     Assert.assertEquals(pPrime[k].value(x), y[k].getPartialDerivative(1), 4.0e-8 * FastMath.abs(p[k].value(x)));
169                 }
170             }
171 
172             PolynomialFunction[] result = interpolator.getPolynomials();
173             for (int k = 0; k < p.length; ++k) {
174                 checkPolynomial(p[k], result[k]);
175             }
176 
177         }
178     }
179 
180     @Test
181     public void testSine() {
182         HermiteInterpolator interpolator = new HermiteInterpolator();
183         for (double x = 0; x < FastMath.PI; x += 0.5) {
184             interpolator.addSamplePoint(x, new double[] { FastMath.sin(x) });
185         }
186         DSFactory factory = new DSFactory(1, 2);
187         for (double x = 0.1; x <= 2.9; x += 0.01) {
188             DerivativeStructure y = interpolator.value(factory.variable(0, x))[0];
189             Assert.assertEquals( FastMath.sin(x), y.getValue(), 3.5e-5);
190             Assert.assertEquals( FastMath.cos(x), y.getPartialDerivative(1), 1.3e-4);
191             Assert.assertEquals(-FastMath.sin(x), y.getPartialDerivative(2), 2.9e-3);
192         }
193     }
194 
195     @Test
196     public void testSquareRoot() {
197         HermiteInterpolator interpolator = new HermiteInterpolator();
198         for (double x = 1.0; x < 3.6; x += 0.5) {
199             interpolator.addSamplePoint(x, new double[] { FastMath.sqrt(x) });
200         }
201         DSFactory factory = new DSFactory(1, 1);
202         for (double x = 1.1; x < 3.5; x += 0.01) {
203             DerivativeStructure y = interpolator.value(factory.variable(0, x))[0];
204             Assert.assertEquals(FastMath.sqrt(x), y.getValue(), 1.5e-4);
205             Assert.assertEquals(0.5 / FastMath.sqrt(x), y.getPartialDerivative(1), 8.5e-4);
206         }
207     }
208 
209     @Test
210     public void testWikipedia() {
211         // this test corresponds to the example from Wikipedia page:
212         // http://en.wikipedia.org/wiki/Hermite_interpolation
213         HermiteInterpolator interpolator = new HermiteInterpolator();
214         interpolator.addSamplePoint(-1, new double[] { 2 }, new double[] { -8 }, new double[] { 56 });
215         interpolator.addSamplePoint( 0, new double[] { 1 }, new double[] {  0 }, new double[] {  0 });
216         interpolator.addSamplePoint( 1, new double[] { 2 }, new double[] {  8 }, new double[] { 56 });
217         DSFactory factory = new DSFactory(1, 1);
218         for (double x = -1.0; x <= 1.0; x += 0.125) {
219             DerivativeStructure y = interpolator.value(factory.variable(0, x))[0];
220             double x2 = x * x;
221             double x4 = x2 * x2;
222             double x8 = x4 * x4;
223             Assert.assertEquals(x8 + 1, y.getValue(), 1.0e-15);
224             Assert.assertEquals(8 * x4 * x2 * x, y.getPartialDerivative(1), 1.0e-15);
225         }
226         checkPolynomial(new PolynomialFunction(new double[] { 1, 0, 0, 0, 0, 0, 0, 0, 1 }),
227                         interpolator.getPolynomials()[0]);
228     }
229 
230     @Test
231     public void testWikipediaGradient() {
232         // this test corresponds to the example from Wikipedia page:
233         // http://en.wikipedia.org/wiki/Hermite_interpolation
234         HermiteInterpolator interpolator = new HermiteInterpolator();
235         interpolator.addSamplePoint(-1, new double[] { 2 }, new double[] { -8 }, new double[] { 56 });
236         interpolator.addSamplePoint( 0, new double[] { 1 }, new double[] {  0 }, new double[] {  0 });
237         interpolator.addSamplePoint( 1, new double[] { 2 }, new double[] {  8 }, new double[] { 56 });
238         for (double x = -1.0; x <= 1.0; x += 0.125) {
239             Gradient y = interpolator.value(Gradient.variable(1, 0, x))[0];
240             double x2 = x * x;
241             double x4 = x2 * x2;
242             double x8 = x4 * x4;
243             Assert.assertEquals(x8 + 1, y.getValue(), 1.0e-15);
244             Assert.assertEquals(8 * x4 * x2 * x, y.getPartialDerivative(0), 1.0e-15);
245         }
246         checkPolynomial(new PolynomialFunction(new double[] { 1, 0, 0, 0, 0, 0, 0, 0, 1 }),
247                         interpolator.getPolynomials()[0]);
248     }
249 
250     @Test
251     public void testOnePointParabola() {
252         HermiteInterpolator interpolator = new HermiteInterpolator();
253         interpolator.addSamplePoint(0, new double[] { 1 }, new double[] { 1 }, new double[] { 2 });
254         DSFactory factory = new DSFactory(1, 1);
255         for (double x = -1.0; x <= 1.0; x += 0.125) {
256             DerivativeStructure y = interpolator.value(factory.variable(0, x))[0];
257             Assert.assertEquals(1 + x * (1 + x), y.getValue(), 1.0e-15);
258             Assert.assertEquals(1 + 2 * x, y.getPartialDerivative(1), 1.0e-15);
259         }
260         checkPolynomial(new PolynomialFunction(new double[] { 1, 1, 1 }),
261                         interpolator.getPolynomials()[0]);
262     }
263 
264     private PolynomialFunction randomPolynomial(int degree, Random random) {
265         double[] coeff = new double[ 1 + degree];
266         for (int j = 0; j < degree; ++j) {
267             coeff[j] = random.nextDouble();
268         }
269         return new PolynomialFunction(coeff);
270     }
271 
272     @Test(expected=MathIllegalArgumentException.class)
273     public void testEmptySample() {
274         new HermiteInterpolator().value(0.0);
275     }
276 
277     @Test(expected=MathIllegalArgumentException.class)
278     public void testDuplicatedAbscissa() {
279         HermiteInterpolator interpolator = new HermiteInterpolator();
280         interpolator.addSamplePoint(1.0, new double[] { 0.0 });
281         interpolator.addSamplePoint(1.0, new double[] { 1.0 });
282     }
283 
284     private void checkPolynomial(PolynomialFunction expected, PolynomialFunction result) {
285         Assert.assertTrue(result.degree() >= expected.degree());
286         double[] cE = expected.getCoefficients();
287         double[] cR = result.getCoefficients();
288         for (int i = 0; i < cE.length; ++i) {
289             Assert.assertEquals(cE[i], cR[i], 1.0e-8 * FastMath.abs(cE[i]));
290         }
291         for (int i = cE.length; i < cR.length; ++i) {
292             Assert.assertEquals(0.0, cR[i], 1.0e-9);
293         }
294     }
295 
296 }
297