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  
23  package org.hipparchus.analysis.differentiation;
24  
25  import org.hipparchus.UnitTestUtils;
26  import org.hipparchus.analysis.QuinticFunction;
27  import org.hipparchus.analysis.UnivariateFunction;
28  import org.hipparchus.analysis.UnivariateMatrixFunction;
29  import org.hipparchus.analysis.UnivariateVectorFunction;
30  import org.hipparchus.analysis.function.Gaussian;
31  import org.hipparchus.analysis.function.Sin;
32  import org.hipparchus.exception.LocalizedCoreFormats;
33  import org.hipparchus.exception.MathIllegalArgumentException;
34  import org.hipparchus.exception.MathRuntimeException;
35  import org.hipparchus.util.FastMath;
36  import org.junit.Assert;
37  import org.junit.Test;
38  
39  /**
40   * Test for class {@link FiniteDifferencesDifferentiator}.
41   */
42  public class FiniteDifferencesDifferentiatorTest {
43  
44      @Test(expected=MathIllegalArgumentException.class)
45      public void testWrongNumberOfPoints() {
46          new FiniteDifferencesDifferentiator(1, 1.0);
47      }
48  
49      @Test(expected=MathIllegalArgumentException.class)
50      public void testWrongStepSize() {
51          new FiniteDifferencesDifferentiator(3, 0.0);
52      }
53  
54      @Test
55      public void testSerialization() {
56          FiniteDifferencesDifferentiator differentiator =
57                  new FiniteDifferencesDifferentiator(3, 1.0e-3);
58          FiniteDifferencesDifferentiator recovered =
59                  (FiniteDifferencesDifferentiator) UnitTestUtils.serializeAndRecover(differentiator);
60          Assert.assertEquals(differentiator.getNbPoints(), recovered.getNbPoints());
61          Assert.assertEquals(differentiator.getStepSize(), recovered.getStepSize(), 1.0e-15);
62      }
63  
64      @Test
65      public void testConstant() {
66          FiniteDifferencesDifferentiator differentiator =
67                  new FiniteDifferencesDifferentiator(5, 0.01);
68          UnivariateDifferentiableFunction f =
69                  differentiator.differentiate(new UnivariateFunction() {
70                      @Override
71                      public double value(double x) {
72                          return 42.0;
73                      }
74                  });
75          DSFactory factory = new DSFactory(1, 2);
76          for (double x = -10; x < 10; x += 0.1) {
77              DerivativeStructure y = f.value(factory.variable(0, x));
78              Assert.assertEquals(42.0, y.getValue(), 1.0e-15);
79              Assert.assertEquals( 0.0, y.getPartialDerivative(1), 1.0e-15);
80              Assert.assertEquals( 0.0, y.getPartialDerivative(2), 1.0e-15);
81          }
82      }
83  
84      @Test
85      public void testLinear() {
86          FiniteDifferencesDifferentiator differentiator =
87                  new FiniteDifferencesDifferentiator(5, 0.01);
88          UnivariateDifferentiableFunction f =
89                  differentiator.differentiate(new UnivariateFunction() {
90                      @Override
91                      public double value(double x) {
92                          return 2 - 3 * x;
93                      }
94                  });
95          DSFactory factory = new DSFactory(1, 2);
96          for (double x = -10; x < 10; x += 0.1) {
97              DerivativeStructure y = f.value(factory.variable(0, x));
98              Assert.assertEquals("" + (2 - 3 * x - y.getValue()), 2 - 3 * x, y.getValue(), 2.0e-15);
99              Assert.assertEquals(-3.0, y.getPartialDerivative(1), 4.0e-13);
100             Assert.assertEquals( 0.0, y.getPartialDerivative(2), 9.0e-11);
101         }
102     }
103 
104     @Test
105     public void testGaussian() {
106         FiniteDifferencesDifferentiator differentiator =
107                 new FiniteDifferencesDifferentiator(9, 0.02);
108         UnivariateDifferentiableFunction gaussian = new Gaussian(1.0, 2.0);
109         UnivariateDifferentiableFunction f =
110                 differentiator.differentiate(gaussian);
111         double[] expectedError = new double[] {
112             6.939e-18, 1.284e-15, 2.477e-13, 1.168e-11, 2.840e-9, 7.971e-8
113         };
114         double[] maxError = new double[expectedError.length];
115         DSFactory factory = new DSFactory(1, maxError.length - 1);
116         for (double x = -10; x < 10; x += 0.1) {
117             DerivativeStructure dsX  = factory.variable(0, x);
118             DerivativeStructure yRef = gaussian.value(dsX);
119             DerivativeStructure y    = f.value(dsX);
120             Assert.assertEquals(f.value(dsX.getValue()), f.value(dsX).getValue(), 1.0e-15);
121             for (int order = 0; order <= yRef.getOrder(); ++order) {
122                 maxError[order] = FastMath.max(maxError[order],
123                                         FastMath.abs(yRef.getPartialDerivative(order) -
124                                                      y.getPartialDerivative(order)));
125             }
126         }
127         for (int i = 0; i < maxError.length; ++i) {
128             Assert.assertEquals(expectedError[i], maxError[i], 0.01 * expectedError[i]);
129         }
130     }
131 
132     @Test
133     public void testStepSizeUnstability() {
134         UnivariateDifferentiableFunction quintic = new QuinticFunction();
135         UnivariateDifferentiableFunction goodStep =
136                 new FiniteDifferencesDifferentiator(7, 0.25).differentiate(quintic);
137         UnivariateDifferentiableFunction badStep =
138                 new FiniteDifferencesDifferentiator(7, 1.0e-6).differentiate(quintic);
139         double[] maxErrorGood = new double[7];
140         double[] maxErrorBad  = new double[7];
141         DSFactory factory = new DSFactory(1, maxErrorGood.length - 1);
142         for (double x = -10; x < 10; x += 0.1) {
143             DerivativeStructure dsX  = factory.variable(0, x);
144             DerivativeStructure yRef  = quintic.value(dsX);
145             DerivativeStructure yGood = goodStep.value(dsX);
146             DerivativeStructure yBad  = badStep.value(dsX);
147             for (int order = 0; order <= 6; ++order) {
148                 maxErrorGood[order] = FastMath.max(maxErrorGood[order],
149                                                    FastMath.abs(yRef.getPartialDerivative(order) -
150                                                                 yGood.getPartialDerivative(order)));
151                 maxErrorBad[order]  = FastMath.max(maxErrorBad[order],
152                                                    FastMath.abs(yRef.getPartialDerivative(order) -
153                                                                 yBad.getPartialDerivative(order)));
154             }
155         }
156 
157         // the 0.25 step size is good for finite differences in the quintic on this abscissa range for 7 points
158         // the errors are fair
159         final double[] expectedGood = new double[] {
160             7.276e-12, 7.276e-11, 9.968e-10, 3.092e-9, 5.432e-8, 8.196e-8, 1.818e-6
161         };
162 
163         // the 1.0e-6 step size is far too small for finite differences in the quintic on this abscissa range for 7 points
164         // the errors are huge!
165         final double[] expectedBad = new double[] {
166             2.910e-11, 2.087e-5, 147.7, 3.820e7, 6.354e14, 6.548e19, 1.543e27
167         };
168 
169         for (int i = 0; i < maxErrorGood.length; ++i) {
170             Assert.assertEquals(expectedGood[i], maxErrorGood[i], 0.01 * expectedGood[i]);
171             Assert.assertEquals(expectedBad[i],  maxErrorBad[i],  0.01 * expectedBad[i]);
172         }
173 
174     }
175 
176     @Test(expected=MathIllegalArgumentException.class)
177     public void testWrongOrder() {
178         UnivariateDifferentiableFunction f =
179                 new FiniteDifferencesDifferentiator(3, 0.01).differentiate(new UnivariateFunction() {
180                     @Override
181                     public double value(double x) {
182                         // this exception should not be thrown because wrong order
183                         // should be detected before function call
184                         throw MathRuntimeException.createInternalError();
185                     }
186                 });
187         f.value(new DSFactory(1, 3).variable(0, 1.0));
188     }
189 
190     @Test(expected=MathIllegalArgumentException.class)
191     public void testWrongOrderVector() {
192         UnivariateDifferentiableVectorFunction f =
193                 new FiniteDifferencesDifferentiator(3, 0.01).differentiate(new UnivariateVectorFunction() {
194                     @Override
195                     public double[] value(double x) {
196                         // this exception should not be thrown because wrong order
197                         // should be detected before function call
198                         throw MathRuntimeException.createInternalError();
199                     }
200                 });
201         f.value(new DSFactory(1, 3).variable(0, 1.0));
202     }
203 
204     @Test(expected=MathIllegalArgumentException.class)
205     public void testWrongOrderMatrix() {
206         UnivariateDifferentiableMatrixFunction f =
207                 new FiniteDifferencesDifferentiator(3, 0.01).differentiate(new UnivariateMatrixFunction() {
208                     @Override
209                     public double[][] value(double x) {
210                         // this exception should not be thrown because wrong order
211                         // should be detected before function call
212                         throw MathRuntimeException.createInternalError();
213                     }
214                 });
215         f.value(new DSFactory(1, 3).variable(0, 1.0));
216     }
217 
218     @Test(expected=MathIllegalArgumentException.class)
219     public void testTooLargeStep() {
220         new FiniteDifferencesDifferentiator(3, 2.5, 0.0, 1.0);
221     }
222 
223     @Test
224     public void testBounds() {
225 
226         final double slope = 2.5;
227         UnivariateFunction f = new UnivariateFunction() {
228             @Override
229             public double value(double x) {
230                 if (x < 0) {
231                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_SMALL,
232                                                            x, 0);
233                 } else if (x > 1) {
234                     throw new MathIllegalArgumentException(LocalizedCoreFormats.NUMBER_TOO_LARGE,
235                                                            x, 1);
236                 } else {
237                     return slope * x;
238                 }
239             }
240         };
241 
242         UnivariateDifferentiableFunction missingBounds =
243                 new FiniteDifferencesDifferentiator(3, 0.1).differentiate(f);
244         UnivariateDifferentiableFunction properlyBounded =
245                 new FiniteDifferencesDifferentiator(3, 0.1, 0.0, 1.0).differentiate(f);
246         DSFactory factory = new DSFactory(1, 1);
247         DerivativeStructure tLow  = factory.variable(0, 0.05);
248         DerivativeStructure tHigh = factory.variable(0, 0.95);
249 
250         try {
251             // here, we did not set the bounds, so the differences are evaluated out of domain
252             // using f(-0.05), f(0.05), f(0.15)
253             missingBounds.value(tLow);
254             Assert.fail("an exception should have been thrown");
255         } catch (MathIllegalArgumentException nse) {
256             Assert.assertEquals(LocalizedCoreFormats.NUMBER_TOO_SMALL, nse.getSpecifier());
257             Assert.assertEquals(-0.05, ((Double) nse.getParts()[0]).doubleValue(), 1.0e-10);
258         } catch (Exception e) {
259             Assert.fail("wrong exception caught: " + e.getClass().getName());
260         }
261 
262         try {
263             // here, we did not set the bounds, so the differences are evaluated out of domain
264             // using f(0.85), f(0.95), f(1.05)
265             missingBounds.value(tHigh);
266             Assert.fail("an exception should have been thrown");
267         } catch (MathIllegalArgumentException nle) {
268             Assert.assertEquals(LocalizedCoreFormats.NUMBER_TOO_LARGE, nle.getSpecifier());
269             Assert.assertEquals(1.05, ((Double) nle.getParts()[0]).doubleValue(), 1.0e-10);
270         } catch (Exception e) {
271             Assert.fail("wrong exception caught: " + e.getClass().getName());
272         }
273 
274         // here, we did set the bounds, so evaluations are done within domain
275         // using f(0.0), f(0.1), f(0.2)
276         Assert.assertEquals(slope, properlyBounded.value(tLow).getPartialDerivative(1), 1.0e-10);
277 
278         // here, we did set the bounds, so evaluations are done within domain
279         // using f(0.8), f(0.9), f(1.0)
280         Assert.assertEquals(slope, properlyBounded.value(tHigh).getPartialDerivative(1), 1.0e-10);
281 
282     }
283 
284     @Test
285     public void testBoundedSqrt() {
286 
287         UnivariateFunctionDifferentiator differentiator =
288                 new FiniteDifferencesDifferentiator(9, 1.0 / 32, 0.0, Double.POSITIVE_INFINITY);
289         UnivariateDifferentiableFunction sqrt = differentiator.differentiate(new UnivariateFunction() {
290             @Override
291             public double value(double x) {
292                 return FastMath.sqrt(x);
293             }
294         });
295 
296         // we are able to compute derivative near 0, but the accuracy is much poorer there
297         DSFactory factory = new DSFactory(1, 1);
298         DerivativeStructure t001 = factory.variable(0, 0.01);
299         Assert.assertEquals(0.5 / FastMath.sqrt(t001.getValue()), sqrt.value(t001).getPartialDerivative(1), 1.6);
300         DerivativeStructure t01 = factory.variable(0, 0.1);
301         Assert.assertEquals(0.5 / FastMath.sqrt(t01.getValue()), sqrt.value(t01).getPartialDerivative(1), 7.0e-3);
302         DerivativeStructure t03 = factory.variable(0, 0.3);
303         Assert.assertEquals(0.5 / FastMath.sqrt(t03.getValue()), sqrt.value(t03).getPartialDerivative(1), 2.1e-7);
304 
305     }
306 
307     @Test
308     public void testVectorFunction() {
309 
310         FiniteDifferencesDifferentiator differentiator =
311                 new FiniteDifferencesDifferentiator(7, 0.01);
312         UnivariateDifferentiableVectorFunction f =
313                 differentiator.differentiate(new UnivariateVectorFunction() {
314 
315             @Override
316             public double[] value(double x) {
317                 return new double[] { FastMath.cos(x), FastMath.sin(x) };
318             }
319 
320         });
321 
322         DSFactory factory = new DSFactory(1, 2);
323         for (double x = -10; x < 10; x += 0.1) {
324             DerivativeStructure dsX = factory.variable(0, x);
325             DerivativeStructure[] y = f.value(dsX);
326             double cos = FastMath.cos(x);
327             double sin = FastMath.sin(x);
328             double[] f1 = f.value(dsX.getValue());
329             DerivativeStructure[] f2 = f.value(dsX);
330             Assert.assertEquals(f1.length, f2.length);
331             for (int i = 0; i < f1.length; ++i) {
332                 Assert.assertEquals(f1[i], f2[i].getValue(), 1.0e-15);
333             }
334             Assert.assertEquals( cos, y[0].getValue(), 7.0e-16);
335             Assert.assertEquals( sin, y[1].getValue(), 7.0e-16);
336             Assert.assertEquals(-sin, y[0].getPartialDerivative(1), 6.0e-14);
337             Assert.assertEquals( cos, y[1].getPartialDerivative(1), 6.0e-14);
338             Assert.assertEquals(-cos, y[0].getPartialDerivative(2), 2.0e-11);
339             Assert.assertEquals(-sin, y[1].getPartialDerivative(2), 2.0e-11);
340         }
341 
342     }
343 
344     @Test
345     public void testMatrixFunction() {
346 
347         FiniteDifferencesDifferentiator differentiator =
348                 new FiniteDifferencesDifferentiator(7, 0.01);
349         UnivariateDifferentiableMatrixFunction f =
350                 differentiator.differentiate(new UnivariateMatrixFunction() {
351 
352             @Override
353             public double[][] value(double x) {
354                 return new double[][] {
355                     { FastMath.cos(x),  FastMath.sin(x)  },
356                     { FastMath.cosh(x), FastMath.sinh(x) }
357                 };
358             }
359 
360         });
361 
362         DSFactory factory = new DSFactory(1, 2);
363         for (double x = -1; x < 1; x += 0.02) {
364             DerivativeStructure dsX = factory.variable(0, x);
365             DerivativeStructure[][] y = f.value(dsX);
366             double cos = FastMath.cos(x);
367             double sin = FastMath.sin(x);
368             double cosh = FastMath.cosh(x);
369             double sinh = FastMath.sinh(x);
370             double[][] f1 = f.value(dsX.getValue());
371             DerivativeStructure[][] f2 = f.value(dsX);
372             Assert.assertEquals(f1.length, f2.length);
373             for (int i = 0; i < f1.length; ++i) {
374                 Assert.assertEquals(f1[i].length, f2[i].length);
375                 for (int j = 0; j < f1[i].length; ++j) {
376                     Assert.assertEquals(f1[i][j], f2[i][j].getValue(), 1.0e-15);
377                 }
378             }
379             Assert.assertEquals(cos,   y[0][0].getValue(), 7.0e-18);
380             Assert.assertEquals(sin,   y[0][1].getValue(), 6.0e-17);
381             Assert.assertEquals(cosh,  y[1][0].getValue(), 3.0e-16);
382             Assert.assertEquals(sinh,  y[1][1].getValue(), 3.0e-16);
383             Assert.assertEquals(-sin,  y[0][0].getPartialDerivative(1), 2.0e-14);
384             Assert.assertEquals( cos,  y[0][1].getPartialDerivative(1), 2.0e-14);
385             Assert.assertEquals( sinh, y[1][0].getPartialDerivative(1), 3.0e-14);
386             Assert.assertEquals( cosh, y[1][1].getPartialDerivative(1), 3.0e-14);
387             Assert.assertEquals(-cos,  y[0][0].getPartialDerivative(2), 3.0e-12);
388             Assert.assertEquals(-sin,  y[0][1].getPartialDerivative(2), 3.0e-12);
389             Assert.assertEquals( cosh, y[1][0].getPartialDerivative(2), 6.0e-12);
390             Assert.assertEquals( sinh, y[1][1].getPartialDerivative(2), 6.0e-12);
391         }
392 
393     }
394 
395     @Test
396     public void testSeveralFreeParameters() {
397         FiniteDifferencesDifferentiator differentiator =
398                 new FiniteDifferencesDifferentiator(5, 0.001);
399         UnivariateDifferentiableFunction sine = new Sin();
400         UnivariateDifferentiableFunction f =
401                 differentiator.differentiate(sine);
402         double[] expectedError = new double[] {
403             6.696e-16, 1.371e-12, 2.007e-8, 1.754e-5
404         };
405         double[] maxError = new double[expectedError.length];
406         DSFactory factory = new DSFactory(2, maxError.length - 1);
407         for (double x = -2; x < 2; x += 0.1) {
408            for (double y = -2; y < 2; y += 0.1) {
409                DerivativeStructure dsX  = factory.variable(0, x);
410                DerivativeStructure dsY  = factory.variable(1, y);
411                DerivativeStructure dsT  = dsX.multiply(3).subtract(dsY.multiply(2));
412                DerivativeStructure sRef = sine.value(dsT);
413                DerivativeStructure s    = f.value(dsT);
414                for (int xOrder = 0; xOrder <= sRef.getOrder(); ++xOrder) {
415                    for (int yOrder = 0; yOrder <= sRef.getOrder(); ++yOrder) {
416                        if (xOrder + yOrder <= sRef.getOrder()) {
417                            maxError[xOrder +yOrder] = FastMath.max(maxError[xOrder + yOrder],
418                                                                     FastMath.abs(sRef.getPartialDerivative(xOrder, yOrder) -
419                                                                                  s.getPartialDerivative(xOrder, yOrder)));
420                        }
421                    }
422                }
423            }
424        }
425        for (int i = 0; i < maxError.length; ++i) {
426            Assert.assertEquals(expectedError[i], maxError[i], 0.01 * expectedError[i]);
427        }
428     }
429 
430 }