001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.differentiation;
018
019 import java.util.ArrayList;
020 import java.util.Arrays;
021 import java.util.List;
022 import java.util.concurrent.atomic.AtomicReference;
023
024 import org.apache.commons.math3.exception.DimensionMismatchException;
025 import org.apache.commons.math3.exception.NumberIsTooLargeException;
026 import org.apache.commons.math3.util.ArithmeticUtils;
027 import org.apache.commons.math3.util.FastMath;
028 import org.apache.commons.math3.util.MathArrays;
029
030 /** Class holding "compiled" computation rules for derivative structures.
031 * <p>This class implements the computation rules described in Dan Kalman's paper <a
032 * href="http://www.math.american.edu/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
033 * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
034 * no. 3, June 2002. However, in order to avoid performances bottlenecks, the recursive
035 * rules are "compiled" once in an unfold form. This class does this recursion unrolling
036 * and stores the computation rules as simple loops with pre-computed indirection arrays.</p>
037 * <p>
038 * This class maps all derivative computation into single dimension arrays that hold the
039 * value and partial derivatives. The class does not hold these arrays, which remains under
040 * the responsibility of the caller. For each combination of number of free parameters and
041 * derivation order, only one compiler is necessary, and this compiler will be used to
042 * perform computations on all arrays provided to it, which can represent hundreds or
043 * thousands of different parameters kept together with all theur partial derivatives.
044 * </p>
045 * <p>
046 * The arrays on which compilers operate contain only the partial derivatives together
047 * with the 0<sup>th</sup> derivative, i.e. the value. The partial derivatives are stored in
048 * a compiler-specific order, which can be retrieved using methods {@link
049 * #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} and {@link
050 * #getPartialDerivativeOrders(int)}. The value is guaranteed to be stored as the first element
051 * (i.e. the {@link #getPartialDerivativeIndex(int...) getPartialDerivativeIndex} method returns
052 * 0 when called with 0 for all derivation orders and {@link #getPartialDerivativeOrders(int)
053 * getPartialDerivativeOrders} returns an array filled with 0 when called with 0 as the index).
054 * </p>
055 * <p>
056 * Note that the ordering changes with number of parameters and derivation order. For example
057 * given 2 parameters x and y, df/dy is stored at index 2 when derivation order is set to 1 (in
058 * this case the array has three elements: f, df/dx and df/dy). If derivation order is set to
059 * 2, then df/dy will be stored at index 3 (in this case the array has six elements: f, df/dx,
060 * df/dxdx, df/dy, df/dxdy and df/dydy).
061 * </p>
062 * <p>
063 * Given this structure, users can perform some simple operations like adding, subtracting
064 * or multiplying constants and negating the elements by themselves, knowing if they want to
065 * mutate their array or create a new array. These simple operations are not provided by
066 * the compiler. The compiler provides only the more complex operations between several arrays.
067 * </p>
068 * <p>This class is mainly used as the engine for scalar variable {@link DerivativeStructure}.
069 * It can also be used directly to hold several variables in arrays for more complex data
070 * structures. User can for example store a vector of n variables depending on three x, y
071 * and z free parameters in one array as follows:
072 * <pre>
073 * // parameter 0 is x, parameter 1 is y, parameter 2 is z
074 * int parameters = 3;
075 * DSCompiler compiler = DSCompiler.getCompiler(parameters, order);
076 * int size = compiler.getSize();
077 *
078 * // pack all elements in a single array
079 * double[] array = new double[n * size];
080 * for (int i = 0; i < n; ++i) {
081 *
082 * // we know value is guaranteed to be the first element
083 * array[i * size] = v[i];
084 *
085 * // we don't know where first derivatives are stored, so we ask the compiler
086 * array[i * size + compiler.getPartialDerivativeIndex(1, 0, 0) = dvOnDx[i][0];
087 * array[i * size + compiler.getPartialDerivativeIndex(0, 1, 0) = dvOnDy[i][0];
088 * array[i * size + compiler.getPartialDerivativeIndex(0, 0, 1) = dvOnDz[i][0];
089 *
090 * // we let all higher order derivatives set to 0
091 *
092 * }
093 * </pre>
094 * Then in another function, user can perform some operations on all elements stored
095 * in the single array, such as a simple product of all variables:
096 * <pre>
097 * // compute the product of all elements
098 * double[] product = new double[size];
099 * prod[0] = 1.0;
100 * for (int i = 0; i < n; ++i) {
101 * double[] tmp = product.clone();
102 * compiler.multiply(tmp, 0, array, i * size, product, 0);
103 * }
104 *
105 * // value
106 * double p = product[0];
107 *
108 * // first derivatives
109 * double dPdX = product[compiler.getPartialDerivativeIndex(1, 0, 0)];
110 * double dPdY = product[compiler.getPartialDerivativeIndex(0, 1, 0)];
111 * double dPdZ = product[compiler.getPartialDerivativeIndex(0, 0, 1)];
112 *
113 * // cross derivatives (assuming order was at least 2)
114 * double dPdXdX = product[compiler.getPartialDerivativeIndex(2, 0, 0)];
115 * double dPdXdY = product[compiler.getPartialDerivativeIndex(1, 1, 0)];
116 * double dPdXdZ = product[compiler.getPartialDerivativeIndex(1, 0, 1)];
117 * double dPdYdY = product[compiler.getPartialDerivativeIndex(0, 2, 0)];
118 * double dPdYdZ = product[compiler.getPartialDerivativeIndex(0, 1, 1)];
119 * double dPdZdZ = product[compiler.getPartialDerivativeIndex(0, 0, 2)];
120 * </p>
121 * @see DerivativeStructure
122 * @version $Id: DSCompiler.java 1421949 2012-12-14 15:53:32Z luc $
123 * @since 3.1
124 */
125 public class DSCompiler {
126
127 /** Array of all compilers created so far. */
128 private static AtomicReference<DSCompiler[][]> compilers =
129 new AtomicReference<DSCompiler[][]>(null);
130
131 /** Number of free parameters. */
132 private final int parameters;
133
134 /** Derivation order. */
135 private final int order;
136
137 /** Number of partial derivatives (including the single 0 order derivative element). */
138 private final int[][] sizes;
139
140 /** Indirection array for partial derivatives. */
141 private final int[][] derivativesIndirection;
142
143 /** Indirection array of the lower derivative elements. */
144 private final int[] lowerIndirection;
145
146 /** Indirection arrays for multiplication. */
147 private final int[][][] multIndirection;
148
149 /** Indirection arrays for function composition. */
150 private final int[][][] compIndirection;
151
152 /** Private constructor, reserved for the factory method {@link #getCompiler(int, int)}.
153 * @param parameters number of free parameters
154 * @param order derivation order
155 * @param valueCompiler compiler for the value part
156 * @param derivativeCompiler compiler for the derivative part
157 */
158 private DSCompiler(final int parameters, final int order,
159 final DSCompiler valueCompiler, final DSCompiler derivativeCompiler) {
160
161 this.parameters = parameters;
162 this.order = order;
163 this.sizes = compileSizes(parameters, order, valueCompiler);
164 this.derivativesIndirection =
165 compileDerivativesIndirection(parameters, order,
166 valueCompiler, derivativeCompiler);
167 this.lowerIndirection =
168 compileLowerIndirection(parameters, order,
169 valueCompiler, derivativeCompiler);
170 this.multIndirection =
171 compileMultiplicationIndirection(parameters, order,
172 valueCompiler, derivativeCompiler, lowerIndirection);
173 this.compIndirection =
174 compileCompositionIndirection(parameters, order,
175 valueCompiler, derivativeCompiler,
176 sizes, derivativesIndirection);
177
178 }
179
180 /** Get the compiler for number of free parameters and order.
181 * @param parameters number of free parameters
182 * @param order derivation order
183 * @return cached rules set
184 */
185 public static DSCompiler getCompiler(int parameters, int order) {
186
187 // get the cached compilers
188 final DSCompiler[][] cache = compilers.get();
189 if (cache != null && cache.length > parameters && cache[parameters].length > order) {
190 if (cache[parameters][order] != null) {
191 // the compiler has already been created
192 return cache[parameters][order];
193 }
194 }
195
196 // we need to create more compilers
197 final int maxParameters = FastMath.max(parameters, cache == null ? 0 : cache.length);
198 final int maxOrder = FastMath.max(order, cache == null ? 0 : cache[0].length);
199 final DSCompiler[][] newCache = new DSCompiler[maxParameters + 1][maxOrder + 1];
200
201 if (cache != null) {
202 // preserve the already created compilers
203 for (int i = 0; i < cache.length; ++i) {
204 System.arraycopy(cache[i], 0, newCache[i], 0, cache[i].length);
205 }
206 }
207
208 // create the array in increasing diagonal order
209 for (int diag = 0; diag <= parameters + order; ++diag) {
210 for (int o = FastMath.max(0, diag - parameters); o <= FastMath.min(order, diag); ++o) {
211 final int p = diag - o;
212 if (newCache[p][o] == null) {
213 final DSCompiler valueCompiler = (p == 0) ? null : newCache[p - 1][o];
214 final DSCompiler derivativeCompiler = (o == 0) ? null : newCache[p][o - 1];
215 newCache[p][o] = new DSCompiler(p, o, valueCompiler, derivativeCompiler);
216 }
217 }
218 }
219
220 // atomically reset the cached compilers array
221 compilers.compareAndSet(cache, newCache);
222
223 return newCache[parameters][order];
224
225 }
226
227 /** Compile the sizes array.
228 * @param parameters number of free parameters
229 * @param order derivation order
230 * @param valueCompiler compiler for the value part
231 * @return sizes array
232 */
233 private static int[][] compileSizes(final int parameters, final int order,
234 final DSCompiler valueCompiler) {
235
236 final int[][] sizes = new int[parameters + 1][order + 1];
237 if (parameters == 0) {
238 Arrays.fill(sizes[0], 1);
239 } else {
240 System.arraycopy(valueCompiler.sizes, 0, sizes, 0, parameters);
241 sizes[parameters][0] = 1;
242 for (int i = 0; i < order; ++i) {
243 sizes[parameters][i + 1] = sizes[parameters][i] + sizes[parameters - 1][i + 1];
244 }
245 }
246
247 return sizes;
248
249 }
250
251 /** Compile the derivatives indirection array.
252 * @param parameters number of free parameters
253 * @param order derivation order
254 * @param valueCompiler compiler for the value part
255 * @param derivativeCompiler compiler for the derivative part
256 * @return derivatives indirection array
257 */
258 private static int[][] compileDerivativesIndirection(final int parameters, final int order,
259 final DSCompiler valueCompiler,
260 final DSCompiler derivativeCompiler) {
261
262 if (parameters == 0 || order == 0) {
263 return new int[1][parameters];
264 }
265
266 final int vSize = valueCompiler.derivativesIndirection.length;
267 final int dSize = derivativeCompiler.derivativesIndirection.length;
268 final int[][] derivativesIndirection = new int[vSize + dSize][parameters];
269
270 // set up the indices for the value part
271 for (int i = 0; i < vSize; ++i) {
272 // copy the first indices, the last one remaining set to 0
273 System.arraycopy(valueCompiler.derivativesIndirection[i], 0,
274 derivativesIndirection[i], 0,
275 parameters - 1);
276 }
277
278 // set up the indices for the derivative part
279 for (int i = 0; i < dSize; ++i) {
280
281 // copy the indices
282 System.arraycopy(derivativeCompiler.derivativesIndirection[i], 0,
283 derivativesIndirection[vSize + i], 0,
284 parameters);
285
286 // increment the derivation order for the last parameter
287 derivativesIndirection[vSize + i][parameters - 1]++;
288
289 }
290
291 return derivativesIndirection;
292
293 }
294
295 /** Compile the lower derivatives indirection array.
296 * <p>
297 * This indirection array contains the indices of all elements
298 * except derivatives for last derivation order.
299 * </p>
300 * @param parameters number of free parameters
301 * @param order derivation order
302 * @param valueCompiler compiler for the value part
303 * @param derivativeCompiler compiler for the derivative part
304 * @return lower derivatives indirection array
305 */
306 private static int[] compileLowerIndirection(final int parameters, final int order,
307 final DSCompiler valueCompiler,
308 final DSCompiler derivativeCompiler) {
309
310 if (parameters == 0 || order <= 1) {
311 return new int[] { 0 };
312 }
313
314 // this is an implementation of definition 6 in Dan Kalman's paper.
315 final int vSize = valueCompiler.lowerIndirection.length;
316 final int dSize = derivativeCompiler.lowerIndirection.length;
317 final int[] lowerIndirection = new int[vSize + dSize];
318 System.arraycopy(valueCompiler.lowerIndirection, 0, lowerIndirection, 0, vSize);
319 for (int i = 0; i < dSize; ++i) {
320 lowerIndirection[vSize + i] = valueCompiler.getSize() + derivativeCompiler.lowerIndirection[i];
321 }
322
323 return lowerIndirection;
324
325 }
326
327 /** Compile the multiplication indirection array.
328 * <p>
329 * This indirection array contains the indices of all pairs of elements
330 * involved when computing a multiplication. This allows a straightforward
331 * loop-based multiplication (see {@link #multiply(double[], int, double[], int, double[], int)}).
332 * </p>
333 * @param parameters number of free parameters
334 * @param order derivation order
335 * @param valueCompiler compiler for the value part
336 * @param derivativeCompiler compiler for the derivative part
337 * @param lowerIndirection lower derivatives indirection array
338 * @return multiplication indirection array
339 */
340 private static int[][][] compileMultiplicationIndirection(final int parameters, final int order,
341 final DSCompiler valueCompiler,
342 final DSCompiler derivativeCompiler,
343 final int[] lowerIndirection) {
344
345 if ((parameters == 0) || (order == 0)) {
346 return new int[][][] { { { 1, 0, 0 } } };
347 }
348
349 // this is an implementation of definition 3 in Dan Kalman's paper.
350 final int vSize = valueCompiler.multIndirection.length;
351 final int dSize = derivativeCompiler.multIndirection.length;
352 final int[][][] multIndirection = new int[vSize + dSize][][];
353
354 System.arraycopy(valueCompiler.multIndirection, 0, multIndirection, 0, vSize);
355
356 for (int i = 0; i < dSize; ++i) {
357 final int[][] dRow = derivativeCompiler.multIndirection[i];
358 List<int[]> row = new ArrayList<int[]>();
359 for (int j = 0; j < dRow.length; ++j) {
360 row.add(new int[] { dRow[j][0], lowerIndirection[dRow[j][1]], vSize + dRow[j][2] });
361 row.add(new int[] { dRow[j][0], vSize + dRow[j][1], lowerIndirection[dRow[j][2]] });
362 }
363
364 // combine terms with similar derivation orders
365 final List<int[]> combined = new ArrayList<int[]>(row.size());
366 for (int j = 0; j < row.size(); ++j) {
367 final int[] termJ = row.get(j);
368 if (termJ[0] > 0) {
369 for (int k = j + 1; k < row.size(); ++k) {
370 final int[] termK = row.get(k);
371 if (termJ[1] == termK[1] && termJ[2] == termK[2]) {
372 // combine termJ and termK
373 termJ[0] += termK[0];
374 // make sure we will skip termK later on in the outer loop
375 termK[0] = 0;
376 }
377 }
378 combined.add(termJ);
379 }
380 }
381
382 multIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
383
384 }
385
386 return multIndirection;
387
388 }
389
390 /** Compile the function composition indirection array.
391 * <p>
392 * This indirection array contains the indices of all sets of elements
393 * involved when computing a composition. This allows a straightforward
394 * loop-based composition (see {@link #compose(double[], int, double[], double[], int)}).
395 * </p>
396 * @param parameters number of free parameters
397 * @param order derivation order
398 * @param valueCompiler compiler for the value part
399 * @param derivativeCompiler compiler for the derivative part
400 * @param sizes sizes array
401 * @param derivativesIndirection derivatives indirection array
402 * @return multiplication indirection array
403 */
404 private static int[][][] compileCompositionIndirection(final int parameters, final int order,
405 final DSCompiler valueCompiler,
406 final DSCompiler derivativeCompiler,
407 final int[][] sizes,
408 final int[][] derivativesIndirection) {
409
410 if ((parameters == 0) || (order == 0)) {
411 return new int[][][] { { { 1, 0 } } };
412 }
413
414 final int vSize = valueCompiler.compIndirection.length;
415 final int dSize = derivativeCompiler.compIndirection.length;
416 final int[][][] compIndirection = new int[vSize + dSize][][];
417
418 // the composition rules from the value part can be reused as is
419 System.arraycopy(valueCompiler.compIndirection, 0, compIndirection, 0, vSize);
420
421 // the composition rules for the derivative part are deduced by
422 // differentiation the rules from the underlying compiler once
423 // with respect to the parameter this compiler handles and the
424 // underlying one did not handle
425 for (int i = 0; i < dSize; ++i) {
426 List<int[]> row = new ArrayList<int[]>();
427 for (int[] term : derivativeCompiler.compIndirection[i]) {
428
429 // handle term p * f_k(g(x)) * g_l1(x) * g_l2(x) * ... * g_lp(x)
430
431 // derive the first factor in the term: f_k with respect to new parameter
432 int[] derivedTermF = new int[term.length + 1];
433 derivedTermF[0] = term[0]; // p
434 derivedTermF[1] = term[1] + 1; // f_(k+1)
435 int[] orders = new int[parameters];
436 orders[parameters - 1] = 1;
437 derivedTermF[term.length] = getPartialDerivativeIndex(parameters, order, sizes, orders); // g_1
438 for (int j = 2; j < term.length; ++j) {
439 // convert the indices as the mapping for the current order
440 // is different from the mapping with one less order
441 derivedTermF[j] = convertIndex(term[j], parameters,
442 derivativeCompiler.derivativesIndirection,
443 parameters, order, sizes);
444 }
445 Arrays.sort(derivedTermF, 2, derivedTermF.length);
446 row.add(derivedTermF);
447
448 // derive the various g_l
449 for (int l = 2; l < term.length; ++l) {
450 int[] derivedTermG = new int[term.length];
451 derivedTermG[0] = term[0];
452 derivedTermG[1] = term[1];
453 for (int j = 2; j < term.length; ++j) {
454 // convert the indices as the mapping for the current order
455 // is different from the mapping with one less order
456 derivedTermG[j] = convertIndex(term[j], parameters,
457 derivativeCompiler.derivativesIndirection,
458 parameters, order, sizes);
459 if (j == l) {
460 // derive this term
461 System.arraycopy(derivativesIndirection[derivedTermG[j]], 0, orders, 0, parameters);
462 orders[parameters - 1]++;
463 derivedTermG[j] = getPartialDerivativeIndex(parameters, order, sizes, orders);
464 }
465 }
466 Arrays.sort(derivedTermG, 2, derivedTermG.length);
467 row.add(derivedTermG);
468 }
469
470 }
471
472 // combine terms with similar derivation orders
473 final List<int[]> combined = new ArrayList<int[]>(row.size());
474 for (int j = 0; j < row.size(); ++j) {
475 final int[] termJ = row.get(j);
476 if (termJ[0] > 0) {
477 for (int k = j + 1; k < row.size(); ++k) {
478 final int[] termK = row.get(k);
479 boolean equals = termJ.length == termK.length;
480 for (int l = 1; equals && l < termJ.length; ++l) {
481 equals &= termJ[l] == termK[l];
482 }
483 if (equals) {
484 // combine termJ and termK
485 termJ[0] += termK[0];
486 // make sure we will skip termK later on in the outer loop
487 termK[0] = 0;
488 }
489 }
490 combined.add(termJ);
491 }
492 }
493
494 compIndirection[vSize + i] = combined.toArray(new int[combined.size()][]);
495
496 }
497
498 return compIndirection;
499
500 }
501
502 /** Get the index of a partial derivative in the array.
503 * <p>
504 * If all orders are set to 0, then the 0<sup>th</sup> order derivative
505 * is returned, which is the value of the function.
506 * </p>
507 * <p>The indices of derivatives are between 0 and {@link #getSize() getSize()} - 1.
508 * Their specific order is fixed for a given compiler, but otherwise not
509 * publicly specified. There are however some simple cases which have guaranteed
510 * indices:
511 * </p>
512 * <ul>
513 * <li>the index of 0<sup>th</sup> order derivative is always 0</li>
514 * <li>if there is only 1 {@link #getFreeParameters() free parameter}, then the
515 * derivatives are sorted in increasing derivation order (i.e. f at index 0, df/dp
516 * at index 1, d<sup>2</sup>f/dp<sup>2</sup> at index 2 ...
517 * d<sup>k</sup>f/dp<sup>k</sup> at index k),</li>
518 * <li>if the {@link #getOrder() derivation order} is 1, then the derivatives
519 * are sorted in incresing free parameter order (i.e. f at index 0, df/dx<sub>1</sub>
520 * at index 1, df/dx<sub>2</sub> at index 2 ... df/dx<sub>k</sub> at index k),</li>
521 * <li>all other cases are not publicly specified</li>
522 * </ul>
523 * <p>
524 * This method is the inverse of method {@link #getPartialDerivativeOrders(int)}
525 * </p>
526 * @param orders derivation orders with respect to each parameter
527 * @return index of the partial derivative
528 * @exception DimensionMismatchException if the numbers of parameters does not
529 * match the instance
530 * @exception NumberIsTooLargeException if sum of derivation orders is larger
531 * than the instance limits
532 * @see #getPartialDerivativeOrders(int)
533 */
534 public int getPartialDerivativeIndex(final int ... orders)
535 throws DimensionMismatchException, NumberIsTooLargeException {
536
537 // safety check
538 if (orders.length != getFreeParameters()) {
539 throw new DimensionMismatchException(orders.length, getFreeParameters());
540 }
541
542 return getPartialDerivativeIndex(parameters, order, sizes, orders);
543
544 }
545
546 /** Get the index of a partial derivative in an array.
547 * @param parameters number of free parameters
548 * @param order derivation order
549 * @param sizes sizes array
550 * @param orders derivation orders with respect to each parameter
551 * (the lenght of this array must match the number of parameters)
552 * @return index of the partial derivative
553 * @exception NumberIsTooLargeException if sum of derivation orders is larger
554 * than the instance limits
555 */
556 private static int getPartialDerivativeIndex(final int parameters, final int order,
557 final int[][] sizes, final int ... orders)
558 throws NumberIsTooLargeException {
559
560 // the value is obtained by diving into the recursive Dan Kalman's structure
561 // this is theorem 2 of his paper, with recursion replaced by iteration
562 int index = 0;
563 int m = order;
564 int ordersSum = 0;
565 for (int i = parameters - 1; i >= 0; --i) {
566
567 // derivative order for current free parameter
568 int derivativeOrder = orders[i];
569
570 // safety check
571 ordersSum += derivativeOrder;
572 if (ordersSum > order) {
573 throw new NumberIsTooLargeException(ordersSum, order, true);
574 }
575
576 while (derivativeOrder-- > 0) {
577 // as long as we differentiate according to current free parameter,
578 // we have to skip the value part and dive into the derivative part
579 // so we add the size of the value part to the base index
580 index += sizes[i][m--];
581 }
582
583 }
584
585 return index;
586
587 }
588
589 /** Convert an index from one (parameters, order) structure to another.
590 * @param index index of a partial derivative in source derivative structure
591 * @param srcP number of free parameters in source derivative structure
592 * @param srcDerivativesIndirection derivatives indirection array for the source
593 * derivative structure
594 * @param destP number of free parameters in destination derivative structure
595 * @param destO derivation order in destination derivative structure
596 * @param destSizes sizes array for the destination derivative structure
597 * @return index of the partial derivative with the <em>same</em> characteristics
598 * in destination derivative structure
599 */
600 private static int convertIndex(final int index,
601 final int srcP, final int[][] srcDerivativesIndirection,
602 final int destP, final int destO, final int[][] destSizes) {
603 int[] orders = new int[destP];
604 System.arraycopy(srcDerivativesIndirection[index], 0, orders, 0, FastMath.min(srcP, destP));
605 return getPartialDerivativeIndex(destP, destO, destSizes, orders);
606 }
607
608 /** Get the derivation orders for a specific index in the array.
609 * <p>
610 * This method is the inverse of {@link #getPartialDerivativeIndex(int...)}.
611 * </p>
612 * @param index of the partial derivative
613 * @return orders derivation orders with respect to each parameter
614 * @see #getPartialDerivativeIndex(int...)
615 */
616 public int[] getPartialDerivativeOrders(final int index) {
617 return derivativesIndirection[index];
618 }
619
620 /** Get the number of free parameters.
621 * @return number of free parameters
622 */
623 public int getFreeParameters() {
624 return parameters;
625 }
626
627 /** Get the derivation order.
628 * @return derivation order
629 */
630 public int getOrder() {
631 return order;
632 }
633
634 /** Get the array size required for holding partial derivatives data.
635 * <p>
636 * This number includes the single 0 order derivative element, which is
637 * guaranteed to be stored in the first element of the array.
638 * </p>
639 * @return array size required for holding partial derivatives data
640 */
641 public int getSize() {
642 return sizes[parameters][order];
643 }
644
645 /** Compute linear combination.
646 * The derivative structure built will be a1 * ds1 + a2 * ds2
647 * @param a1 first scale factor
648 * @param c1 first base (unscaled) component
649 * @param offset1 offset of first operand in its array
650 * @param a2 second scale factor
651 * @param c2 second base (unscaled) component
652 * @param offset2 offset of second operand in its array
653 * @param result array where result must be stored (it may be
654 * one of the input arrays)
655 * @param resultOffset offset of the result in its array
656 */
657 public void linearCombination(final double a1, final double[] c1, final int offset1,
658 final double a2, final double[] c2, final int offset2,
659 final double[] result, final int resultOffset) {
660 for (int i = 0; i < getSize(); ++i) {
661 result[resultOffset + i] =
662 MathArrays.linearCombination(a1, c1[offset1 + i], a2, c2[offset2 + i]);
663 }
664 }
665
666 /** Compute linear combination.
667 * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
668 * @param a1 first scale factor
669 * @param c1 first base (unscaled) component
670 * @param offset1 offset of first operand in its array
671 * @param a2 second scale factor
672 * @param c2 second base (unscaled) component
673 * @param offset2 offset of second operand in its array
674 * @param a3 third scale factor
675 * @param c3 third base (unscaled) component
676 * @param offset3 offset of third operand in its array
677 * @param result array where result must be stored (it may be
678 * one of the input arrays)
679 * @param resultOffset offset of the result in its array
680 */
681 public void linearCombination(final double a1, final double[] c1, final int offset1,
682 final double a2, final double[] c2, final int offset2,
683 final double a3, final double[] c3, final int offset3,
684 final double[] result, final int resultOffset) {
685 for (int i = 0; i < getSize(); ++i) {
686 result[resultOffset + i] =
687 MathArrays.linearCombination(a1, c1[offset1 + i],
688 a2, c2[offset2 + i],
689 a3, c3[offset3 + i]);
690 }
691 }
692
693 /** Compute linear combination.
694 * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
695 * @param a1 first scale factor
696 * @param c1 first base (unscaled) component
697 * @param offset1 offset of first operand in its array
698 * @param a2 second scale factor
699 * @param c2 second base (unscaled) component
700 * @param offset2 offset of second operand in its array
701 * @param a3 third scale factor
702 * @param c3 third base (unscaled) component
703 * @param offset3 offset of third operand in its array
704 * @param a4 fourth scale factor
705 * @param c4 fourth base (unscaled) component
706 * @param offset4 offset of fourth operand in its array
707 * @param result array where result must be stored (it may be
708 * one of the input arrays)
709 * @param resultOffset offset of the result in its array
710 */
711 public void linearCombination(final double a1, final double[] c1, final int offset1,
712 final double a2, final double[] c2, final int offset2,
713 final double a3, final double[] c3, final int offset3,
714 final double a4, final double[] c4, final int offset4,
715 final double[] result, final int resultOffset) {
716 for (int i = 0; i < getSize(); ++i) {
717 result[resultOffset + i] =
718 MathArrays.linearCombination(a1, c1[offset1 + i],
719 a2, c2[offset2 + i],
720 a3, c3[offset3 + i],
721 a4, c4[offset4 + i]);
722 }
723 }
724
725 /** Perform addition of two derivative structures.
726 * @param lhs array holding left hand side of addition
727 * @param lhsOffset offset of the left hand side in its array
728 * @param rhs array right hand side of addition
729 * @param rhsOffset offset of the right hand side in its array
730 * @param result array where result must be stored (it may be
731 * one of the input arrays)
732 * @param resultOffset offset of the result in its array
733 */
734 public void add(final double[] lhs, final int lhsOffset,
735 final double[] rhs, final int rhsOffset,
736 final double[] result, final int resultOffset) {
737 for (int i = 0; i < getSize(); ++i) {
738 result[resultOffset + i] = lhs[lhsOffset + i] + rhs[rhsOffset + i];
739 }
740 }
741 /** Perform subtraction of two derivative structures.
742 * @param lhs array holding left hand side of subtraction
743 * @param lhsOffset offset of the left hand side in its array
744 * @param rhs array right hand side of subtraction
745 * @param rhsOffset offset of the right hand side in its array
746 * @param result array where result must be stored (it may be
747 * one of the input arrays)
748 * @param resultOffset offset of the result in its array
749 */
750 public void subtract(final double[] lhs, final int lhsOffset,
751 final double[] rhs, final int rhsOffset,
752 final double[] result, final int resultOffset) {
753 for (int i = 0; i < getSize(); ++i) {
754 result[resultOffset + i] = lhs[lhsOffset + i] - rhs[rhsOffset + i];
755 }
756 }
757
758 /** Perform multiplication of two derivative structures.
759 * @param lhs array holding left hand side of multiplication
760 * @param lhsOffset offset of the left hand side in its array
761 * @param rhs array right hand side of multiplication
762 * @param rhsOffset offset of the right hand side in its array
763 * @param result array where result must be stored (for
764 * multiplication the result array <em>cannot</em> be one of
765 * the input arrays)
766 * @param resultOffset offset of the result in its array
767 */
768 public void multiply(final double[] lhs, final int lhsOffset,
769 final double[] rhs, final int rhsOffset,
770 final double[] result, final int resultOffset) {
771 for (int i = 0; i < multIndirection.length; ++i) {
772 final int[][] mappingI = multIndirection[i];
773 double r = 0;
774 for (int j = 0; j < mappingI.length; ++j) {
775 r += mappingI[j][0] *
776 lhs[lhsOffset + mappingI[j][1]] *
777 rhs[rhsOffset + mappingI[j][2]];
778 }
779 result[resultOffset + i] = r;
780 }
781 }
782
783 /** Perform division of two derivative structures.
784 * @param lhs array holding left hand side of division
785 * @param lhsOffset offset of the left hand side in its array
786 * @param rhs array right hand side of division
787 * @param rhsOffset offset of the right hand side in its array
788 * @param result array where result must be stored (for
789 * division the result array <em>cannot</em> be one of
790 * the input arrays)
791 * @param resultOffset offset of the result in its array
792 */
793 public void divide(final double[] lhs, final int lhsOffset,
794 final double[] rhs, final int rhsOffset,
795 final double[] result, final int resultOffset) {
796 final double[] reciprocal = new double[getSize()];
797 pow(rhs, lhsOffset, -1, reciprocal, 0);
798 multiply(lhs, lhsOffset, reciprocal, 0, result, resultOffset);
799 }
800
801 /** Perform remainder of two derivative structures.
802 * @param lhs array holding left hand side of remainder
803 * @param lhsOffset offset of the left hand side in its array
804 * @param rhs array right hand side of remainder
805 * @param rhsOffset offset of the right hand side in its array
806 * @param result array where result must be stored (it may be
807 * one of the input arrays)
808 * @param resultOffset offset of the result in its array
809 */
810 public void remainder(final double[] lhs, final int lhsOffset,
811 final double[] rhs, final int rhsOffset,
812 final double[] result, final int resultOffset) {
813
814 // compute k such that lhs % rhs = lhs - k rhs
815 final double rem = lhs[lhsOffset] % rhs[rhsOffset];
816 final double k = FastMath.rint((lhs[lhsOffset] - rem) / rhs[rhsOffset]);
817
818 // set up value
819 result[resultOffset] = rem;
820
821 // set up partial derivatives
822 for (int i = 1; i < getSize(); ++i) {
823 result[resultOffset + i] = lhs[lhsOffset + i] - k * rhs[rhsOffset + i];
824 }
825
826 }
827
828 /** Compute power of a derivative structure.
829 * @param operand array holding the operand
830 * @param operandOffset offset of the operand in its array
831 * @param p power to apply
832 * @param result array where result must be stored (for
833 * power the result array <em>cannot</em> be the input
834 * array)
835 * @param resultOffset offset of the result in its array
836 */
837 public void pow(final double[] operand, final int operandOffset, final double p,
838 final double[] result, final int resultOffset) {
839
840 // create the function value and derivatives
841 // [x^p, px^(p-1), p(p-1)x^(p-2), ... ]
842 double[] function = new double[1 + order];
843 double xk = FastMath.pow(operand[operandOffset], p - order);
844 for (int i = order; i > 0; --i) {
845 function[i] = xk;
846 xk *= operand[operandOffset];
847 }
848 function[0] = xk;
849 double coefficient = p;
850 for (int i = 1; i <= order; ++i) {
851 function[i] *= coefficient;
852 coefficient *= p - i;
853 }
854
855 // apply function composition
856 compose(operand, operandOffset, function, result, resultOffset);
857
858 }
859
860 /** Compute integer power of a derivative structure.
861 * @param operand array holding the operand
862 * @param operandOffset offset of the operand in its array
863 * @param n power to apply
864 * @param result array where result must be stored (for
865 * power the result array <em>cannot</em> be the input
866 * array)
867 * @param resultOffset offset of the result in its array
868 */
869 public void pow(final double[] operand, final int operandOffset, final int n,
870 final double[] result, final int resultOffset) {
871
872 if (n == 0) {
873 // special case, x^0 = 1 for all x
874 result[resultOffset] = 1.0;
875 Arrays.fill(result, resultOffset + 1, resultOffset + getSize(), 0);
876 return;
877 }
878
879 // create the power function value and derivatives
880 // [x^n, nx^(n-1), n(n-1)x^(n-2), ... ]
881 double[] function = new double[1 + order];
882
883 if (n > 0) {
884 // strictly positive power
885 final int maxOrder = FastMath.min(order, n);
886 double xk = FastMath.pow(operand[operandOffset], n - maxOrder);
887 for (int i = maxOrder; i > 0; --i) {
888 function[i] = xk;
889 xk *= operand[operandOffset];
890 }
891 function[0] = xk;
892 } else {
893 // strictly negative power
894 final double inv = 1.0 / operand[operandOffset];
895 double xk = FastMath.pow(inv, -n);
896 for (int i = 0; i <= order; ++i) {
897 function[i] = xk;
898 xk *= inv;
899 }
900 }
901
902 double coefficient = n;
903 for (int i = 1; i <= order; ++i) {
904 function[i] *= coefficient;
905 coefficient *= n - i;
906 }
907
908 // apply function composition
909 compose(operand, operandOffset, function, result, resultOffset);
910
911 }
912
913 /** Compute power of a derivative structure.
914 * @param x array holding the base
915 * @param xOffset offset of the base in its array
916 * @param y array holding the exponent
917 * @param yOffset offset of the exponent in its array
918 * @param result array where result must be stored (for
919 * power the result array <em>cannot</em> be the input
920 * array)
921 * @param resultOffset offset of the result in its array
922 */
923 public void pow(final double[] x, final int xOffset,
924 final double[] y, final int yOffset,
925 final double[] result, final int resultOffset) {
926 final double[] logX = new double[getSize()];
927 log(x, xOffset, logX, 0);
928 final double[] yLogX = new double[getSize()];
929 multiply(logX, 0, y, yOffset, yLogX, 0);
930 exp(yLogX, 0, result, resultOffset);
931 }
932
933 /** Compute n<sup>th</sup> root of a derivative structure.
934 * @param operand array holding the operand
935 * @param operandOffset offset of the operand in its array
936 * @param n order of the root
937 * @param result array where result must be stored (for
938 * n<sup>th</sup> root the result array <em>cannot</em> be the input
939 * array)
940 * @param resultOffset offset of the result in its array
941 */
942 public void rootN(final double[] operand, final int operandOffset, final int n,
943 final double[] result, final int resultOffset) {
944
945 // create the function value and derivatives
946 // [x^(1/n), (1/n)x^((1/n)-1), (1-n)/n^2x^((1/n)-2), ... ]
947 double[] function = new double[1 + order];
948 double xk;
949 if (n == 2) {
950 function[0] = FastMath.sqrt(operand[operandOffset]);
951 xk = 0.5 / function[0];
952 } else if (n == 3) {
953 function[0] = FastMath.cbrt(operand[operandOffset]);
954 xk = 1.0 / (3.0 * function[0] * function[0]);
955 } else {
956 function[0] = FastMath.pow(operand[operandOffset], 1.0 / n);
957 xk = 1.0 / (n * FastMath.pow(function[0], n - 1));
958 }
959 final double nReciprocal = 1.0 / n;
960 final double xReciprocal = 1.0 / operand[operandOffset];
961 for (int i = 1; i <= order; ++i) {
962 function[i] = xk;
963 xk *= xReciprocal * (nReciprocal - i);
964 }
965
966 // apply function composition
967 compose(operand, operandOffset, function, result, resultOffset);
968
969 }
970
971 /** Compute exponential of a derivative structure.
972 * @param operand array holding the operand
973 * @param operandOffset offset of the operand in its array
974 * @param result array where result must be stored (for
975 * exponential the result array <em>cannot</em> be the input
976 * array)
977 * @param resultOffset offset of the result in its array
978 */
979 public void exp(final double[] operand, final int operandOffset,
980 final double[] result, final int resultOffset) {
981
982 // create the function value and derivatives
983 double[] function = new double[1 + order];
984 Arrays.fill(function, FastMath.exp(operand[operandOffset]));
985
986 // apply function composition
987 compose(operand, operandOffset, function, result, resultOffset);
988
989 }
990
991 /** Compute exp(x) - 1 of a derivative structure.
992 * @param operand array holding the operand
993 * @param operandOffset offset of the operand in its array
994 * @param result array where result must be stored (for
995 * exponential the result array <em>cannot</em> be the input
996 * array)
997 * @param resultOffset offset of the result in its array
998 */
999 public void expm1(final double[] operand, final int operandOffset,
1000 final double[] result, final int resultOffset) {
1001
1002 // create the function value and derivatives
1003 double[] function = new double[1 + order];
1004 function[0] = FastMath.expm1(operand[operandOffset]);
1005 Arrays.fill(function, 1, 1 + order, FastMath.exp(operand[operandOffset]));
1006
1007 // apply function composition
1008 compose(operand, operandOffset, function, result, resultOffset);
1009
1010 }
1011
1012 /** Compute natural logarithm of a derivative structure.
1013 * @param operand array holding the operand
1014 * @param operandOffset offset of the operand in its array
1015 * @param result array where result must be stored (for
1016 * logarithm the result array <em>cannot</em> be the input
1017 * array)
1018 * @param resultOffset offset of the result in its array
1019 */
1020 public void log(final double[] operand, final int operandOffset,
1021 final double[] result, final int resultOffset) {
1022
1023 // create the function value and derivatives
1024 double[] function = new double[1 + order];
1025 function[0] = FastMath.log(operand[operandOffset]);
1026 if (order > 0) {
1027 double inv = 1.0 / operand[operandOffset];
1028 double xk = inv;
1029 for (int i = 1; i <= order; ++i) {
1030 function[i] = xk;
1031 xk *= -i * inv;
1032 }
1033 }
1034
1035 // apply function composition
1036 compose(operand, operandOffset, function, result, resultOffset);
1037
1038 }
1039
1040 /** Computes shifted logarithm of a derivative structure.
1041 * @param operand array holding the operand
1042 * @param operandOffset offset of the operand in its array
1043 * @param result array where result must be stored (for
1044 * shifted logarithm the result array <em>cannot</em> be the input array)
1045 * @param resultOffset offset of the result in its array
1046 */
1047 public void log1p(final double[] operand, final int operandOffset,
1048 final double[] result, final int resultOffset) {
1049
1050 // create the function value and derivatives
1051 double[] function = new double[1 + order];
1052 function[0] = FastMath.log1p(operand[operandOffset]);
1053 if (order > 0) {
1054 double inv = 1.0 / (1.0 + operand[operandOffset]);
1055 double xk = inv;
1056 for (int i = 1; i <= order; ++i) {
1057 function[i] = xk;
1058 xk *= -i * inv;
1059 }
1060 }
1061
1062 // apply function composition
1063 compose(operand, operandOffset, function, result, resultOffset);
1064
1065 }
1066
1067 /** Computes base 10 logarithm of a derivative structure.
1068 * @param operand array holding the operand
1069 * @param operandOffset offset of the operand in its array
1070 * @param result array where result must be stored (for
1071 * base 10 logarithm the result array <em>cannot</em> be the input array)
1072 * @param resultOffset offset of the result in its array
1073 */
1074 public void log10(final double[] operand, final int operandOffset,
1075 final double[] result, final int resultOffset) {
1076
1077 // create the function value and derivatives
1078 double[] function = new double[1 + order];
1079 function[0] = FastMath.log10(operand[operandOffset]);
1080 if (order > 0) {
1081 double inv = 1.0 / operand[operandOffset];
1082 double xk = inv / FastMath.log(10.0);
1083 for (int i = 1; i <= order; ++i) {
1084 function[i] = xk;
1085 xk *= -i * inv;
1086 }
1087 }
1088
1089 // apply function composition
1090 compose(operand, operandOffset, function, result, resultOffset);
1091
1092 }
1093
1094 /** Compute cosine of a derivative structure.
1095 * @param operand array holding the operand
1096 * @param operandOffset offset of the operand in its array
1097 * @param result array where result must be stored (for
1098 * cosine the result array <em>cannot</em> be the input
1099 * array)
1100 * @param resultOffset offset of the result in its array
1101 */
1102 public void cos(final double[] operand, final int operandOffset,
1103 final double[] result, final int resultOffset) {
1104
1105 // create the function value and derivatives
1106 double[] function = new double[1 + order];
1107 function[0] = FastMath.cos(operand[operandOffset]);
1108 if (order > 0) {
1109 function[1] = -FastMath.sin(operand[operandOffset]);
1110 for (int i = 2; i <= order; ++i) {
1111 function[i] = -function[i - 2];
1112 }
1113 }
1114
1115 // apply function composition
1116 compose(operand, operandOffset, function, result, resultOffset);
1117
1118 }
1119
1120 /** Compute sine of a derivative structure.
1121 * @param operand array holding the operand
1122 * @param operandOffset offset of the operand in its array
1123 * @param result array where result must be stored (for
1124 * sine the result array <em>cannot</em> be the input
1125 * array)
1126 * @param resultOffset offset of the result in its array
1127 */
1128 public void sin(final double[] operand, final int operandOffset,
1129 final double[] result, final int resultOffset) {
1130
1131 // create the function value and derivatives
1132 double[] function = new double[1 + order];
1133 function[0] = FastMath.sin(operand[operandOffset]);
1134 if (order > 0) {
1135 function[1] = FastMath.cos(operand[operandOffset]);
1136 for (int i = 2; i <= order; ++i) {
1137 function[i] = -function[i - 2];
1138 }
1139 }
1140
1141 // apply function composition
1142 compose(operand, operandOffset, function, result, resultOffset);
1143
1144 }
1145
1146 /** Compute tangent of a derivative structure.
1147 * @param operand array holding the operand
1148 * @param operandOffset offset of the operand in its array
1149 * @param result array where result must be stored (for
1150 * tangent the result array <em>cannot</em> be the input
1151 * array)
1152 * @param resultOffset offset of the result in its array
1153 */
1154 public void tan(final double[] operand, final int operandOffset,
1155 final double[] result, final int resultOffset) {
1156
1157 // create the function value and derivatives
1158 final double[] function = new double[1 + order];
1159 final double t = FastMath.tan(operand[operandOffset]);
1160 function[0] = t;
1161
1162 if (order > 0) {
1163
1164 // the nth order derivative of tan has the form:
1165 // dn(tan(x)/dxn = P_n(tan(x))
1166 // where P_n(t) is a degree n+1 polynomial with same parity as n+1
1167 // P_0(t) = t, P_1(t) = 1 + t^2, P_2(t) = 2 t (1 + t^2) ...
1168 // the general recurrence relation for P_n is:
1169 // P_n(x) = (1+t^2) P_(n-1)'(t)
1170 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1171 final double[] p = new double[order + 2];
1172 p[1] = 1;
1173 final double t2 = t * t;
1174 for (int n = 1; n <= order; ++n) {
1175
1176 // update and evaluate polynomial P_n(t)
1177 double v = 0;
1178 p[n + 1] = n * p[n];
1179 for (int k = n + 1; k >= 0; k -= 2) {
1180 v = v * t2 + p[k];
1181 if (k > 2) {
1182 p[k - 2] = (k - 1) * p[k - 1] + (k - 3) * p[k - 3];
1183 } else if (k == 2) {
1184 p[0] = p[1];
1185 }
1186 }
1187 if ((n & 0x1) == 0) {
1188 v *= t;
1189 }
1190
1191 function[n] = v;
1192
1193 }
1194 }
1195
1196 // apply function composition
1197 compose(operand, operandOffset, function, result, resultOffset);
1198
1199 }
1200
1201 /** Compute arc cosine of a derivative structure.
1202 * @param operand array holding the operand
1203 * @param operandOffset offset of the operand in its array
1204 * @param result array where result must be stored (for
1205 * arc cosine the result array <em>cannot</em> be the input
1206 * array)
1207 * @param resultOffset offset of the result in its array
1208 */
1209 public void acos(final double[] operand, final int operandOffset,
1210 final double[] result, final int resultOffset) {
1211
1212 // create the function value and derivatives
1213 double[] function = new double[1 + order];
1214 final double x = operand[operandOffset];
1215 function[0] = FastMath.acos(x);
1216 if (order > 0) {
1217 // the nth order derivative of acos has the form:
1218 // dn(acos(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
1219 // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1220 // P_1(x) = -1, P_2(x) = -x, P_3(x) = -2x^2 - 1 ...
1221 // the general recurrence relation for P_n is:
1222 // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
1223 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1224 final double[] p = new double[order];
1225 p[0] = -1;
1226 final double x2 = x * x;
1227 final double f = 1.0 / (1 - x2);
1228 double coeff = FastMath.sqrt(f);
1229 function[1] = coeff * p[0];
1230 for (int n = 2; n <= order; ++n) {
1231
1232 // update and evaluate polynomial P_n(x)
1233 double v = 0;
1234 p[n - 1] = (n - 1) * p[n - 2];
1235 for (int k = n - 1; k >= 0; k -= 2) {
1236 v = v * x2 + p[k];
1237 if (k > 2) {
1238 p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
1239 } else if (k == 2) {
1240 p[0] = p[1];
1241 }
1242 }
1243 if ((n & 0x1) == 0) {
1244 v *= x;
1245 }
1246
1247 coeff *= f;
1248 function[n] = coeff * v;
1249
1250 }
1251 }
1252
1253 // apply function composition
1254 compose(operand, operandOffset, function, result, resultOffset);
1255
1256 }
1257
1258 /** Compute arc sine of a derivative structure.
1259 * @param operand array holding the operand
1260 * @param operandOffset offset of the operand in its array
1261 * @param result array where result must be stored (for
1262 * arc sine the result array <em>cannot</em> be the input
1263 * array)
1264 * @param resultOffset offset of the result in its array
1265 */
1266 public void asin(final double[] operand, final int operandOffset,
1267 final double[] result, final int resultOffset) {
1268
1269 // create the function value and derivatives
1270 double[] function = new double[1 + order];
1271 final double x = operand[operandOffset];
1272 function[0] = FastMath.asin(x);
1273 if (order > 0) {
1274 // the nth order derivative of asin has the form:
1275 // dn(asin(x)/dxn = P_n(x) / [1 - x^2]^((2n-1)/2)
1276 // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1277 // P_1(x) = 1, P_2(x) = x, P_3(x) = 2x^2 + 1 ...
1278 // the general recurrence relation for P_n is:
1279 // P_n(x) = (1-x^2) P_(n-1)'(x) + (2n-3) x P_(n-1)(x)
1280 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1281 final double[] p = new double[order];
1282 p[0] = 1;
1283 final double x2 = x * x;
1284 final double f = 1.0 / (1 - x2);
1285 double coeff = FastMath.sqrt(f);
1286 function[1] = coeff * p[0];
1287 for (int n = 2; n <= order; ++n) {
1288
1289 // update and evaluate polynomial P_n(x)
1290 double v = 0;
1291 p[n - 1] = (n - 1) * p[n - 2];
1292 for (int k = n - 1; k >= 0; k -= 2) {
1293 v = v * x2 + p[k];
1294 if (k > 2) {
1295 p[k - 2] = (k - 1) * p[k - 1] + (2 * n - k) * p[k - 3];
1296 } else if (k == 2) {
1297 p[0] = p[1];
1298 }
1299 }
1300 if ((n & 0x1) == 0) {
1301 v *= x;
1302 }
1303
1304 coeff *= f;
1305 function[n] = coeff * v;
1306
1307 }
1308 }
1309
1310 // apply function composition
1311 compose(operand, operandOffset, function, result, resultOffset);
1312
1313 }
1314
1315 /** Compute arc tangent of a derivative structure.
1316 * @param operand array holding the operand
1317 * @param operandOffset offset of the operand in its array
1318 * @param result array where result must be stored (for
1319 * arc tangent the result array <em>cannot</em> be the input
1320 * array)
1321 * @param resultOffset offset of the result in its array
1322 */
1323 public void atan(final double[] operand, final int operandOffset,
1324 final double[] result, final int resultOffset) {
1325
1326 // create the function value and derivatives
1327 double[] function = new double[1 + order];
1328 final double x = operand[operandOffset];
1329 function[0] = FastMath.atan(x);
1330 if (order > 0) {
1331 // the nth order derivative of atan has the form:
1332 // dn(atan(x)/dxn = Q_n(x) / (1 + x^2)^n
1333 // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
1334 // Q_1(x) = 1, Q_2(x) = -2x, Q_3(x) = 6x^2 - 2 ...
1335 // the general recurrence relation for Q_n is:
1336 // Q_n(x) = (1+x^2) Q_(n-1)'(x) - 2(n-1) x Q_(n-1)(x)
1337 // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
1338 final double[] q = new double[order];
1339 q[0] = 1;
1340 final double x2 = x * x;
1341 final double f = 1.0 / (1 + x2);
1342 double coeff = f;
1343 function[1] = coeff * q[0];
1344 for (int n = 2; n <= order; ++n) {
1345
1346 // update and evaluate polynomial Q_n(x)
1347 double v = 0;
1348 q[n - 1] = -n * q[n - 2];
1349 for (int k = n - 1; k >= 0; k -= 2) {
1350 v = v * x2 + q[k];
1351 if (k > 2) {
1352 q[k - 2] = (k - 1) * q[k - 1] + (k - 1 - 2 * n) * q[k - 3];
1353 } else if (k == 2) {
1354 q[0] = q[1];
1355 }
1356 }
1357 if ((n & 0x1) == 0) {
1358 v *= x;
1359 }
1360
1361 coeff *= f;
1362 function[n] = coeff * v;
1363
1364 }
1365 }
1366
1367 // apply function composition
1368 compose(operand, operandOffset, function, result, resultOffset);
1369
1370 }
1371
1372 /** Compute two arguments arc tangent of a derivative structure.
1373 * @param y array holding the first operand
1374 * @param yOffset offset of the first operand in its array
1375 * @param x array holding the second operand
1376 * @param xOffset offset of the second operand in its array
1377 * @param result array where result must be stored (for
1378 * two arguments arc tangent the result array <em>cannot</em>
1379 * be the input array)
1380 * @param resultOffset offset of the result in its array
1381 */
1382 public void atan2(final double[] y, final int yOffset,
1383 final double[] x, final int xOffset,
1384 final double[] result, final int resultOffset) {
1385
1386 // compute r = sqrt(x^2+y^2)
1387 double[] tmp1 = new double[getSize()];
1388 multiply(x, xOffset, x, xOffset, tmp1, 0); // x^2
1389 double[] tmp2 = new double[getSize()];
1390 multiply(y, yOffset, y, yOffset, tmp2, 0); // y^2
1391 add(tmp1, 0, tmp2, 0, tmp2, 0); // x^2 + y^2
1392 rootN(tmp2, 0, 2, tmp1, 0); // r = sqrt(x^2 + y^2)
1393
1394 if (x[xOffset] >= 0) {
1395
1396 // compute atan2(y, x) = 2 atan(y / (r + x))
1397 add(tmp1, 0, x, xOffset, tmp2, 0); // r + x
1398 divide(y, yOffset, tmp2, 0, tmp1, 0); // y /(r + x)
1399 atan(tmp1, 0, tmp2, 0); // atan(y / (r + x))
1400 for (int i = 0; i < tmp2.length; ++i) {
1401 result[resultOffset + i] = 2 * tmp2[i]; // 2 * atan(y / (r + x))
1402 }
1403
1404 } else {
1405
1406 // compute atan2(y, x) = +/- pi - 2 atan(y / (r - x))
1407 subtract(tmp1, 0, x, xOffset, tmp2, 0); // r - x
1408 divide(y, yOffset, tmp2, 0, tmp1, 0); // y /(r - x)
1409 atan(tmp1, 0, tmp2, 0); // atan(y / (r - x))
1410 result[resultOffset] =
1411 ((tmp2[0] <= 0) ? -FastMath.PI : FastMath.PI) - 2 * tmp2[0]; // +/-pi - 2 * atan(y / (r - x))
1412 for (int i = 1; i < tmp2.length; ++i) {
1413 result[resultOffset + i] = -2 * tmp2[i]; // +/-pi - 2 * atan(y / (r - x))
1414 }
1415
1416 }
1417
1418 }
1419
1420 /** Compute hyperbolic cosine of a derivative structure.
1421 * @param operand array holding the operand
1422 * @param operandOffset offset of the operand in its array
1423 * @param result array where result must be stored (for
1424 * hyperbolic cosine the result array <em>cannot</em> be the input
1425 * array)
1426 * @param resultOffset offset of the result in its array
1427 */
1428 public void cosh(final double[] operand, final int operandOffset,
1429 final double[] result, final int resultOffset) {
1430
1431 // create the function value and derivatives
1432 double[] function = new double[1 + order];
1433 function[0] = FastMath.cosh(operand[operandOffset]);
1434 if (order > 0) {
1435 function[1] = FastMath.sinh(operand[operandOffset]);
1436 for (int i = 2; i <= order; ++i) {
1437 function[i] = function[i - 2];
1438 }
1439 }
1440
1441 // apply function composition
1442 compose(operand, operandOffset, function, result, resultOffset);
1443
1444 }
1445
1446 /** Compute hyperbolic sine of a derivative structure.
1447 * @param operand array holding the operand
1448 * @param operandOffset offset of the operand in its array
1449 * @param result array where result must be stored (for
1450 * hyperbolic sine the result array <em>cannot</em> be the input
1451 * array)
1452 * @param resultOffset offset of the result in its array
1453 */
1454 public void sinh(final double[] operand, final int operandOffset,
1455 final double[] result, final int resultOffset) {
1456
1457 // create the function value and derivatives
1458 double[] function = new double[1 + order];
1459 function[0] = FastMath.sinh(operand[operandOffset]);
1460 if (order > 0) {
1461 function[1] = FastMath.cosh(operand[operandOffset]);
1462 for (int i = 2; i <= order; ++i) {
1463 function[i] = function[i - 2];
1464 }
1465 }
1466
1467 // apply function composition
1468 compose(operand, operandOffset, function, result, resultOffset);
1469
1470 }
1471
1472 /** Compute hyperbolic tangent of a derivative structure.
1473 * @param operand array holding the operand
1474 * @param operandOffset offset of the operand in its array
1475 * @param result array where result must be stored (for
1476 * hyperbolic tangent the result array <em>cannot</em> be the input
1477 * array)
1478 * @param resultOffset offset of the result in its array
1479 */
1480 public void tanh(final double[] operand, final int operandOffset,
1481 final double[] result, final int resultOffset) {
1482
1483 // create the function value and derivatives
1484 final double[] function = new double[1 + order];
1485 final double t = FastMath.tanh(operand[operandOffset]);
1486 function[0] = t;
1487
1488 if (order > 0) {
1489
1490 // the nth order derivative of tanh has the form:
1491 // dn(tanh(x)/dxn = P_n(tanh(x))
1492 // where P_n(t) is a degree n+1 polynomial with same parity as n+1
1493 // P_0(t) = t, P_1(t) = 1 - t^2, P_2(t) = -2 t (1 - t^2) ...
1494 // the general recurrence relation for P_n is:
1495 // P_n(x) = (1-t^2) P_(n-1)'(t)
1496 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1497 final double[] p = new double[order + 2];
1498 p[1] = 1;
1499 final double t2 = t * t;
1500 for (int n = 1; n <= order; ++n) {
1501
1502 // update and evaluate polynomial P_n(t)
1503 double v = 0;
1504 p[n + 1] = -n * p[n];
1505 for (int k = n + 1; k >= 0; k -= 2) {
1506 v = v * t2 + p[k];
1507 if (k > 2) {
1508 p[k - 2] = (k - 1) * p[k - 1] - (k - 3) * p[k - 3];
1509 } else if (k == 2) {
1510 p[0] = p[1];
1511 }
1512 }
1513 if ((n & 0x1) == 0) {
1514 v *= t;
1515 }
1516
1517 function[n] = v;
1518
1519 }
1520 }
1521
1522 // apply function composition
1523 compose(operand, operandOffset, function, result, resultOffset);
1524
1525 }
1526
1527 /** Compute inverse hyperbolic cosine of a derivative structure.
1528 * @param operand array holding the operand
1529 * @param operandOffset offset of the operand in its array
1530 * @param result array where result must be stored (for
1531 * inverse hyperbolic cosine the result array <em>cannot</em> be the input
1532 * array)
1533 * @param resultOffset offset of the result in its array
1534 */
1535 public void acosh(final double[] operand, final int operandOffset,
1536 final double[] result, final int resultOffset) {
1537
1538 // create the function value and derivatives
1539 double[] function = new double[1 + order];
1540 final double x = operand[operandOffset];
1541 function[0] = FastMath.acosh(x);
1542 if (order > 0) {
1543 // the nth order derivative of acosh has the form:
1544 // dn(acosh(x)/dxn = P_n(x) / [x^2 - 1]^((2n-1)/2)
1545 // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1546 // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 + 1 ...
1547 // the general recurrence relation for P_n is:
1548 // P_n(x) = (x^2-1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
1549 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1550 final double[] p = new double[order];
1551 p[0] = 1;
1552 final double x2 = x * x;
1553 final double f = 1.0 / (x2 - 1);
1554 double coeff = FastMath.sqrt(f);
1555 function[1] = coeff * p[0];
1556 for (int n = 2; n <= order; ++n) {
1557
1558 // update and evaluate polynomial P_n(x)
1559 double v = 0;
1560 p[n - 1] = (1 - n) * p[n - 2];
1561 for (int k = n - 1; k >= 0; k -= 2) {
1562 v = v * x2 + p[k];
1563 if (k > 2) {
1564 p[k - 2] = (1 - k) * p[k - 1] + (k - 2 * n) * p[k - 3];
1565 } else if (k == 2) {
1566 p[0] = -p[1];
1567 }
1568 }
1569 if ((n & 0x1) == 0) {
1570 v *= x;
1571 }
1572
1573 coeff *= f;
1574 function[n] = coeff * v;
1575
1576 }
1577 }
1578
1579 // apply function composition
1580 compose(operand, operandOffset, function, result, resultOffset);
1581
1582 }
1583
1584 /** Compute inverse hyperbolic sine of a derivative structure.
1585 * @param operand array holding the operand
1586 * @param operandOffset offset of the operand in its array
1587 * @param result array where result must be stored (for
1588 * inverse hyperbolic sine the result array <em>cannot</em> be the input
1589 * array)
1590 * @param resultOffset offset of the result in its array
1591 */
1592 public void asinh(final double[] operand, final int operandOffset,
1593 final double[] result, final int resultOffset) {
1594
1595 // create the function value and derivatives
1596 double[] function = new double[1 + order];
1597 final double x = operand[operandOffset];
1598 function[0] = FastMath.asinh(x);
1599 if (order > 0) {
1600 // the nth order derivative of asinh has the form:
1601 // dn(asinh(x)/dxn = P_n(x) / [x^2 + 1]^((2n-1)/2)
1602 // where P_n(x) is a degree n-1 polynomial with same parity as n-1
1603 // P_1(x) = 1, P_2(x) = -x, P_3(x) = 2x^2 - 1 ...
1604 // the general recurrence relation for P_n is:
1605 // P_n(x) = (x^2+1) P_(n-1)'(x) - (2n-3) x P_(n-1)(x)
1606 // as per polynomial parity, we can store coefficients of both P_(n-1) and P_n in the same array
1607 final double[] p = new double[order];
1608 p[0] = 1;
1609 final double x2 = x * x;
1610 final double f = 1.0 / (1 + x2);
1611 double coeff = FastMath.sqrt(f);
1612 function[1] = coeff * p[0];
1613 for (int n = 2; n <= order; ++n) {
1614
1615 // update and evaluate polynomial P_n(x)
1616 double v = 0;
1617 p[n - 1] = (1 - n) * p[n - 2];
1618 for (int k = n - 1; k >= 0; k -= 2) {
1619 v = v * x2 + p[k];
1620 if (k > 2) {
1621 p[k - 2] = (k - 1) * p[k - 1] + (k - 2 * n) * p[k - 3];
1622 } else if (k == 2) {
1623 p[0] = p[1];
1624 }
1625 }
1626 if ((n & 0x1) == 0) {
1627 v *= x;
1628 }
1629
1630 coeff *= f;
1631 function[n] = coeff * v;
1632
1633 }
1634 }
1635
1636 // apply function composition
1637 compose(operand, operandOffset, function, result, resultOffset);
1638
1639 }
1640
1641 /** Compute inverse hyperbolic tangent of a derivative structure.
1642 * @param operand array holding the operand
1643 * @param operandOffset offset of the operand in its array
1644 * @param result array where result must be stored (for
1645 * inverse hyperbolic tangent the result array <em>cannot</em> be the input
1646 * array)
1647 * @param resultOffset offset of the result in its array
1648 */
1649 public void atanh(final double[] operand, final int operandOffset,
1650 final double[] result, final int resultOffset) {
1651
1652 // create the function value and derivatives
1653 double[] function = new double[1 + order];
1654 final double x = operand[operandOffset];
1655 function[0] = FastMath.atanh(x);
1656 if (order > 0) {
1657 // the nth order derivative of atanh has the form:
1658 // dn(atanh(x)/dxn = Q_n(x) / (1 - x^2)^n
1659 // where Q_n(x) is a degree n-1 polynomial with same parity as n-1
1660 // Q_1(x) = 1, Q_2(x) = 2x, Q_3(x) = 6x^2 + 2 ...
1661 // the general recurrence relation for Q_n is:
1662 // Q_n(x) = (1-x^2) Q_(n-1)'(x) + 2(n-1) x Q_(n-1)(x)
1663 // as per polynomial parity, we can store coefficients of both Q_(n-1) and Q_n in the same array
1664 final double[] q = new double[order];
1665 q[0] = 1;
1666 final double x2 = x * x;
1667 final double f = 1.0 / (1 - x2);
1668 double coeff = f;
1669 function[1] = coeff * q[0];
1670 for (int n = 2; n <= order; ++n) {
1671
1672 // update and evaluate polynomial Q_n(x)
1673 double v = 0;
1674 q[n - 1] = n * q[n - 2];
1675 for (int k = n - 1; k >= 0; k -= 2) {
1676 v = v * x2 + q[k];
1677 if (k > 2) {
1678 q[k - 2] = (k - 1) * q[k - 1] + (2 * n - k + 1) * q[k - 3];
1679 } else if (k == 2) {
1680 q[0] = q[1];
1681 }
1682 }
1683 if ((n & 0x1) == 0) {
1684 v *= x;
1685 }
1686
1687 coeff *= f;
1688 function[n] = coeff * v;
1689
1690 }
1691 }
1692
1693 // apply function composition
1694 compose(operand, operandOffset, function, result, resultOffset);
1695
1696 }
1697
1698 /** Compute composition of a derivative structure by a function.
1699 * @param operand array holding the operand
1700 * @param operandOffset offset of the operand in its array
1701 * @param f array of value and derivatives of the function at
1702 * the current point (i.e. at {@code operand[operandOffset]}).
1703 * @param result array where result must be stored (for
1704 * composition the result array <em>cannot</em> be the input
1705 * array)
1706 * @param resultOffset offset of the result in its array
1707 */
1708 public void compose(final double[] operand, final int operandOffset, final double[] f,
1709 final double[] result, final int resultOffset) {
1710 for (int i = 0; i < compIndirection.length; ++i) {
1711 final int[][] mappingI = compIndirection[i];
1712 double r = 0;
1713 for (int j = 0; j < mappingI.length; ++j) {
1714 final int[] mappingIJ = mappingI[j];
1715 double product = mappingIJ[0] * f[mappingIJ[1]];
1716 for (int k = 2; k < mappingIJ.length; ++k) {
1717 product *= operand[operandOffset + mappingIJ[k]];
1718 }
1719 r += product;
1720 }
1721 result[resultOffset + i] = r;
1722 }
1723 }
1724
1725 /** Evaluate Taylor expansion of a derivative structure.
1726 * @param ds array holding the derivative structure
1727 * @param dsOffset offset of the derivative structure in its array
1728 * @param delta parameters offsets (Δx, Δy, ...)
1729 * @return value of the Taylor expansion at x + Δx, y + Δy, ...
1730 */
1731 public double taylor(final double[] ds, final int dsOffset, final double ... delta) {
1732 double value = 0;
1733 for (int i = getSize() - 1; i >= 0; --i) {
1734 final int[] orders = getPartialDerivativeOrders(i);
1735 double term = ds[dsOffset + i];
1736 for (int k = 0; k < orders.length; ++k) {
1737 if (orders[k] > 0) {
1738 term *= FastMath.pow(delta[k], orders[k]) / ArithmeticUtils.factorial(orders[k]);
1739 }
1740 }
1741 value += term;
1742 }
1743 return value;
1744 }
1745
1746 /** Check rules set compatibility.
1747 * @param compiler other compiler to check against instance
1748 * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
1749 */
1750 public void checkCompatibility(final DSCompiler compiler)
1751 throws DimensionMismatchException {
1752 if (parameters != compiler.parameters) {
1753 throw new DimensionMismatchException(parameters, compiler.parameters);
1754 }
1755 if (order != compiler.order) {
1756 throw new DimensionMismatchException(order, compiler.order);
1757 }
1758 }
1759
1760 }