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