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.nabla.differences;
18  
19  import java.util.Random;
20  
21  import org.apache.commons.nabla.Polynomial;
22  import org.apache.commons.nabla.ReferenceFunction;
23  import org.apache.commons.nabla.core.DifferentialPair;
24  import org.apache.commons.nabla.core.DifferentiationException;
25  import org.apache.commons.nabla.core.UnivariateDerivative;
26  import org.apache.commons.nabla.differences.FiniteDifferencesDifferentiator;
27  
28  import junit.framework.TestCase;
29  
30  public abstract class AbstractFiniteDifferencesTest extends TestCase {
31  
32      protected abstract FiniteDifferencesDifferentiator buildDifferentiator(double h);
33      protected abstract int getOrder();
34  
35      /** Check the reference implementation for polynomials. */
36      public void testReference() {
37  
38          Polynomial legendre5  = new Polynomial(new double[] { 63, 0, -70, 0, 15, 0 });
39          assertEquals(17155.92466365559104, legendre5.f(Math.PI), 1.0e-9);
40          assertEquals(28626.24675148200242, legendre5.fPrime(Math.PI), 1.0e-8);
41  
42          Polynomial legendre4  = new Polynomial(new double[] { 35, 0, -30, 0, 3 });
43          assertEquals(3116.230054157404545, legendre4.f(Math.PI), 1.0e-10);
44          assertEquals(4152.383176026587230, legendre4.fPrime(Math.PI), 1.0e-9);
45  
46          Polynomial chebyshev7 = new Polynomial(new double[] { 64, 0, -112, 0, 56, 0, -7, 0 });
47          assertEquals(160738.9222272848309, chebyshev7.f(Math.PI), 1.0e-8);
48          assertEquals(377804.3612820780352, chebyshev7.fPrime(Math.PI), 1.0e-7);
49  
50          Polynomial laguerre3  = new Polynomial(new double[] { -1, 9, -18, 6 });
51          assertEquals(7.271495164888129102, laguerre3.f(Math.PI), 1.0e-13);
52          assertEquals(8.939854561348202436, laguerre3.fPrime(Math.PI), 1.0e-12);
53  
54      }
55  
56      /** The differentiation schemes of order p should compute exact differentials
57       * for polynomials up to degree p.
58       */
59      public void testExactPolynomialDifferentiation() throws DifferentiationException {
60          for (int k = 0; k <= getOrder(); ++k) {
61              ReferenceFunction reference = Polynomial.randomPolynomial(random, k);
62              UnivariateDerivative derivative = buildDifferentiator(0.1).differentiate(reference);
63              for (DifferentialPair t : transcendentals) {
64                  double error = derivative.f(t).getFirstDerivative() - reference.fPrime(t.getValue());
65                  assertEquals(0, error, 1.0e-13 * Math.abs(reference.fPrime(t.getValue())));
66              }
67          }
68      }
69  
70      /** The derivative instances are bound to their primitive
71       * and should follow their updates transparently. */
72      public void testPrimitiveBinding() throws DifferentiationException {
73  
74          // compute the differential once only
75          Polynomial reference = Polynomial.randomPolynomial(random, 3);
76          FiniteDifferencesDifferentiator differentiator = buildDifferentiator(0.1);
77          UnivariateDerivative derivative = differentiator.differentiate(reference);
78          double errorScaleFactor = differentiator.getSignedErrorScaleFactor();
79  
80  
81          // compute the differential values from this initial state
82          DifferentialPair[] beforeChange =
83              new DifferentialPair[transcendentals.length];
84          for (int i = 0; i < transcendentals.length; ++i) {
85              DifferentialPair t = transcendentals[i];
86              beforeChange[i] = derivative.f(t);
87              assertEquals(0, beforeChange[i].getFirstDerivative() - reference.fPrime(t.getValue()), 100 * Math.abs(errorScaleFactor));
88          }
89  
90          // change the primitive WITHOUT recomputing the differential
91          ((Polynomial) derivative.getPrimitive()).addToSelf(Polynomial.randomPolynomial(random, 4));
92  
93          // compute the new values reusing the original differential
94          DifferentialPair[] afterChange =
95              new DifferentialPair[transcendentals.length];
96          for (int i = 0; i < transcendentals.length; ++i) {
97              DifferentialPair t = transcendentals[i];
98              afterChange[i] = derivative.f(t);
99              assertEquals(0, afterChange[i].getFirstDerivative() - reference.fPrime(t.getValue()), 100 * Math.abs(errorScaleFactor));
100         }
101 
102         // check the change was important
103         for (int i = 0; i < transcendentals.length; ++i) {
104             DifferentialPair change =
105                 DifferentialPair.subtract(beforeChange[i], afterChange[i]);
106             assertTrue(Math.abs(change.getValue()) > 10000.0 * Math.abs(errorScaleFactor));
107             assertTrue(Math.abs(change.getFirstDerivative()) > 10000.0 * Math.abs(errorScaleFactor));
108         }
109 
110     }
111 
112     /** Check the differentials values in a range.
113      * @param reference reference differential
114      * @param h finite differences step
115      * @param lower lower bound of the test range
116      * @param upper upper bound of the test range
117      * @param n number of test points in the range
118      * @param epsilon tolerance on the differentials values
119      */
120     protected void checkValues(ReferenceFunction reference, double h,
121                                double lower, double upper, int n,
122                                double epsilon) {
123         try {
124             UnivariateDerivative derivative = buildDifferentiator(h).differentiate(reference);
125             int n1 = n - 1;
126             for (int i = 0; i < n; ++i) {
127                 double t = ((n1 - i) * lower + i * upper) / n1;
128                 assertEquals(reference.fPrime(t),
129                              derivative.f(DifferentialPair.newVariable(t)).getFirstDerivative(),
130                              epsilon);
131             }
132         } catch (DifferentiationException de) {
133             fail(de.getMessage());
134         }
135     }
136 
137     /** Check the error order at some point.
138      * @param reference reference differential
139      * @param t abscissa of the test point
140      * @param hMin minimal step
141      * @param hMax minimal step
142      * @param n number of steps to try
143      * @param error expected order-normalized error
144      * @param epsilon tolerance on the expected extremal values
145      */
146     protected void checkOrder(ReferenceFunction reference,
147                               double t, double hMin, double hMax, int n,
148                               double error, double epsilon) {
149         try {
150             int n1 = n - 1;
151             for (int i = 0; i < n; ++i) {
152                 double h = ((n1 - i) * hMin + i * hMax) / n1;
153                 FiniteDifferencesDifferentiator fds = buildDifferentiator(h);
154                 UnivariateDerivative derivative = fds.differentiate(reference);
155                 double rawError =
156                     derivative.f(DifferentialPair.newVariable(t)).getFirstDerivative() - reference.fPrime(t);
157                 double orderNormalizedError = Math.abs(rawError / fds.getSignedErrorScaleFactor());
158                 assertEquals(error, orderNormalizedError, epsilon);
159             }
160         } catch (DifferentiationException de) {
161             fail(de.getMessage());
162         }
163     }
164 
165     protected ReferenceFunction hugeDifferentialsFunction() {
166         return new ReferenceFunction() {
167             public double f(double t) {
168                 return Math.sin(3 * Math.tan(t * t) - Math.log(t));
169             }
170             public double fPrime(double t) {
171                 double t2 = t * t;
172                 double cosT2 = Math.cos(t2);
173                 double tanT2 = Math.tan(t2);
174                 return (6 * t / (cosT2 * cosT2) - 1 / t) * Math.cos(3 * tanT2 - Math.log(t));
175             }
176         };
177     }
178 
179     public void setUp() {
180         random = new Random(0x6c05e9f92868d664l);
181         transcendentals = new DifferentialPair[] {
182                                                   DifferentialPair.newVariable(Math.PI),
183                                                   DifferentialPair.newVariable(Math.E),
184                                                   DifferentialPair.newVariable(Math.pow(2, Math.sqrt(2)))
185         };
186 
187     }
188 
189     public void tearDown() {
190         random = null;
191         transcendentals = null;
192     }
193 
194     protected Random random;
195     private DifferentialPair[] transcendentals;
196 
197 }