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    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  package org.apache.commons.math4.legacy.analysis.differentiation;
19  
20  import java.lang.reflect.Field;
21  import java.util.HashMap;
22  import java.util.Map;
23  
24  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
25  import org.apache.commons.numbers.combinatorics.BinomialCoefficient;
26  import org.junit.Assert;
27  import org.junit.Test;
28  
29  
30  /**
31   * Test for class {@link DSCompiler}.
32   */
33  public class DSCompilerTest {
34  
35      @Test
36      public void testSize() {
37          for (int i = 0; i < 6; ++i) {
38              for (int j = 0; j < 6; ++j) {
39                  long expected = BinomialCoefficient.value(i + j, i);
40                  Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
41                  Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
42              }
43          }
44      }
45  
46      @Test
47      public void testIndices() {
48  
49          DSCompiler c = DSCompiler.getCompiler(0, 0);
50          checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
51  
52          c = DSCompiler.getCompiler(0, 1);
53          checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
54  
55          c = DSCompiler.getCompiler(1, 0);
56          checkIndices(c.getPartialDerivativeOrders(0), 0);
57  
58          c = DSCompiler.getCompiler(1, 1);
59          checkIndices(c.getPartialDerivativeOrders(0), 0);
60          checkIndices(c.getPartialDerivativeOrders(1), 1);
61  
62          c = DSCompiler.getCompiler(1, 2);
63          checkIndices(c.getPartialDerivativeOrders(0), 0);
64          checkIndices(c.getPartialDerivativeOrders(1), 1);
65          checkIndices(c.getPartialDerivativeOrders(2), 2);
66  
67          c = DSCompiler.getCompiler(2, 1);
68          checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
69          checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
70          checkIndices(c.getPartialDerivativeOrders(2), 0, 1);
71  
72          c = DSCompiler.getCompiler(1, 3);
73          checkIndices(c.getPartialDerivativeOrders(0), 0);
74          checkIndices(c.getPartialDerivativeOrders(1), 1);
75          checkIndices(c.getPartialDerivativeOrders(2), 2);
76          checkIndices(c.getPartialDerivativeOrders(3), 3);
77  
78          c = DSCompiler.getCompiler(2, 2);
79          checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
80          checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
81          checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
82          checkIndices(c.getPartialDerivativeOrders(3), 0, 1);
83          checkIndices(c.getPartialDerivativeOrders(4), 1, 1);
84          checkIndices(c.getPartialDerivativeOrders(5), 0, 2);
85  
86          c = DSCompiler.getCompiler(3, 1);
87          checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
88          checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
89          checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0);
90          checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1);
91  
92          c = DSCompiler.getCompiler(1, 4);
93          checkIndices(c.getPartialDerivativeOrders(0), 0);
94          checkIndices(c.getPartialDerivativeOrders(1), 1);
95          checkIndices(c.getPartialDerivativeOrders(2), 2);
96          checkIndices(c.getPartialDerivativeOrders(3), 3);
97          checkIndices(c.getPartialDerivativeOrders(4), 4);
98  
99          c = DSCompiler.getCompiler(2, 3);
100         checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
101         checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
102         checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
103         checkIndices(c.getPartialDerivativeOrders(3), 3, 0);
104         checkIndices(c.getPartialDerivativeOrders(4), 0, 1);
105         checkIndices(c.getPartialDerivativeOrders(5), 1, 1);
106         checkIndices(c.getPartialDerivativeOrders(6), 2, 1);
107         checkIndices(c.getPartialDerivativeOrders(7), 0, 2);
108         checkIndices(c.getPartialDerivativeOrders(8), 1, 2);
109         checkIndices(c.getPartialDerivativeOrders(9), 0, 3);
110 
111         c = DSCompiler.getCompiler(3, 2);
112         checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
113         checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
114         checkIndices(c.getPartialDerivativeOrders(2), 2, 0, 0);
115         checkIndices(c.getPartialDerivativeOrders(3), 0, 1, 0);
116         checkIndices(c.getPartialDerivativeOrders(4), 1, 1, 0);
117         checkIndices(c.getPartialDerivativeOrders(5), 0, 2, 0);
118         checkIndices(c.getPartialDerivativeOrders(6), 0, 0, 1);
119         checkIndices(c.getPartialDerivativeOrders(7), 1, 0, 1);
120         checkIndices(c.getPartialDerivativeOrders(8), 0, 1, 1);
121         checkIndices(c.getPartialDerivativeOrders(9), 0, 0, 2);
122 
123         c = DSCompiler.getCompiler(4, 1);
124         checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0, 0);
125         checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0, 0);
126         checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0, 0);
127         checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1, 0);
128         checkIndices(c.getPartialDerivativeOrders(4), 0, 0, 0, 1);
129     }
130 
131     @Test(expected=DimensionMismatchException.class)
132     public void testIncompatibleParams() {
133         DSCompiler.getCompiler(3, 2).checkCompatibility(DSCompiler.getCompiler(4, 2));
134     }
135 
136     @Test(expected=DimensionMismatchException.class)
137     public void testIncompatibleOrder() {
138         DSCompiler.getCompiler(3, 3).checkCompatibility(DSCompiler.getCompiler(3, 2));
139     }
140 
141     @Test
142     public void testSymmetry() {
143         for (int i = 0; i < 6; ++i) {
144             for (int j = 0; j < 6; ++j) {
145                 DSCompiler c = DSCompiler.getCompiler(i, j);
146                 for (int k = 0; k < c.getSize(); ++k) {
147                     Assert.assertEquals(k, c.getPartialDerivativeIndex(c.getPartialDerivativeOrders(k)));
148                 }
149             }
150         }
151     }
152 
153     @Test public void testMultiplicationRules()
154         throws SecurityException, NoSuchFieldException, IllegalArgumentException, IllegalAccessException {
155 
156         Map<String,String> referenceRules = new HashMap<>();
157         referenceRules.put("(f*g)",          "f * g");
158         referenceRules.put("d(f*g)/dx",      "f * dg/dx + df/dx * g");
159         referenceRules.put("d(f*g)/dy",      referenceRules.get("d(f*g)/dx").replaceAll("x", "y"));
160         referenceRules.put("d(f*g)/dz",      referenceRules.get("d(f*g)/dx").replaceAll("x", "z"));
161         referenceRules.put("d(f*g)/dt",      referenceRules.get("d(f*g)/dx").replaceAll("x", "t"));
162         referenceRules.put("d2(f*g)/dx2",    "f * d2g/dx2 + 2 * df/dx * dg/dx + d2f/dx2 * g");
163         referenceRules.put("d2(f*g)/dy2",    referenceRules.get("d2(f*g)/dx2").replaceAll("x", "y"));
164         referenceRules.put("d2(f*g)/dz2",    referenceRules.get("d2(f*g)/dx2").replaceAll("x", "z"));
165         referenceRules.put("d2(f*g)/dt2",    referenceRules.get("d2(f*g)/dx2").replaceAll("x", "t"));
166         referenceRules.put("d2(f*g)/dxdy",   "f * d2g/dxdy + df/dy * dg/dx + df/dx * dg/dy + d2f/dxdy * g");
167         referenceRules.put("d2(f*g)/dxdz",   referenceRules.get("d2(f*g)/dxdy").replaceAll("y", "z"));
168         referenceRules.put("d2(f*g)/dxdt",   referenceRules.get("d2(f*g)/dxdy").replaceAll("y", "t"));
169         referenceRules.put("d2(f*g)/dydz",   referenceRules.get("d2(f*g)/dxdz").replaceAll("x", "y"));
170         referenceRules.put("d2(f*g)/dydt",   referenceRules.get("d2(f*g)/dxdt").replaceAll("x", "y"));
171         referenceRules.put("d2(f*g)/dzdt",   referenceRules.get("d2(f*g)/dxdt").replaceAll("x", "z"));
172         referenceRules.put("d3(f*g)/dx3",    "f * d3g/dx3 +" +
173                                              " 3 * df/dx * d2g/dx2 +" +
174                                              " 3 * d2f/dx2 * dg/dx +" +
175                                              " d3f/dx3 * g");
176         referenceRules.put("d3(f*g)/dy3",   referenceRules.get("d3(f*g)/dx3").replaceAll("x", "y"));
177         referenceRules.put("d3(f*g)/dz3",   referenceRules.get("d3(f*g)/dx3").replaceAll("x", "z"));
178         referenceRules.put("d3(f*g)/dt3",   referenceRules.get("d3(f*g)/dx3").replaceAll("x", "t"));
179         referenceRules.put("d3(f*g)/dx2dy",  "f * d3g/dx2dy +" +
180                                              " df/dy * d2g/dx2 +" +
181                                              " 2 * df/dx * d2g/dxdy +" +
182                                              " 2 * d2f/dxdy * dg/dx +" +
183                                              " d2f/dx2 * dg/dy +" +
184                                              " d3f/dx2dy * g");
185         referenceRules.put("d3(f*g)/dxdy2",  "f * d3g/dxdy2 +" +
186                                              " 2 * df/dy * d2g/dxdy +" +
187                                              " d2f/dy2 * dg/dx +" +
188                                              " df/dx * d2g/dy2 +" +
189                                              " 2 * d2f/dxdy * dg/dy +" +
190                                              " d3f/dxdy2 * g");
191         referenceRules.put("d3(f*g)/dx2dz",   referenceRules.get("d3(f*g)/dx2dy").replaceAll("y", "z"));
192         referenceRules.put("d3(f*g)/dy2dz",   referenceRules.get("d3(f*g)/dx2dz").replaceAll("x", "y"));
193         referenceRules.put("d3(f*g)/dxdz2",   referenceRules.get("d3(f*g)/dxdy2").replaceAll("y", "z"));
194         referenceRules.put("d3(f*g)/dydz2",   referenceRules.get("d3(f*g)/dxdz2").replaceAll("x", "y"));
195         referenceRules.put("d3(f*g)/dx2dt",   referenceRules.get("d3(f*g)/dx2dz").replaceAll("z", "t"));
196         referenceRules.put("d3(f*g)/dy2dt",   referenceRules.get("d3(f*g)/dx2dt").replaceAll("x", "y"));
197         referenceRules.put("d3(f*g)/dz2dt",   referenceRules.get("d3(f*g)/dx2dt").replaceAll("x", "z"));
198         referenceRules.put("d3(f*g)/dxdt2",   referenceRules.get("d3(f*g)/dxdy2").replaceAll("y", "t"));
199         referenceRules.put("d3(f*g)/dydt2",   referenceRules.get("d3(f*g)/dxdt2").replaceAll("x", "y"));
200         referenceRules.put("d3(f*g)/dzdt2",   referenceRules.get("d3(f*g)/dxdt2").replaceAll("x", "z"));
201         referenceRules.put("d3(f*g)/dxdydz", "f * d3g/dxdydz +" +
202                                              " df/dz * d2g/dxdy +" +
203                                              " df/dy * d2g/dxdz +" +
204                                              " d2f/dydz * dg/dx +" +
205                                              " df/dx * d2g/dydz +" +
206                                              " d2f/dxdz * dg/dy +" +
207                                              " d2f/dxdy * dg/dz +" +
208                                              " d3f/dxdydz * g");
209         referenceRules.put("d3(f*g)/dxdydt", referenceRules.get("d3(f*g)/dxdydz").replaceAll("z", "t"));
210         referenceRules.put("d3(f*g)/dxdzdt", referenceRules.get("d3(f*g)/dxdydt").replaceAll("y", "z"));
211         referenceRules.put("d3(f*g)/dydzdt", referenceRules.get("d3(f*g)/dxdzdt").replaceAll("x", "y"));
212 
213         Field multFieldArrayField = DSCompiler.class.getDeclaredField("multIndirection");
214         multFieldArrayField.setAccessible(true);
215         for (int i = 0; i < 5; ++i) {
216             for (int j = 0; j < 4; ++j) {
217                 DSCompiler compiler = DSCompiler.getCompiler(i, j);
218                 int[][][] multIndirection = (int[][][]) multFieldArrayField.get(compiler);
219                 for (int k = 0; k < multIndirection.length; ++k) {
220                     String product = ordersToString(compiler.getPartialDerivativeOrders(k),
221                                                     "(f*g)", "x", "y", "z", "t");
222                     StringBuilder rule = new StringBuilder();
223                     for (int[] term : multIndirection[k]) {
224                         if (rule.length() > 0) {
225                             rule.append(" + ");
226                         }
227                         if (term[0] > 1) {
228                             rule.append(term[0]).append(" * ");
229                         }
230                         rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[1]),
231                                                    "f", "x", "y", "z", "t"));
232                         rule.append(" * ");
233                         rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[2]),
234                                                    "g", "x", "y", "z", "t"));
235                     }
236                     Assert.assertEquals(product, referenceRules.get(product), rule.toString());
237                 }
238             }
239         }
240     }
241 
242     @Test public void testCompositionRules()
243         throws SecurityException, NoSuchFieldException, IllegalArgumentException, IllegalAccessException {
244 
245         // the following reference rules have all been computed independently from the library,
246         // using only pencil and paper and some search and replace to handle symmetries
247         Map<String,String> referenceRules = new HashMap<>();
248         referenceRules.put("(f(g))",              "(f(g))");
249         referenceRules.put("d(f(g))/dx",          "d(f(g))/dg * dg/dx");
250         referenceRules.put("d(f(g))/dy",          referenceRules.get("d(f(g))/dx").replaceAll("x", "y"));
251         referenceRules.put("d(f(g))/dz",          referenceRules.get("d(f(g))/dx").replaceAll("x", "z"));
252         referenceRules.put("d(f(g))/dt",          referenceRules.get("d(f(g))/dx").replaceAll("x", "t"));
253         referenceRules.put("d2(f(g))/dx2",        "d2(f(g))/dg2 * dg/dx * dg/dx + d(f(g))/dg * d2g/dx2");
254         referenceRules.put("d2(f(g))/dy2",        referenceRules.get("d2(f(g))/dx2").replaceAll("x", "y"));
255         referenceRules.put("d2(f(g))/dz2",        referenceRules.get("d2(f(g))/dx2").replaceAll("x", "z"));
256         referenceRules.put("d2(f(g))/dt2",        referenceRules.get("d2(f(g))/dx2").replaceAll("x", "t"));
257         referenceRules.put("d2(f(g))/dxdy",       "d2(f(g))/dg2 * dg/dx * dg/dy + d(f(g))/dg * d2g/dxdy");
258         referenceRules.put("d2(f(g))/dxdz",       referenceRules.get("d2(f(g))/dxdy").replaceAll("y", "z"));
259         referenceRules.put("d2(f(g))/dxdt",       referenceRules.get("d2(f(g))/dxdy").replaceAll("y", "t"));
260         referenceRules.put("d2(f(g))/dydz",       referenceRules.get("d2(f(g))/dxdz").replaceAll("x", "y"));
261         referenceRules.put("d2(f(g))/dydt",       referenceRules.get("d2(f(g))/dxdt").replaceAll("x", "y"));
262         referenceRules.put("d2(f(g))/dzdt",       referenceRules.get("d2(f(g))/dxdt").replaceAll("x", "z"));
263         referenceRules.put("d3(f(g))/dx3",        "d3(f(g))/dg3 * dg/dx * dg/dx * dg/dx +" +
264                                                   " 3 * d2(f(g))/dg2 * dg/dx * d2g/dx2 +" +
265                                                   " d(f(g))/dg * d3g/dx3");
266         referenceRules.put("d3(f(g))/dy3",        referenceRules.get("d3(f(g))/dx3").replaceAll("x", "y"));
267         referenceRules.put("d3(f(g))/dz3",        referenceRules.get("d3(f(g))/dx3").replaceAll("x", "z"));
268         referenceRules.put("d3(f(g))/dt3",        referenceRules.get("d3(f(g))/dx3").replaceAll("x", "t"));
269         referenceRules.put("d3(f(g))/dxdy2",      "d3(f(g))/dg3 * dg/dx * dg/dy * dg/dy +" +
270                                                   " 2 * d2(f(g))/dg2 * dg/dy * d2g/dxdy +" +
271                                                   " d2(f(g))/dg2 * dg/dx * d2g/dy2 +" +
272                                                   " d(f(g))/dg * d3g/dxdy2");
273         referenceRules.put("d3(f(g))/dxdz2",      referenceRules.get("d3(f(g))/dxdy2").replaceAll("y", "z"));
274         referenceRules.put("d3(f(g))/dxdt2",      referenceRules.get("d3(f(g))/dxdy2").replaceAll("y", "t"));
275         referenceRules.put("d3(f(g))/dydz2",      referenceRules.get("d3(f(g))/dxdz2").replaceAll("x", "y"));
276         referenceRules.put("d3(f(g))/dydt2",      referenceRules.get("d3(f(g))/dxdt2").replaceAll("x", "y"));
277         referenceRules.put("d3(f(g))/dzdt2",      referenceRules.get("d3(f(g))/dxdt2").replaceAll("x", "z"));
278         referenceRules.put("d3(f(g))/dx2dy",      "d3(f(g))/dg3 * dg/dx * dg/dx * dg/dy +" +
279                                                   " 2 * d2(f(g))/dg2 * dg/dx * d2g/dxdy +" +
280                                                   " d2(f(g))/dg2 * d2g/dx2 * dg/dy +" +
281                                                   " d(f(g))/dg * d3g/dx2dy");
282         referenceRules.put("d3(f(g))/dx2dz",      referenceRules.get("d3(f(g))/dx2dy").replaceAll("y", "z"));
283         referenceRules.put("d3(f(g))/dx2dt",      referenceRules.get("d3(f(g))/dx2dy").replaceAll("y", "t"));
284         referenceRules.put("d3(f(g))/dy2dz",      referenceRules.get("d3(f(g))/dx2dz").replaceAll("x", "y"));
285         referenceRules.put("d3(f(g))/dy2dt",      referenceRules.get("d3(f(g))/dx2dt").replaceAll("x", "y"));
286         referenceRules.put("d3(f(g))/dz2dt",      referenceRules.get("d3(f(g))/dx2dt").replaceAll("x", "z"));
287         referenceRules.put("d3(f(g))/dxdydz",     "d3(f(g))/dg3 * dg/dx * dg/dy * dg/dz +" +
288                                                   " d2(f(g))/dg2 * dg/dy * d2g/dxdz +" +
289                                                   " d2(f(g))/dg2 * dg/dx * d2g/dydz +" +
290                                                   " d2(f(g))/dg2 * d2g/dxdy * dg/dz +" +
291                                                   " d(f(g))/dg * d3g/dxdydz");
292         referenceRules.put("d3(f(g))/dxdydt",     referenceRules.get("d3(f(g))/dxdydz").replaceAll("z", "t"));
293         referenceRules.put("d3(f(g))/dxdzdt",     referenceRules.get("d3(f(g))/dxdydt").replaceAll("y", "z"));
294         referenceRules.put("d3(f(g))/dydzdt",     referenceRules.get("d3(f(g))/dxdzdt").replaceAll("x", "y"));
295         referenceRules.put("d4(f(g))/dx4",        "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dx * dg/dx +" +
296                                                   " 6 * d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dx2 +" +
297                                                   " 3 * d2(f(g))/dg2 * d2g/dx2 * d2g/dx2 +" +
298                                                   " 4 * d2(f(g))/dg2 * dg/dx * d3g/dx3 +" +
299                                                   " d(f(g))/dg * d4g/dx4");
300         referenceRules.put("d4(f(g))/dy4",        referenceRules.get("d4(f(g))/dx4").replaceAll("x", "y"));
301         referenceRules.put("d4(f(g))/dz4",        referenceRules.get("d4(f(g))/dx4").replaceAll("x", "z"));
302         referenceRules.put("d4(f(g))/dt4",        referenceRules.get("d4(f(g))/dx4").replaceAll("x", "t"));
303         referenceRules.put("d4(f(g))/dx3dy",      "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dx * dg/dy +" +
304                                                   " 3 * d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dxdy +" +
305                                                   " 3 * d3(f(g))/dg3 * dg/dx * d2g/dx2 * dg/dy +" +
306                                                   " 3 * d2(f(g))/dg2 * d2g/dx2 * d2g/dxdy +" +
307                                                   " 3 * d2(f(g))/dg2 * dg/dx * d3g/dx2dy +" +
308                                                   " d2(f(g))/dg2 * d3g/dx3 * dg/dy +" +
309                                                   " d(f(g))/dg * d4g/dx3dy");
310         referenceRules.put("d4(f(g))/dx3dz",      referenceRules.get("d4(f(g))/dx3dy").replaceAll("y", "z"));
311         referenceRules.put("d4(f(g))/dx3dt",      referenceRules.get("d4(f(g))/dx3dy").replaceAll("y", "t"));
312         referenceRules.put("d4(f(g))/dxdy3",      "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dy * dg/dy +" +
313                                                   " 3 * d3(f(g))/dg3 * dg/dy * dg/dy * d2g/dxdy +" +
314                                                   " 3 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dy2 +" +
315                                                   " 3 * d2(f(g))/dg2 * d2g/dxdy * d2g/dy2 +" +
316                                                   " 3 * d2(f(g))/dg2 * dg/dy * d3g/dxdy2 +" +
317                                                   " d2(f(g))/dg2 * dg/dx * d3g/dy3 +" +
318                                                   " d(f(g))/dg * d4g/dxdy3");
319         referenceRules.put("d4(f(g))/dxdz3",      referenceRules.get("d4(f(g))/dxdy3").replaceAll("y", "z"));
320         referenceRules.put("d4(f(g))/dxdt3",      referenceRules.get("d4(f(g))/dxdy3").replaceAll("y", "t"));
321         referenceRules.put("d4(f(g))/dy3dz",      referenceRules.get("d4(f(g))/dx3dz").replaceAll("x", "y"));
322         referenceRules.put("d4(f(g))/dy3dt",      referenceRules.get("d4(f(g))/dx3dt").replaceAll("x", "y"));
323         referenceRules.put("d4(f(g))/dydz3",      referenceRules.get("d4(f(g))/dxdz3").replaceAll("x", "y"));
324         referenceRules.put("d4(f(g))/dydt3",      referenceRules.get("d4(f(g))/dxdt3").replaceAll("x", "y"));
325         referenceRules.put("d4(f(g))/dz3dt",      referenceRules.get("d4(f(g))/dx3dt").replaceAll("x", "z"));
326         referenceRules.put("d4(f(g))/dzdt3",      referenceRules.get("d4(f(g))/dxdt3").replaceAll("x", "z"));
327         referenceRules.put("d4(f(g))/dx2dy2",     "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dy * dg/dy +" +
328                                                   " 4 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dxdy +" +
329                                                   " d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dy2 +" +
330                                                   " 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dxdy +" +
331                                                   " 2 * d2(f(g))/dg2 * dg/dx * d3g/dxdy2 +" +
332                                                   " d3(f(g))/dg3 * d2g/dx2 * dg/dy * dg/dy +" +
333                                                   " 2 * d2(f(g))/dg2 * dg/dy * d3g/dx2dy +" +
334                                                   " d2(f(g))/dg2 * d2g/dx2 * d2g/dy2 +" +
335                                                   " d(f(g))/dg * d4g/dx2dy2");
336         referenceRules.put("d4(f(g))/dx2dz2",     referenceRules.get("d4(f(g))/dx2dy2").replaceAll("y", "z"));
337         referenceRules.put("d4(f(g))/dx2dt2",     referenceRules.get("d4(f(g))/dx2dy2").replaceAll("y", "t"));
338         referenceRules.put("d4(f(g))/dy2dz2",     referenceRules.get("d4(f(g))/dx2dz2").replaceAll("x", "y"));
339         referenceRules.put("d4(f(g))/dy2dt2",     referenceRules.get("d4(f(g))/dx2dt2").replaceAll("x", "y"));
340         referenceRules.put("d4(f(g))/dz2dt2",     referenceRules.get("d4(f(g))/dx2dt2").replaceAll("x", "z"));
341 
342         referenceRules.put("d4(f(g))/dx2dydz",    "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dy * dg/dz +" +
343                                                   " 2 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dxdz +" +
344                                                   " d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dydz +" +
345                                                   " 2 * d3(f(g))/dg3 * dg/dx * d2g/dxdy * dg/dz +" +
346                                                   " 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dxdz +" +
347                                                   " 2 * d2(f(g))/dg2 * dg/dx * d3g/dxdydz +" +
348                                                   " d3(f(g))/dg3 * d2g/dx2 * dg/dy * dg/dz +" +
349                                                   " d2(f(g))/dg2 * dg/dy * d3g/dx2dz +" +
350                                                   " d2(f(g))/dg2 * d2g/dx2 * d2g/dydz +" +
351                                                   " d2(f(g))/dg2 * d3g/dx2dy * dg/dz +" +
352                                                   " d(f(g))/dg * d4g/dx2dydz");
353         referenceRules.put("d4(f(g))/dx2dydt",    referenceRules.get("d4(f(g))/dx2dydz").replaceAll("z", "t"));
354         referenceRules.put("d4(f(g))/dx2dzdt",    referenceRules.get("d4(f(g))/dx2dydt").replaceAll("y", "z"));
355         referenceRules.put("d4(f(g))/dxdy2dz",    "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dy * dg/dz +" +
356                                                   " d3(f(g))/dg3 * dg/dy * dg/dy * d2g/dxdz +" +
357                                                   " 2 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dydz +" +
358                                                   " 2 * d3(f(g))/dg3 * dg/dy * d2g/dxdy * dg/dz +" +
359                                                   " 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dydz +" +
360                                                   " 2 * d2(f(g))/dg2 * dg/dy * d3g/dxdydz +" +
361                                                   " d3(f(g))/dg3 * dg/dx * d2g/dy2 * dg/dz +" +
362                                                   " d2(f(g))/dg2 * d2g/dy2 * d2g/dxdz +" +
363                                                   " d2(f(g))/dg2 * dg/dx * d3g/dy2dz +" +
364                                                   " d2(f(g))/dg2 * d3g/dxdy2 * dg/dz +" +
365                                                   " d(f(g))/dg * d4g/dxdy2dz");
366         referenceRules.put("d4(f(g))/dxdy2dt",    referenceRules.get("d4(f(g))/dxdy2dz").replaceAll("z", "t"));
367         referenceRules.put("d4(f(g))/dy2dzdt",    referenceRules.get("d4(f(g))/dx2dzdt").replaceAll("x", "y"));
368         referenceRules.put("d4(f(g))/dxdydz2",    "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dz * dg/dz +" +
369                                                   " 2 * d3(f(g))/dg3 * dg/dy * dg/dz * d2g/dxdz +" +
370                                                   " 2 * d3(f(g))/dg3 * dg/dx * dg/dz * d2g/dydz +" +
371                                                   " d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dz2 +" +
372                                                   " 2 * d2(f(g))/dg2 * d2g/dxdz * d2g/dydz +" +
373                                                   " d2(f(g))/dg2 * dg/dy * d3g/dxdz2 +" +
374                                                   " d2(f(g))/dg2 * dg/dx * d3g/dydz2 +" +
375                                                   " d3(f(g))/dg3 * d2g/dxdy * dg/dz * dg/dz +" +
376                                                   " 2 * d2(f(g))/dg2 * dg/dz * d3g/dxdydz +" +
377                                                   " d2(f(g))/dg2 * d2g/dxdy * d2g/dz2 +" +
378                                                   " d(f(g))/dg * d4g/dxdydz2");
379         referenceRules.put("d4(f(g))/dxdz2dt",    referenceRules.get("d4(f(g))/dxdy2dt").replaceAll("y", "z"));
380         referenceRules.put("d4(f(g))/dydz2dt",    referenceRules.get("d4(f(g))/dxdz2dt").replaceAll("x", "y"));
381         referenceRules.put("d4(f(g))/dxdydt2",    referenceRules.get("d4(f(g))/dxdydz2").replaceAll("z", "t"));
382         referenceRules.put("d4(f(g))/dxdzdt2",    referenceRules.get("d4(f(g))/dxdydt2").replaceAll("y", "z"));
383         referenceRules.put("d4(f(g))/dydzdt2",    referenceRules.get("d4(f(g))/dxdzdt2").replaceAll("x", "y"));
384         referenceRules.put("d4(f(g))/dxdydzdt",   "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dz * dg/dt +" +
385                                                   " d3(f(g))/dg3 * dg/dy * dg/dz * d2g/dxdt +" +
386                                                   " d3(f(g))/dg3 * dg/dx * dg/dz * d2g/dydt +" +
387                                                   " d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dzdt +" +
388                                                   " d3(f(g))/dg3 * dg/dy * d2g/dxdz * dg/dt +" +
389                                                   " d2(f(g))/dg2 * d2g/dxdz * d2g/dydt +" +
390                                                   " d2(f(g))/dg2 * dg/dy * d3g/dxdzdt +" +
391                                                   " d3(f(g))/dg3 * dg/dx * d2g/dydz * dg/dt +" +
392                                                   " d2(f(g))/dg2 * d2g/dydz * d2g/dxdt +" +
393                                                   " d2(f(g))/dg2 * dg/dx * d3g/dydzdt +" +
394                                                   " d3(f(g))/dg3 * d2g/dxdy * dg/dz * dg/dt +" +
395                                                   " d2(f(g))/dg2 * dg/dz * d3g/dxdydt +" +
396                                                   " d2(f(g))/dg2 * d2g/dxdy * d2g/dzdt +" +
397                                                   " d2(f(g))/dg2 * d3g/dxdydz * dg/dt +" +
398                                                   " d(f(g))/dg * d4g/dxdydzdt");
399 
400         Field compFieldArrayField = DSCompiler.class.getDeclaredField("compIndirection");
401         compFieldArrayField.setAccessible(true);
402         for (int i = 0; i < 5; ++i) {
403             for (int j = 0; j < 5; ++j) {
404                 DSCompiler compiler = DSCompiler.getCompiler(i, j);
405                 int[][][] compIndirection = (int[][][]) compFieldArrayField.get(compiler);
406                 for (int k = 0; k < compIndirection.length; ++k) {
407                     String product = ordersToString(compiler.getPartialDerivativeOrders(k),
408                                                     "(f(g))", "x", "y", "z", "t");
409                     StringBuilder rule = new StringBuilder();
410                     for (int[] term : compIndirection[k]) {
411                         if (rule.length() > 0) {
412                             rule.append(" + ");
413                         }
414                         if (term[0] > 1) {
415                             rule.append(term[0]).append(" * ");
416                         }
417                         rule.append(orderToString(term[1], "(f(g))", "g"));
418                         for (int l = 2; l < term.length; ++l) {
419                             rule.append(" * ");
420                             rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[l]),
421                                                        "g", "x", "y", "z", "t"));
422                         }
423                     }
424                     Assert.assertEquals(product, referenceRules.get(product), rule.toString());
425                 }
426             }
427         }
428     }
429 
430     private void checkIndices(int[] indices, int ... expected) {
431         Assert.assertEquals(expected.length, indices.length);
432         for (int i = 0; i < expected.length; ++i) {
433             Assert.assertEquals(expected[i], indices[i]);
434         }
435     }
436 
437     private String orderToString(int order, String functionName, String parameterName) {
438         if (order == 0) {
439             return functionName;
440         } else if (order == 1) {
441             return "d" + functionName + "/d" + parameterName;
442         } else {
443             return "d" + order + functionName + "/d" + parameterName + order;
444         }
445     }
446 
447     private String ordersToString(int[] orders, String functionName, String ... parametersNames) {
448 
449         int sumOrders = 0;
450         for (int order : orders) {
451             sumOrders += order;
452         }
453 
454         if (sumOrders == 0) {
455             return functionName;
456         }
457 
458         StringBuilder builder = new StringBuilder();
459         builder.append('d');
460         if (sumOrders > 1) {
461             builder.append(sumOrders);
462         }
463         builder.append(functionName).append('/');
464         for (int i = 0; i < orders.length; ++i) {
465             if (orders[i] > 0) {
466                 builder.append('d').append(parametersNames[i]);
467                 if (orders[i] > 1) {
468                     builder.append(orders[i]);
469                 }
470             }
471         }
472         return builder.toString();
473     }
474 }