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