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 java.lang.reflect.Field;
26  import java.lang.reflect.InvocationTargetException;
27  import java.lang.reflect.Method;
28  import java.util.HashMap;
29  import java.util.Map;
30  import java.util.stream.Stream;
31  
32  import org.hipparchus.exception.MathIllegalArgumentException;
33  import org.hipparchus.util.CombinatoricsUtils;
34  import org.junit.Assert;
35  import org.junit.Test;
36  
37  
38  /**
39   * Test for class {@link DSCompiler}.
40   */
41  public class DSCompilerTest {
42  
43      @Test
44      public void testSize() {
45          for (int i = 0; i < 6; ++i) {
46              for (int j = 0; j < 6; ++j) {
47                  long expected = CombinatoricsUtils.binomialCoefficient(i + j, i);
48                  Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
49                  Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
50              }
51          }
52      }
53  
54      @Test
55      public void testIndices() {
56  
57          DSCompiler c = DSCompiler.getCompiler(0, 0);
58          checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
59  
60          c = DSCompiler.getCompiler(0, 1);
61          checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
62  
63          c = DSCompiler.getCompiler(1, 0);
64          checkIndices(c.getPartialDerivativeOrders(0), 0);
65  
66          c = DSCompiler.getCompiler(1, 1);
67          checkIndices(c.getPartialDerivativeOrders(0), 0);
68          checkIndices(c.getPartialDerivativeOrders(1), 1);
69  
70          c = DSCompiler.getCompiler(1, 2);
71          checkIndices(c.getPartialDerivativeOrders(0), 0);
72          checkIndices(c.getPartialDerivativeOrders(1), 1);
73          checkIndices(c.getPartialDerivativeOrders(2), 2);
74  
75          c = DSCompiler.getCompiler(2, 1);
76          checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
77          checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
78          checkIndices(c.getPartialDerivativeOrders(2), 0, 1);
79  
80          c = DSCompiler.getCompiler(1, 3);
81          checkIndices(c.getPartialDerivativeOrders(0), 0);
82          checkIndices(c.getPartialDerivativeOrders(1), 1);
83          checkIndices(c.getPartialDerivativeOrders(2), 2);
84          checkIndices(c.getPartialDerivativeOrders(3), 3);
85  
86          c = DSCompiler.getCompiler(2, 2);
87          checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
88          checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
89          checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
90          checkIndices(c.getPartialDerivativeOrders(3), 0, 1);
91          checkIndices(c.getPartialDerivativeOrders(4), 1, 1);
92          checkIndices(c.getPartialDerivativeOrders(5), 0, 2);
93  
94          c = DSCompiler.getCompiler(3, 1);
95          checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
96          checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
97          checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0);
98          checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1);
99  
100         c = DSCompiler.getCompiler(1, 4);
101         checkIndices(c.getPartialDerivativeOrders(0), 0);
102         checkIndices(c.getPartialDerivativeOrders(1), 1);
103         checkIndices(c.getPartialDerivativeOrders(2), 2);
104         checkIndices(c.getPartialDerivativeOrders(3), 3);
105         checkIndices(c.getPartialDerivativeOrders(4), 4);
106 
107         c = DSCompiler.getCompiler(2, 3);
108         checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
109         checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
110         checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
111         checkIndices(c.getPartialDerivativeOrders(3), 3, 0);
112         checkIndices(c.getPartialDerivativeOrders(4), 0, 1);
113         checkIndices(c.getPartialDerivativeOrders(5), 1, 1);
114         checkIndices(c.getPartialDerivativeOrders(6), 2, 1);
115         checkIndices(c.getPartialDerivativeOrders(7), 0, 2);
116         checkIndices(c.getPartialDerivativeOrders(8), 1, 2);
117         checkIndices(c.getPartialDerivativeOrders(9), 0, 3);
118 
119         c = DSCompiler.getCompiler(3, 2);
120         checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
121         checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
122         checkIndices(c.getPartialDerivativeOrders(2), 2, 0, 0);
123         checkIndices(c.getPartialDerivativeOrders(3), 0, 1, 0);
124         checkIndices(c.getPartialDerivativeOrders(4), 1, 1, 0);
125         checkIndices(c.getPartialDerivativeOrders(5), 0, 2, 0);
126         checkIndices(c.getPartialDerivativeOrders(6), 0, 0, 1);
127         checkIndices(c.getPartialDerivativeOrders(7), 1, 0, 1);
128         checkIndices(c.getPartialDerivativeOrders(8), 0, 1, 1);
129         checkIndices(c.getPartialDerivativeOrders(9), 0, 0, 2);
130 
131         c = DSCompiler.getCompiler(4, 1);
132         checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0, 0);
133         checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0, 0);
134         checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0, 0);
135         checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1, 0);
136         checkIndices(c.getPartialDerivativeOrders(4), 0, 0, 0, 1);
137 
138     }
139 
140     @Test(expected=MathIllegalArgumentException.class)
141     public void testIncompatibleParams() {
142         DSCompiler.getCompiler(3, 2).checkCompatibility(DSCompiler.getCompiler(4, 2));
143     }
144 
145     @Test(expected=MathIllegalArgumentException.class)
146     public void testIncompatibleOrder() {
147         DSCompiler.getCompiler(3, 3).checkCompatibility(DSCompiler.getCompiler(3, 2));
148     }
149 
150     @Test
151     public void testSymmetry() {
152         for (int i = 0; i < 6; ++i) {
153             for (int j = 0; j < 6; ++j) {
154                 DSCompiler c = DSCompiler.getCompiler(i, j);
155                 for (int k = 0; k < c.getSize(); ++k) {
156                     Assert.assertEquals(k, c.getPartialDerivativeIndex(c.getPartialDerivativeOrders(k)));
157                 }
158             }
159         }
160     }
161 
162     @Test
163     public void testMultiplicationRules()
164         throws SecurityException, NoSuchFieldException, IllegalArgumentException,
165                IllegalAccessException, NoSuchMethodException, InvocationTargetException {
166 
167         Map<String,String> referenceRules = new HashMap<String, String>();
168         referenceRules.put("(f*g)",             "f * g");
169         referenceRules.put("∂(f*g)/∂p₀",        "f * ∂g/∂p₀ + ∂f/∂p₀ * g");
170         referenceRules.put("∂(f*g)/∂p₁",        referenceRules.get("∂(f*g)/∂p₀").replaceAll("p₀", "p₁"));
171         referenceRules.put("∂(f*g)/∂p₂",        referenceRules.get("∂(f*g)/∂p₀").replaceAll("p₀", "p₂"));
172         referenceRules.put("∂(f*g)/∂p₃",        referenceRules.get("∂(f*g)/∂p₀").replaceAll("p₀", "p₃"));
173         referenceRules.put("∂²(f*g)/∂p₀²",      "f * ∂²g/∂p₀² + 2 * ∂f/∂p₀ * ∂g/∂p₀ + ∂²f/∂p₀² * g");
174         referenceRules.put("∂²(f*g)/∂p₁²",      referenceRules.get("∂²(f*g)/∂p₀²").replaceAll("p₀", "p₁"));
175         referenceRules.put("∂²(f*g)/∂p₂²",      referenceRules.get("∂²(f*g)/∂p₀²").replaceAll("p₀", "p₂"));
176         referenceRules.put("∂²(f*g)/∂p₃²",      referenceRules.get("∂²(f*g)/∂p₀²").replaceAll("p₀", "p₃"));
177         referenceRules.put("∂²(f*g)/∂p₀∂p₁",    "f * ∂²g/∂p₀∂p₁ + ∂f/∂p₁ * ∂g/∂p₀ + ∂f/∂p₀ * ∂g/∂p₁ + ∂²f/∂p₀∂p₁ * g");
178         referenceRules.put("∂²(f*g)/∂p₀∂p₂",    referenceRules.get("∂²(f*g)/∂p₀∂p₁").replaceAll("p₁", "p₂"));
179         referenceRules.put("∂²(f*g)/∂p₀∂p₃",    referenceRules.get("∂²(f*g)/∂p₀∂p₁").replaceAll("p₁", "p₃"));
180         referenceRules.put("∂²(f*g)/∂p₁∂p₂",    referenceRules.get("∂²(f*g)/∂p₀∂p₂").replaceAll("p₀", "p₁"));
181         referenceRules.put("∂²(f*g)/∂p₁∂p₃",    referenceRules.get("∂²(f*g)/∂p₀∂p₃").replaceAll("p₀", "p₁"));
182         referenceRules.put("∂²(f*g)/∂p₂∂p₃",    referenceRules.get("∂²(f*g)/∂p₀∂p₃").replaceAll("p₀", "p₂"));
183         referenceRules.put("∂³(f*g)/∂p₀³",      "f * ∂³g/∂p₀³ +" +
184                                                 " 3 * ∂f/∂p₀ * ∂²g/∂p₀² +" +
185                                                 " 3 * ∂²f/∂p₀² * ∂g/∂p₀ +" +
186                                                 " ∂³f/∂p₀³ * g");
187         referenceRules.put("∂³(f*g)/∂p₁³",      referenceRules.get("∂³(f*g)/∂p₀³").replaceAll("p₀", "p₁"));
188         referenceRules.put("∂³(f*g)/∂p₂³",      referenceRules.get("∂³(f*g)/∂p₀³").replaceAll("p₀", "p₂"));
189         referenceRules.put("∂³(f*g)/∂p₃³",      referenceRules.get("∂³(f*g)/∂p₀³").replaceAll("p₀", "p₃"));
190         referenceRules.put("∂³(f*g)/∂p₀²∂p₁",   "f * ∂³g/∂p₀²∂p₁ +" +
191                                                 " ∂f/∂p₁ * ∂²g/∂p₀² +" +
192                                                 " 2 * ∂f/∂p₀ * ∂²g/∂p₀∂p₁ +" +
193                                                 " 2 * ∂²f/∂p₀∂p₁ * ∂g/∂p₀ +" +
194                                                 " ∂²f/∂p₀² * ∂g/∂p₁ +" +
195                                                 " ∂³f/∂p₀²∂p₁ * g");
196         referenceRules.put("∂³(f*g)/∂p₀∂p₁²",   "f * ∂³g/∂p₀∂p₁² +" +
197                                                 " 2 * ∂f/∂p₁ * ∂²g/∂p₀∂p₁ +" +
198                                                 " ∂²f/∂p₁² * ∂g/∂p₀ +" +
199                                                 " ∂f/∂p₀ * ∂²g/∂p₁² +" +
200                                                 " 2 * ∂²f/∂p₀∂p₁ * ∂g/∂p₁ +" +
201                                                 " ∂³f/∂p₀∂p₁² * g");
202         referenceRules.put("∂³(f*g)/∂p₀²∂p₂",   referenceRules.get("∂³(f*g)/∂p₀²∂p₁").replaceAll("p₁", "p₂"));
203         referenceRules.put("∂³(f*g)/∂p₁²∂p₂",   referenceRules.get("∂³(f*g)/∂p₀²∂p₂").replaceAll("p₀", "p₁"));
204         referenceRules.put("∂³(f*g)/∂p₀∂p₂²",   referenceRules.get("∂³(f*g)/∂p₀∂p₁²").replaceAll("p₁", "p₂"));
205         referenceRules.put("∂³(f*g)/∂p₁∂p₂²",   referenceRules.get("∂³(f*g)/∂p₀∂p₂²").replaceAll("p₀", "p₁"));
206         referenceRules.put("∂³(f*g)/∂p₀²∂p₃",   referenceRules.get("∂³(f*g)/∂p₀²∂p₂").replaceAll("p₂", "p₃"));
207         referenceRules.put("∂³(f*g)/∂p₁²∂p₃",   referenceRules.get("∂³(f*g)/∂p₀²∂p₃").replaceAll("p₀", "p₁"));
208         referenceRules.put("∂³(f*g)/∂p₂²∂p₃",   referenceRules.get("∂³(f*g)/∂p₀²∂p₃").replaceAll("p₀", "p₂"));
209         referenceRules.put("∂³(f*g)/∂p₀∂p₃²",   referenceRules.get("∂³(f*g)/∂p₀∂p₁²").replaceAll("p₁", "p₃"));
210         referenceRules.put("∂³(f*g)/∂p₁∂p₃²",   referenceRules.get("∂³(f*g)/∂p₀∂p₃²").replaceAll("p₀", "p₁"));
211         referenceRules.put("∂³(f*g)/∂p₂∂p₃²",   referenceRules.get("∂³(f*g)/∂p₀∂p₃²").replaceAll("p₀", "p₂"));
212         referenceRules.put("∂³(f*g)/∂p₀∂p₁∂p₂", "f * ∂³g/∂p₀∂p₁∂p₂ +" +
213                                                 " ∂f/∂p₂ * ∂²g/∂p₀∂p₁ +" +
214                                                 " ∂f/∂p₁ * ∂²g/∂p₀∂p₂ +" +
215                                                 " ∂²f/∂p₁∂p₂ * ∂g/∂p₀ +" +
216                                                 " ∂f/∂p₀ * ∂²g/∂p₁∂p₂ +" +
217                                                 " ∂²f/∂p₀∂p₂ * ∂g/∂p₁ +" +
218                                                 " ∂²f/∂p₀∂p₁ * ∂g/∂p₂ +" +
219                                                 " ∂³f/∂p₀∂p₁∂p₂ * g");
220         referenceRules.put("∂³(f*g)/∂p₀∂p₁∂p₃", referenceRules.get("∂³(f*g)/∂p₀∂p₁∂p₂").replaceAll("p₂", "p₃"));
221         referenceRules.put("∂³(f*g)/∂p₀∂p₂∂p₃", referenceRules.get("∂³(f*g)/∂p₀∂p₁∂p₃").replaceAll("p₁", "p₂"));
222         referenceRules.put("∂³(f*g)/∂p₁∂p₂∂p₃", referenceRules.get("∂³(f*g)/∂p₀∂p₂∂p₃").replaceAll("p₀", "p₁"));
223 
224         Field multFieldArrayField = DSCompiler.class.getDeclaredField("multIndirection");
225         multFieldArrayField.setAccessible(true);
226         Class<?> abstractMapperClass = Stream.
227                         of(DSCompiler.class.getDeclaredClasses()).
228                         filter(c -> c.getName().endsWith("AbstractMapper")).
229                         findAny().
230                         get();
231         Class<?> multiplicationMapperClass = Stream.
232                         of(DSCompiler.class.getDeclaredClasses()).
233                         filter(c -> c.getName().endsWith("MultiplicationMapper")).
234                         findAny().
235                         get();
236         Method coeffMethod = abstractMapperClass.getDeclaredMethod("getCoeff");
237         Field lhsField = multiplicationMapperClass.getDeclaredField("lhsIndex");
238         lhsField.setAccessible(true);
239         Field rhsField = multiplicationMapperClass.getDeclaredField("rhsIndex");
240         rhsField.setAccessible(true);
241         for (int i = 0; i < 5; ++i) {
242             for (int j = 0; j < 4; ++j) {
243                 DSCompiler compiler = DSCompiler.getCompiler(i, j);
244                 Object[][] multIndirection = (Object[][]) multFieldArrayField.get(compiler);
245                 for (int k = 0; k < multIndirection.length; ++k) {
246                     String product = ordersToString(compiler.getPartialDerivativeOrders(k),
247                                                     "(f*g)", variables("p"));
248                     StringBuilder rule = new StringBuilder();
249                     for (Object term : multIndirection[k]) {
250                         if (rule.length() > 0) {
251                             rule.append(" + ");
252                         }
253                         if (((Integer) coeffMethod.invoke(term)).intValue() > 1) {
254                             rule.append(((Integer) coeffMethod.invoke(term)).intValue()).append(" * ");
255                         }
256                         rule.append(ordersToString(compiler.getPartialDerivativeOrders(((Integer) lhsField.get(term)).intValue()),
257                                                    "f", variables("p")));
258                         rule.append(" * ");
259                         rule.append(ordersToString(compiler.getPartialDerivativeOrders(((Integer) rhsField.get(term)).intValue()),
260                                                    "g", variables("p")));
261                     }
262                     Assert.assertEquals(product, referenceRules.get(product), rule.toString());
263                 }
264             }
265         }
266     }
267 
268     @Test
269     public void testCompositionRules()
270         throws SecurityException, NoSuchFieldException, IllegalArgumentException,
271                IllegalAccessException, InvocationTargetException, NoSuchMethodException {
272 
273         // the following reference rules have all been computed independently from the library,
274         // using only pencil and paper and some search and replace to handle symmetries
275         Map<String,String> referenceRules = new HashMap<String, String>();
276         referenceRules.put("(f(g))",              "(f(g))");
277         referenceRules.put("∂(f(g))/∂p₀",          "∂(f(g))/∂g * ∂g/∂p₀");
278         referenceRules.put("∂(f(g))/∂p₁",          referenceRules.get("∂(f(g))/∂p₀").replaceAll("p₀", "p₁"));
279         referenceRules.put("∂(f(g))/∂p₂",          referenceRules.get("∂(f(g))/∂p₀").replaceAll("p₀", "p₂"));
280         referenceRules.put("∂(f(g))/∂p₃",          referenceRules.get("∂(f(g))/∂p₀").replaceAll("p₀", "p₃"));
281         referenceRules.put("∂²(f(g))/∂p₀²",        "∂²(f(g))/∂g² * ∂g/∂p₀ * ∂g/∂p₀ + ∂(f(g))/∂g * ∂²g/∂p₀²");
282         referenceRules.put("∂²(f(g))/∂p₁²",        referenceRules.get("∂²(f(g))/∂p₀²").replaceAll("p₀", "p₁"));
283         referenceRules.put("∂²(f(g))/∂p₂²",        referenceRules.get("∂²(f(g))/∂p₀²").replaceAll("p₀", "p₂"));
284         referenceRules.put("∂²(f(g))/∂p₃²",        referenceRules.get("∂²(f(g))/∂p₀²").replaceAll("p₀", "p₃"));
285         referenceRules.put("∂²(f(g))/∂p₀∂p₁",       "∂²(f(g))/∂g² * ∂g/∂p₀ * ∂g/∂p₁ + ∂(f(g))/∂g * ∂²g/∂p₀∂p₁");
286         referenceRules.put("∂²(f(g))/∂p₀∂p₂",       referenceRules.get("∂²(f(g))/∂p₀∂p₁").replaceAll("p₁", "p₂"));
287         referenceRules.put("∂²(f(g))/∂p₀∂p₃",       referenceRules.get("∂²(f(g))/∂p₀∂p₁").replaceAll("p₁", "p₃"));
288         referenceRules.put("∂²(f(g))/∂p₁∂p₂",       referenceRules.get("∂²(f(g))/∂p₀∂p₂").replaceAll("p₀", "p₁"));
289         referenceRules.put("∂²(f(g))/∂p₁∂p₃",       referenceRules.get("∂²(f(g))/∂p₀∂p₃").replaceAll("p₀", "p₁"));
290         referenceRules.put("∂²(f(g))/∂p₂∂p₃",       referenceRules.get("∂²(f(g))/∂p₀∂p₃").replaceAll("p₀", "p₂"));
291         referenceRules.put("∂³(f(g))/∂p₀³",        "∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₀ +" +
292                                                   " 3 * ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂²g/∂p₀² +" +
293                                                   " ∂(f(g))/∂g * ∂³g/∂p₀³");
294         referenceRules.put("∂³(f(g))/∂p₁³",        referenceRules.get("∂³(f(g))/∂p₀³").replaceAll("p₀", "p₁"));
295         referenceRules.put("∂³(f(g))/∂p₂³",        referenceRules.get("∂³(f(g))/∂p₀³").replaceAll("p₀", "p₂"));
296         referenceRules.put("∂³(f(g))/∂p₃³",        referenceRules.get("∂³(f(g))/∂p₀³").replaceAll("p₀", "p₃"));
297         referenceRules.put("∂³(f(g))/∂p₀∂p₁²",      "∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₁ +" +
298                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂²g/∂p₀∂p₁ +" +
299                                                   " ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂²g/∂p₁² +" +
300                                                   " ∂(f(g))/∂g * ∂³g/∂p₀∂p₁²");
301         referenceRules.put("∂³(f(g))/∂p₀∂p₂²",      referenceRules.get("∂³(f(g))/∂p₀∂p₁²").replaceAll("p₁", "p₂"));
302         referenceRules.put("∂³(f(g))/∂p₀∂p₃²",      referenceRules.get("∂³(f(g))/∂p₀∂p₁²").replaceAll("p₁", "p₃"));
303         referenceRules.put("∂³(f(g))/∂p₁∂p₂²",      referenceRules.get("∂³(f(g))/∂p₀∂p₂²").replaceAll("p₀", "p₁"));
304         referenceRules.put("∂³(f(g))/∂p₁∂p₃²",      referenceRules.get("∂³(f(g))/∂p₀∂p₃²").replaceAll("p₀", "p₁"));
305         referenceRules.put("∂³(f(g))/∂p₂∂p₃²",      referenceRules.get("∂³(f(g))/∂p₀∂p₃²").replaceAll("p₀", "p₂"));
306         referenceRules.put("∂³(f(g))/∂p₀²∂p₁",      "∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₁ +" +
307                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂²g/∂p₀∂p₁ +" +
308                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀² * ∂g/∂p₁ +" +
309                                                   " ∂(f(g))/∂g * ∂³g/∂p₀²∂p₁");
310         referenceRules.put("∂³(f(g))/∂p₀²∂p₂",      referenceRules.get("∂³(f(g))/∂p₀²∂p₁").replaceAll("p₁", "p₂"));
311         referenceRules.put("∂³(f(g))/∂p₀²∂p₃",      referenceRules.get("∂³(f(g))/∂p₀²∂p₁").replaceAll("p₁", "p₃"));
312         referenceRules.put("∂³(f(g))/∂p₁²∂p₂",      referenceRules.get("∂³(f(g))/∂p₀²∂p₂").replaceAll("p₀", "p₁"));
313         referenceRules.put("∂³(f(g))/∂p₁²∂p₃",      referenceRules.get("∂³(f(g))/∂p₀²∂p₃").replaceAll("p₀", "p₁"));
314         referenceRules.put("∂³(f(g))/∂p₂²∂p₃",      referenceRules.get("∂³(f(g))/∂p₀²∂p₃").replaceAll("p₀", "p₂"));
315         referenceRules.put("∂³(f(g))/∂p₀∂p₁∂p₂",     "∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₂ +" +
316                                                   " ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂²g/∂p₀∂p₂ +" +
317                                                   " ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂²g/∂p₁∂p₂ +" +
318                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂g/∂p₂ +" +
319                                                   " ∂(f(g))/∂g * ∂³g/∂p₀∂p₁∂p₂");
320         referenceRules.put("∂³(f(g))/∂p₀∂p₁∂p₃",     referenceRules.get("∂³(f(g))/∂p₀∂p₁∂p₂").replaceAll("p₂", "p₃"));
321         referenceRules.put("∂³(f(g))/∂p₀∂p₂∂p₃",     referenceRules.get("∂³(f(g))/∂p₀∂p₁∂p₃").replaceAll("p₁", "p₂"));
322         referenceRules.put("∂³(f(g))/∂p₁∂p₂∂p₃",     referenceRules.get("∂³(f(g))/∂p₀∂p₂∂p₃").replaceAll("p₀", "p₁"));
323         referenceRules.put("∂⁴(f(g))/∂p₀⁴",        "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₀ +" +
324                                                   " 6 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₀ * ∂²g/∂p₀² +" +
325                                                   " 3 * ∂²(f(g))/∂g² * ∂²g/∂p₀² * ∂²g/∂p₀² +" +
326                                                   " 4 * ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₀³ +" +
327                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀⁴");
328         referenceRules.put("∂⁴(f(g))/∂p₁⁴",        referenceRules.get("∂⁴(f(g))/∂p₀⁴").replaceAll("p₀", "p₁"));
329         referenceRules.put("∂⁴(f(g))/∂p₂⁴",        referenceRules.get("∂⁴(f(g))/∂p₀⁴").replaceAll("p₀", "p₂"));
330         referenceRules.put("∂⁴(f(g))/∂p₃⁴",        referenceRules.get("∂⁴(f(g))/∂p₀⁴").replaceAll("p₀", "p₃"));
331         referenceRules.put("∂⁴(f(g))/∂p₀³∂p₁",      "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₁ +" +
332                                                   " 3 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₀ * ∂²g/∂p₀∂p₁ +" +
333                                                   " 3 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂²g/∂p₀² * ∂g/∂p₁ +" +
334                                                   " 3 * ∂²(f(g))/∂g² * ∂²g/∂p₀² * ∂²g/∂p₀∂p₁ +" +
335                                                   " 3 * ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₀²∂p₁ +" +
336                                                   " ∂²(f(g))/∂g² * ∂³g/∂p₀³ * ∂g/∂p₁ +" +
337                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀³∂p₁");
338         referenceRules.put("∂⁴(f(g))/∂p₀³∂p₂",      referenceRules.get("∂⁴(f(g))/∂p₀³∂p₁").replaceAll("p₁", "p₂"));
339         referenceRules.put("∂⁴(f(g))/∂p₀³∂p₃",      referenceRules.get("∂⁴(f(g))/∂p₀³∂p₁").replaceAll("p₁", "p₃"));
340         referenceRules.put("∂⁴(f(g))/∂p₀∂p₁³",      "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₁ * ∂g/∂p₁ +" +
341                                                   " 3 * ∂³(f(g))/∂g³ * ∂g/∂p₁ * ∂g/∂p₁ * ∂²g/∂p₀∂p₁ +" +
342                                                   " 3 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂²g/∂p₁² +" +
343                                                   " 3 * ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂²g/∂p₁² +" +
344                                                   " 3 * ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂³g/∂p₀∂p₁² +" +
345                                                   " ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₁³ +" +
346                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀∂p₁³");
347         referenceRules.put("∂⁴(f(g))/∂p₀∂p₂³",      referenceRules.get("∂⁴(f(g))/∂p₀∂p₁³").replaceAll("p₁", "p₂"));
348         referenceRules.put("∂⁴(f(g))/∂p₀∂p₃³",      referenceRules.get("∂⁴(f(g))/∂p₀∂p₁³").replaceAll("p₁", "p₃"));
349         referenceRules.put("∂⁴(f(g))/∂p₁³∂p₂",      referenceRules.get("∂⁴(f(g))/∂p₀³∂p₂").replaceAll("p₀", "p₁"));
350         referenceRules.put("∂⁴(f(g))/∂p₁³∂p₃",      referenceRules.get("∂⁴(f(g))/∂p₀³∂p₃").replaceAll("p₀", "p₁"));
351         referenceRules.put("∂⁴(f(g))/∂p₁∂p₂³",      referenceRules.get("∂⁴(f(g))/∂p₀∂p₂³").replaceAll("p₀", "p₁"));
352         referenceRules.put("∂⁴(f(g))/∂p₁∂p₃³",      referenceRules.get("∂⁴(f(g))/∂p₀∂p₃³").replaceAll("p₀", "p₁"));
353         referenceRules.put("∂⁴(f(g))/∂p₂³∂p₃",      referenceRules.get("∂⁴(f(g))/∂p₀³∂p₃").replaceAll("p₀", "p₂"));
354         referenceRules.put("∂⁴(f(g))/∂p₂∂p₃³",      referenceRules.get("∂⁴(f(g))/∂p₀∂p₃³").replaceAll("p₀", "p₂"));
355         referenceRules.put("∂⁴(f(g))/∂p₀²∂p₁²",     "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₁ +" +
356                                                   " 4 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂²g/∂p₀∂p₁ +" +
357                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₀ * ∂²g/∂p₁² +" +
358                                                   " 2 * ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂²g/∂p₀∂p₁ +" +
359                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₀∂p₁² +" +
360                                                   " ∂³(f(g))/∂g³ * ∂²g/∂p₀² * ∂g/∂p₁ * ∂g/∂p₁ +" +
361                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂³g/∂p₀²∂p₁ +" +
362                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀² * ∂²g/∂p₁² +" +
363                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀²∂p₁²");
364         referenceRules.put("∂⁴(f(g))/∂p₀²∂p₂²",     referenceRules.get("∂⁴(f(g))/∂p₀²∂p₁²").replaceAll("p₁", "p₂"));
365         referenceRules.put("∂⁴(f(g))/∂p₀²∂p₃²",     referenceRules.get("∂⁴(f(g))/∂p₀²∂p₁²").replaceAll("p₁", "p₃"));
366         referenceRules.put("∂⁴(f(g))/∂p₁²∂p₂²",     referenceRules.get("∂⁴(f(g))/∂p₀²∂p₂²").replaceAll("p₀", "p₁"));
367         referenceRules.put("∂⁴(f(g))/∂p₁²∂p₃²",     referenceRules.get("∂⁴(f(g))/∂p₀²∂p₃²").replaceAll("p₀", "p₁"));
368         referenceRules.put("∂⁴(f(g))/∂p₂²∂p₃²",     referenceRules.get("∂⁴(f(g))/∂p₀²∂p₃²").replaceAll("p₀", "p₂"));
369 
370         referenceRules.put("∂⁴(f(g))/∂p₀²∂p₁∂p₂",    "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₂ +" +
371                                                   " 2 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂²g/∂p₀∂p₂ +" +
372                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₀ * ∂²g/∂p₁∂p₂ +" +
373                                                   " 2 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂²g/∂p₀∂p₁ * ∂g/∂p₂ +" +
374                                                   " 2 * ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂²g/∂p₀∂p₂ +" +
375                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₀∂p₁∂p₂ +" +
376                                                   " ∂³(f(g))/∂g³ * ∂²g/∂p₀² * ∂g/∂p₁ * ∂g/∂p₂ +" +
377                                                   " ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂³g/∂p₀²∂p₂ +" +
378                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀² * ∂²g/∂p₁∂p₂ +" +
379                                                   " ∂²(f(g))/∂g² * ∂³g/∂p₀²∂p₁ * ∂g/∂p₂ +" +
380                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀²∂p₁∂p₂");
381         referenceRules.put("∂⁴(f(g))/∂p₀²∂p₁∂p₃",    referenceRules.get("∂⁴(f(g))/∂p₀²∂p₁∂p₂").replaceAll("p₂", "p₃"));
382         referenceRules.put("∂⁴(f(g))/∂p₀²∂p₂∂p₃",    referenceRules.get("∂⁴(f(g))/∂p₀²∂p₁∂p₃").replaceAll("p₁", "p₂"));
383         referenceRules.put("∂⁴(f(g))/∂p₀∂p₁²∂p₂",    "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₁ * ∂g/∂p₂ +" +
384                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₁ * ∂g/∂p₁ * ∂²g/∂p₀∂p₂ +" +
385                                                   " 2 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂²g/∂p₁∂p₂ +" +
386                                                   " 2 * ∂³(f(g))/∂g³ * ∂g/∂p₁ * ∂²g/∂p₀∂p₁ * ∂g/∂p₂ +" +
387                                                   " 2 * ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂²g/∂p₁∂p₂ +" +
388                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂³g/∂p₀∂p₁∂p₂ +" +
389                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂²g/∂p₁² * ∂g/∂p₂ +" +
390                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₁² * ∂²g/∂p₀∂p₂ +" +
391                                                   " ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₁²∂p₂ +" +
392                                                   " ∂²(f(g))/∂g² * ∂³g/∂p₀∂p₁² * ∂g/∂p₂ +" +
393                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀∂p₁²∂p₂");
394         referenceRules.put("∂⁴(f(g))/∂p₀∂p₁²∂p₃",    referenceRules.get("∂⁴(f(g))/∂p₀∂p₁²∂p₂").replaceAll("p₂", "p₃"));
395         referenceRules.put("∂⁴(f(g))/∂p₁²∂p₂∂p₃",    referenceRules.get("∂⁴(f(g))/∂p₀²∂p₂∂p₃").replaceAll("p₀", "p₁"));
396         referenceRules.put("∂⁴(f(g))/∂p₀∂p₁∂p₂²",    "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₂ * ∂g/∂p₂ +" +
397                                                   " 2 * ∂³(f(g))/∂g³ * ∂g/∂p₁ * ∂g/∂p₂ * ∂²g/∂p₀∂p₂ +" +
398                                                   " 2 * ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₂ * ∂²g/∂p₁∂p₂ +" +
399                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂²g/∂p₂² +" +
400                                                   " 2 * ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₂ * ∂²g/∂p₁∂p₂ +" +
401                                                   " ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂³g/∂p₀∂p₂² +" +
402                                                   " ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₁∂p₂² +" +
403                                                   " ∂³(f(g))/∂g³ * ∂²g/∂p₀∂p₁ * ∂g/∂p₂ * ∂g/∂p₂ +" +
404                                                   " 2 * ∂²(f(g))/∂g² * ∂g/∂p₂ * ∂³g/∂p₀∂p₁∂p₂ +" +
405                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂²g/∂p₂² +" +
406                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀∂p₁∂p₂²");
407         referenceRules.put("∂⁴(f(g))/∂p₀∂p₂²∂p₃",    referenceRules.get("∂⁴(f(g))/∂p₀∂p₁²∂p₃").replaceAll("p₁", "p₂"));
408         referenceRules.put("∂⁴(f(g))/∂p₁∂p₂²∂p₃",    referenceRules.get("∂⁴(f(g))/∂p₀∂p₂²∂p₃").replaceAll("p₀", "p₁"));
409         referenceRules.put("∂⁴(f(g))/∂p₀∂p₁∂p₃²",    referenceRules.get("∂⁴(f(g))/∂p₀∂p₁∂p₂²").replaceAll("p₂", "p₃"));
410         referenceRules.put("∂⁴(f(g))/∂p₀∂p₂∂p₃²",    referenceRules.get("∂⁴(f(g))/∂p₀∂p₁∂p₃²").replaceAll("p₁", "p₂"));
411         referenceRules.put("∂⁴(f(g))/∂p₁∂p₂∂p₃²",    referenceRules.get("∂⁴(f(g))/∂p₀∂p₂∂p₃²").replaceAll("p₀", "p₁"));
412         referenceRules.put("∂⁴(f(g))/∂p₀∂p₁∂p₂∂p₃",   "∂⁴(f(g))/∂g⁴ * ∂g/∂p₀ * ∂g/∂p₁ * ∂g/∂p₂ * ∂g/∂p₃ +" +
413                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₁ * ∂g/∂p₂ * ∂²g/∂p₀∂p₃ +" +
414                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₂ * ∂²g/∂p₁∂p₃ +" +
415                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂g/∂p₁ * ∂²g/∂p₂∂p₃ +" +
416                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₁ * ∂²g/∂p₀∂p₂ * ∂g/∂p₃ +" +
417                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₂ * ∂²g/∂p₁∂p₃ +" +
418                                                   " ∂²(f(g))/∂g² * ∂g/∂p₁ * ∂³g/∂p₀∂p₂∂p₃ +" +
419                                                   " ∂³(f(g))/∂g³ * ∂g/∂p₀ * ∂²g/∂p₁∂p₂ * ∂g/∂p₃ +" +
420                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₁∂p₂ * ∂²g/∂p₀∂p₃ +" +
421                                                   " ∂²(f(g))/∂g² * ∂g/∂p₀ * ∂³g/∂p₁∂p₂∂p₃ +" +
422                                                   " ∂³(f(g))/∂g³ * ∂²g/∂p₀∂p₁ * ∂g/∂p₂ * ∂g/∂p₃ +" +
423                                                   " ∂²(f(g))/∂g² * ∂g/∂p₂ * ∂³g/∂p₀∂p₁∂p₃ +" +
424                                                   " ∂²(f(g))/∂g² * ∂²g/∂p₀∂p₁ * ∂²g/∂p₂∂p₃ +" +
425                                                   " ∂²(f(g))/∂g² * ∂³g/∂p₀∂p₁∂p₂ * ∂g/∂p₃ +" +
426                                                   " ∂(f(g))/∂g * ∂⁴g/∂p₀∂p₁∂p₂∂p₃");
427 
428         Field compFieldArrayField = DSCompiler.class.getDeclaredField("compIndirection");
429         compFieldArrayField.setAccessible(true);
430 
431         for (int i = 0; i < 5; ++i) {
432             for (int j = 0; j < 5; ++j) {
433                 DSCompiler compiler = DSCompiler.getCompiler(i, j);
434                 Object[][] compIndirection = (Object[][]) compFieldArrayField.get(compiler);
435                 for (int k = 0; k < compIndirection.length; ++k) {
436                     String product = ordersToString(compiler.getPartialDerivativeOrders(k),
437                                                     "(f(g))", variables("p"));
438                     String rule = univariateCompositionMappersToString(compiler, compIndirection[k]);
439                     Assert.assertEquals(product, referenceRules.get(product), rule.toString());
440                 }
441             }
442         }
443 
444     }
445 
446     @Test
447     public void testRebaserRules()
448         throws IllegalAccessException, IllegalArgumentException, InvocationTargetException,
449                NoSuchMethodException, SecurityException, NoSuchFieldException {
450 
451         // the following reference rules have all been computed independently from the library,
452         // using only pencil and paper (which was really tedious) and using search and replace to handle symmetries
453         Map<String,String> referenceRules = new HashMap<String, String>();
454         referenceRules.put("f",              "f");
455         referenceRules.put("∂f/∂q₀",         "∂f/∂p₀ ∂p₀/∂q₀ + ∂f/∂p₁ ∂p₁/∂q₀");
456         referenceRules.put("∂f/∂q₁",         referenceRules.get("∂f/∂q₀").replaceAll("q₀", "q₁"));
457         referenceRules.put("∂f/∂q₂",         referenceRules.get("∂f/∂q₀").replaceAll("q₀", "q₂"));
458         referenceRules.put("∂²f/∂q₀²",       "∂²f/∂p₀² (∂p₀/∂q₀)² + 2 ∂²f/∂p₀∂p₁ ∂p₀/∂q₀ ∂p₁/∂q₀" +
459                                              " + ∂f/∂p₀ ∂²p₀/∂q₀² + ∂²f/∂p₁² (∂p₁/∂q₀)² + ∂f/∂p₁ ∂²p₁/∂q₀²");
460         referenceRules.put("∂²f/∂q₁²",       referenceRules.get("∂²f/∂q₀²").replaceAll("q₀", "q₁"));
461         referenceRules.put("∂²f/∂q₂²",       referenceRules.get("∂²f/∂q₀²").replaceAll("q₀", "q₂"));
462         referenceRules.put("∂²f/∂q₀∂q₁",     "∂²f/∂p₀² ∂p₀/∂q₀ ∂p₀/∂q₁ + ∂²f/∂p₀∂p₁ ∂p₀/∂q₁ ∂p₁/∂q₀ + ∂f/∂p₀ ∂²p₀/∂q₀∂q₁" +
463                                              " + ∂²f/∂p₀∂p₁ ∂p₀/∂q₀ ∂p₁/∂q₁ + ∂²f/∂p₁² ∂p₁/∂q₀ ∂p₁/∂q₁ + ∂f/∂p₁ ∂²p₁/∂q₀∂q₁");
464         referenceRules.put("∂²f/∂q₀∂q₂",     referenceRules.get("∂²f/∂q₀∂q₁").replaceAll("q₁", "q₂"));
465         referenceRules.put("∂²f/∂q₁∂q₂",     referenceRules.get("∂²f/∂q₀∂q₂").replaceAll("q₀", "q₁"));
466         referenceRules.put("∂³f/∂q₀³",       "∂³f/∂p₀³ (∂p₀/∂q₀)³ + 3 ∂³f/∂p₀²∂p₁ (∂p₀/∂q₀)² ∂p₁/∂q₀" +
467                                              " + 3 ∂²f/∂p₀² ∂p₀/∂q₀ ∂²p₀/∂q₀² + 3 ∂³f/∂p₀∂p₁² ∂p₀/∂q₀ (∂p₁/∂q₀)²" +
468                                              " + 3 ∂²f/∂p₀∂p₁ ∂²p₀/∂q₀² ∂p₁/∂q₀ + 3 ∂²f/∂p₀∂p₁ ∂p₀/∂q₀ ∂²p₁/∂q₀²" +
469                                              " + ∂f/∂p₀ ∂³p₀/∂q₀³ + ∂³f/∂p₁³ (∂p₁/∂q₀)³" +
470                                              " + 3 ∂²f/∂p₁² ∂p₁/∂q₀ ∂²p₁/∂q₀² + ∂f/∂p₁ ∂³p₁/∂q₀³");
471         referenceRules.put("∂³f/∂q₁³",       referenceRules.get("∂³f/∂q₀³").replaceAll("q₀", "q₁"));
472         referenceRules.put("∂³f/∂q₂³",       referenceRules.get("∂³f/∂q₀³").replaceAll("q₀", "q₂"));
473         referenceRules.put("∂³f/∂q₀²∂q₁",    "∂³f/∂p₀³ (∂p₀/∂q₀)² ∂p₀/∂q₁ + 2 ∂³f/∂p₀²∂p₁ ∂p₀/∂q₀ ∂p₀/∂q₁ ∂p₁/∂q₀" +
474                                              " + ∂²f/∂p₀² ∂²p₀/∂q₀² ∂p₀/∂q₁ + 2 ∂²f/∂p₀² ∂p₀/∂q₀ ∂²p₀/∂q₀∂q₁" +
475                                              " + ∂³f/∂p₀∂p₁² ∂p₀/∂q₁ (∂p₁/∂q₀)² + 2 ∂²f/∂p₀∂p₁ ∂²p₀/∂q₀∂q₁ ∂p₁/∂q₀" +
476                                              " + ∂²f/∂p₀∂p₁ ∂p₀/∂q₁ ∂²p₁/∂q₀² + ∂f/∂p₀ ∂³p₀/∂q₀²∂q₁" +
477                                              " + ∂³f/∂p₀²∂p₁ (∂p₀/∂q₀)² ∂p₁/∂q₁ + 2 ∂³f/∂p₀∂p₁² ∂p₀/∂q₀ ∂p₁/∂q₀ ∂p₁/∂q₁" +
478                                              " + ∂²f/∂p₀∂p₁ ∂²p₀/∂q₀² ∂p₁/∂q₁ + 2 ∂²f/∂p₀∂p₁ ∂p₀/∂q₀ ∂²p₁/∂q₀∂q₁" +
479                                              " + ∂³f/∂p₁³ (∂p₁/∂q₀)² ∂p₁/∂q₁ + ∂²f/∂p₁² ∂²p₁/∂q₀² ∂p₁/∂q₁" +
480                                              " + 2 ∂²f/∂p₁² ∂p₁/∂q₀ ∂²p₁/∂q₀∂q₁ + ∂f/∂p₁ ∂³p₁/∂q₀²∂q₁");
481         referenceRules.put("∂³f/∂q₀²∂q₂",    referenceRules.get("∂³f/∂q₀²∂q₁").replaceAll("q₁", "q₂"));
482         referenceRules.put("∂³f/∂q₁²∂q₂",    referenceRules.get("∂³f/∂q₀²∂q₂").replaceAll("q₀", "q₁"));
483         referenceRules.put("∂³f/∂q₁²∂q₀",    referenceRules.get("∂³f/∂q₁²∂q₂").replaceAll("q₂", "q₀"));
484         referenceRules.put("∂³f/∂q₀∂q₁²",    "∂³f/∂p₀³ ∂p₀/∂q₀ (∂p₀/∂q₁)² + ∂³f/∂p₀²∂p₁ (∂p₀/∂q₁)² ∂p₁/∂q₀" +
485                                              " + 2 ∂²f/∂p₀² ∂p₀/∂q₁ ∂²p₀/∂q₀∂q₁ + 2 ∂³f/∂p₀²∂p₁ ∂p₀/∂q₀ ∂p₀/∂q₁ ∂p₁/∂q₁" +
486                                              " + 2 ∂³f/∂p₀∂p₁² ∂p₀/∂q₁ ∂p₁/∂q₀ ∂p₁/∂q₁ + 2 ∂²f/∂p₀∂p₁ ∂²p₀/∂q₀∂q₁ ∂p₁/∂q₁" +
487                                              " + 2 ∂²f/∂p₀∂p₁ ∂p₀/∂q₁ ∂²p₁/∂q₀∂q₁ + ∂²f/∂p₀² ∂p₀/∂q₀ ∂²p₀/∂q₁²" +
488                                              " + ∂²f/∂p₀∂p₁ ∂²p₀/∂q₁² ∂p₁/∂q₀ + ∂f/∂p₀ ∂³p₀/∂q₀∂q₁²" +
489                                              " + ∂³f/∂p₀∂p₁² ∂p₀/∂q₀ (∂p₁/∂q₁)² + ∂³f/∂p₁³ ∂p₁/∂q₀ (∂p₁/∂q₁)²" +
490                                              " + 2 ∂²f/∂p₁² ∂p₁/∂q₁ ∂²p₁/∂q₀∂q₁ + ∂²f/∂p₀∂p₁ ∂p₀/∂q₀ ∂²p₁/∂q₁²" +
491                                              " + ∂²f/∂p₁² ∂p₁/∂q₀ ∂²p₁/∂q₁²" +
492                                              " + ∂f/∂p₁ ∂³p₁/∂q₀∂q₁²");
493         referenceRules.put("∂³f/∂q₀∂q₂²",   referenceRules.get("∂³f/∂q₀∂q₁²").replaceAll("q₁", "q₂"));
494         referenceRules.put("∂³f/∂q₁∂q₂²",   referenceRules.get("∂³f/∂q₀∂q₂²").replaceAll("q₀", "q₁"));
495         referenceRules.put("∂³f/∂q₀∂q₁∂q₂", "∂³f/∂p₀³ ∂p₀/∂q₀ ∂p₀/∂q₁ ∂p₀/∂q₂ + ∂³f/∂p₀²∂p₁ ∂p₀/∂q₁ ∂p₀/∂q₂ ∂p₁/∂q₀" +
496                                             " + ∂²f/∂p₀² ∂²p₀/∂q₀∂q₁ ∂p₀/∂q₂ + ∂²f/∂p₀² ∂p₀/∂q₁ ∂²p₀/∂q₀∂q₂" +
497                                             " + ∂³f/∂p₀²∂p₁ ∂p₀/∂q₀ ∂p₀/∂q₂ ∂p₁/∂q₁ + ∂³f/∂p₀∂p₁² ∂p₀/∂q₂ ∂p₁/∂q₀ ∂p₁/∂q₁" +
498                                             " + ∂²f/∂p₀∂p₁ ∂²p₀/∂q₀∂q₂ ∂p₁/∂q₁ + ∂²f/∂p₀∂p₁ ∂p₀/∂q₂ ∂²p₁/∂q₀∂q₁" +
499                                             " + ∂²f/∂p₀² ∂p₀/∂q₀ ∂²p₀/∂q₁∂q₂ + ∂²f/∂p₀∂p₁ ∂²p₀/∂q₁∂q₂ ∂p₁/∂q₀" +
500                                             " + ∂f/∂p₀ ∂³p₀/∂q₀∂q₁∂q₂ + ∂³f/∂p₀²∂p₁ ∂p₀/∂q₀ ∂p₀/∂q₁ ∂p₁/∂q₂" +
501                                             " + ∂³f/∂p₀∂p₁² ∂p₀/∂q₁ ∂p₁/∂q₀ ∂p₁/∂q₂ + ∂²f/∂p₀∂p₁ ∂²p₀/∂q₀∂q₁ ∂p₁/∂q₂" +
502                                             " + ∂²f/∂p₀∂p₁ ∂p₀/∂q₁ ∂²p₁/∂q₀∂q₂ + ∂³f/∂p₀∂p₁² ∂p₀/∂q₀ ∂p₁/∂q₁ ∂p₁/∂q₂" +
503                                             " + ∂³f/∂p₁³ ∂p₁/∂q₀ ∂p₁/∂q₁ ∂p₁/∂q₂ + ∂²f/∂p₁² ∂²p₁/∂q₀∂q₁ ∂p₁/∂q₂" +
504                                             " + ∂²f/∂p₁² ∂p₁/∂q₁ ∂²p₁/∂q₀∂q₂ + ∂²f/∂p₀∂p₁ ∂p₀/∂q₀ ∂²p₁/∂q₁∂q₂" +
505                                             " + ∂²f/∂p₁² ∂p₁/∂q₀ ∂²p₁/∂q₁∂q₂ + ∂f/∂p₁ ∂³p₁/∂q₀∂q₁∂q₂");
506 
507         Method getterMethod = DSCompiler.class.getDeclaredMethod("getRebaser", DSCompiler.class);
508         getterMethod.setAccessible(true);
509 
510         for (int order = 0; order < 4; ++order) {
511 
512             // assuming f = f(p₀, p₁)
513             //          p₀ = p₀(q₀, q₁, q₂)
514             //          p₁ = p₁(q₀, q₁, q₂)
515             DSCompiler c2 = DSCompiler.getCompiler(2, order);
516             DSCompiler c3 = DSCompiler.getCompiler(3, order);
517             Object[][] rebaser = (Object[][]) getterMethod.invoke(c2, c3);
518 
519             Assert.assertEquals(c3.getSize(), rebaser.length);
520             for (int k = 0; k < rebaser.length; ++k) {
521                 String key  = ordersToString(c3.getPartialDerivativeOrders(k), "f", variables("q"));
522                 String rule = multivariateCompositionMappersToString(c2, c3, rebaser[k]);
523                 Assert.assertEquals(referenceRules.get(key), rule);
524             }
525 
526         }
527 
528     }
529 
530     private void checkIndices(int[] indices, int ... expected) {
531         Assert.assertEquals(expected.length, indices.length);
532         for (int i = 0; i < expected.length; ++i) {
533             Assert.assertEquals(expected[i], indices[i]);
534         }
535     }
536 
537     private String orderToString(int order, String functionName, String parameterName) {
538         if (order == 0) {
539             return functionName;
540         } else if (order == 1) {
541             return "∂" + functionName + "/∂" + parameterName;
542         } else {
543             return "∂" + exponent(order) + functionName + "/∂" + parameterName + exponent(order);
544         }
545     }
546 
547     private String ordersToString(int[] orders, String functionName, String ... parametersNames) {
548 
549         int sumOrders = 0;
550         for (int order : orders) {
551             sumOrders += order;
552         }
553 
554         if (sumOrders == 0) {
555             return functionName;
556         }
557 
558         StringBuilder builder = new StringBuilder();
559         builder.append('∂');
560         if (sumOrders > 1) {
561             builder.append(exponent(sumOrders));
562         }
563         builder.append(functionName).append('/');
564         for (int i = 0; i < orders.length; ++i) {
565             if (orders[i] > 0) {
566                 builder.append('∂').append(parametersNames[i]);
567                 if (orders[i] > 1) {
568                     builder.append(exponent(orders[i]));
569                 }
570             }
571         }
572         return builder.toString();
573 
574     }
575 
576     private String univariateCompositionMappersToString(final DSCompiler compiler,  final Object[] mappers) {
577          try {
578              
579              Class<?> abstractMapperClass = Stream.
580                              of(DSCompiler.class.getDeclaredClasses()).
581                              filter(c -> c.getName().endsWith("AbstractMapper")).
582                              findAny().
583                              get();
584              Class<?> univariateCompositionMapperClass = Stream.
585                              of(DSCompiler.class.getDeclaredClasses()).
586                              filter(c -> c.getName().endsWith("UnivariateCompositionMapper")).
587                              findAny().
588                              get();
589              Method coeffMethod = abstractMapperClass.getDeclaredMethod("getCoeff");
590              Field fIndexField = univariateCompositionMapperClass.getDeclaredField("fIndex");
591              fIndexField.setAccessible(true);
592              Field dsIndicesField = univariateCompositionMapperClass.getDeclaredField("dsIndices");
593              dsIndicesField.setAccessible(true);
594 
595              StringBuilder rule = new StringBuilder();
596              for (Object term : mappers) {
597                  if (rule.length() > 0) {
598                      rule.append(" + ");
599                  }
600                  if (((Integer) coeffMethod.invoke(term)).intValue() > 1) {
601                      rule.append(((Integer) coeffMethod.invoke(term)).intValue()).append(" * ");
602                  }
603                  rule.append(orderToString(((Integer) fIndexField.get(term)).intValue(), "(f(g))", "g"));
604                  int[] dsIndex = (int[]) dsIndicesField.get(term);
605                  for (int l = 0; l < dsIndex.length; ++l) {
606                      rule.append(" * ");
607                      rule.append(ordersToString(compiler.getPartialDerivativeOrders(dsIndex[l]),
608                                                 "g", "p₀", "p₁", "p₂", "p₃"));
609                  }
610              }
611              return rule.toString();
612 
613          } catch (NoSuchMethodException | SecurityException | NoSuchFieldException | IllegalAccessException |
614                   IllegalArgumentException | InvocationTargetException e) {
615              Assert.fail(e.getLocalizedMessage());
616              return null;
617          }
618      }
619 
620     private String multivariateCompositionMappersToString(final DSCompiler compiler, final DSCompiler baseCompiler,
621                                                           final Object[] mappers) {
622          try {
623              Class<?> abstractMapperClass = Stream.
624                              of(DSCompiler.class.getDeclaredClasses()).
625                              filter(c -> c.getName().endsWith("AbstractMapper")).
626                              findAny().
627                              get();
628              Class<?> multivariateCompositionMapperClass = Stream.
629                              of(DSCompiler.class.getDeclaredClasses()).
630                              filter(c -> c.getName().endsWith("MultivariateCompositionMapper")).
631                              findAny().
632                              get();
633              Method coeffMethod = abstractMapperClass.getDeclaredMethod("getCoeff");
634              Field dsIndexField = multivariateCompositionMapperClass.getDeclaredField("dsIndex");
635              dsIndexField.setAccessible(true);
636              Field productIndicesField = multivariateCompositionMapperClass.getDeclaredField("productIndices");
637              productIndicesField.setAccessible(true);
638 
639              StringBuilder rule = new StringBuilder();
640              for (int i = 0; i < mappers.length; ++i) {
641                  if (i > 0) {
642                      rule.append(" + ");
643                  }
644                  final int coeff = ((Integer) coeffMethod.invoke(mappers[i])).intValue();
645                  if (coeff > 1) {
646                      rule.append(coeff);
647                      rule.append(' ');
648                  }
649                  final int dsIndex = dsIndexField.getInt(mappers[i]);
650                  rule.append(ordersToString(compiler.getPartialDerivativeOrders(dsIndex),
651                                             "f", variables("p")));
652                  final int[] productIndices = (int[]) productIndicesField.get(mappers[i]);
653                  int j = 0;
654                  while (j < productIndices.length) {
655                      int count = 1;
656                      while (j + count < productIndices.length && productIndices[j + count] == productIndices[j]) {
657                          ++count;
658                      }
659                      final int varIndex   = productIndices[j] / baseCompiler.getSize();
660                      final int varDSIndex = productIndices[j] % baseCompiler.getSize();
661                      rule.append(' ');
662                      if (count > 1) {
663                          rule.append('(');
664                      }
665                      rule.append(ordersToString(baseCompiler.getPartialDerivativeOrders(varDSIndex),
666                                                 variables("p")[varIndex], variables("q")));
667                      if (count > 1) {
668                          rule.append(')');
669                          rule.append(exponent(count));
670                      }
671                      j += count;
672                  }
673              }
674              return rule.toString();
675 
676          } catch (NoSuchMethodException | SecurityException | NoSuchFieldException | IllegalAccessException |
677                   IllegalArgumentException | InvocationTargetException e) {
678              Assert.fail(e.getLocalizedMessage());
679              return null;
680          }
681      }
682 
683      private String[] variables(final String baseName) {
684          return new String[] {
685              baseName + "₀", baseName + "₁", baseName + "₂", baseName + "₃", baseName + "₄",
686              baseName + "₅", baseName + "₆", baseName + "₇", baseName + "₈", baseName + "₉"
687          };
688      }
689 
690     private String exponent(int e) {
691         switch (e) {
692             case 0 : return "";
693             case 1 : return "";
694             case 2 : return "²";
695             case 3 : return "³";
696             case 4 : return "⁴";
697             case 5 : return "⁵";
698             case 6 : return "⁶";
699             case 7 : return "⁷";
700             case 8 : return "⁸";
701             case 9 : return "⁹";
702             default:
703                 Assert.fail("exponent out of range");
704                 return null;
705         }
706     }
707 
708 }