View Javadoc
1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      https://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  
18  /*
19   * This is not the original file distributed by the Apache Software Foundation
20   * It has been modified by the Hipparchus project
21   */
22  package org.hipparchus.analysis.differentiation;
23  
24  import java.io.Serializable;
25  
26  import org.hipparchus.Field;
27  import org.hipparchus.exception.MathIllegalArgumentException;
28  import org.hipparchus.exception.MathRuntimeException;
29  import org.hipparchus.util.FastMath;
30  import org.hipparchus.util.FieldSinCos;
31  import org.hipparchus.util.FieldSinhCosh;
32  import org.hipparchus.util.MathArrays;
33  import org.hipparchus.util.MathUtils;
34  
35  /** Class representing both the value and the differentials of a function.
36   * <p>This class is the workhorse of the differentiation package.</p>
37   * <p>This class is an implementation of the extension to Rall's
38   * numbers described in Dan Kalman's paper <a
39   * href="http://www.dankalman.net/AUhome/pdffiles/mmgautodiff.pdf">Doubly
40   * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
41   * no. 3, June 2002. Rall's numbers are an extension to the real numbers used
42   * throughout mathematical expressions; they hold the derivative together with the
43   * value of a function. Dan Kalman's derivative structures hold all partial derivatives
44   * up to any specified order, with respect to any number of free parameters. Rall's
45   * numbers therefore can be seen as derivative structures for order one derivative and
46   * one free parameter, and real numbers can be seen as derivative structures with zero
47   * order derivative and no free parameters.</p>
48   * <p>{@link DerivativeStructure} instances can be used directly thanks to
49   * the arithmetic operators to the mathematical functions provided as
50   * methods by this class (+, -, *, /, %, sin, cos ...).</p>
51   * <p>Implementing complex expressions by hand using {@link Derivative}-based
52   * classes (or in fact any {@link org.hipparchus.CalculusFieldElement} class) is
53   * a tedious and error-prone task but has the advantage of not requiring users
54   * to compute the derivatives by themselves and allowing to switch for one
55   * derivative implementation to another as they all share the same filed API.</p>
56   * <p>Implementing complex expression can also be done by developing computation
57   * code using standard primitive double values and to use {@link
58   * UnivariateFunctionDifferentiator differentiators} to create the {@link
59   * DerivativeStructure}-based instances. This method is simpler but may be limited in
60   * the accuracy and derivation orders and may be computationally intensive (this is
61   * typically the case for {@link FiniteDifferencesDifferentiator finite differences
62   * differentiator}.</p>
63   * <p>Instances of this class are guaranteed to be immutable.</p>
64   * @see DSCompiler
65   * @see FieldDerivativeStructure
66   */
67  public class DerivativeStructure implements Derivative<DerivativeStructure>, Serializable {
68  
69      /** Serializable UID. */
70      private static final long serialVersionUID = 20161220L;
71  
72      /** Factory that built the instance. */
73      private final DSFactory factory;
74  
75      /** Combined array holding all values. */
76      private final double[] data;
77  
78      /** Build an instance with all values and derivatives set to 0.
79       * @param factory factory that built the instance
80       * @param data combined array holding all values
81       */
82      DerivativeStructure(final DSFactory factory, final double[] data) {
83          this.factory = factory;
84          this.data    = data.clone();
85      }
86  
87      /** Build an instance with all values and derivatives set to 0.
88       * @param factory factory that built the instance
89       * @since 1.4
90       */
91      DerivativeStructure(final DSFactory factory) {
92          this.factory = factory;
93          this.data    = new double[factory.getCompiler().getSize()];
94      }
95  
96      /** {@inheritDoc} */
97      @Override
98      public DerivativeStructure newInstance(final double value) {
99          return factory.constant(value);
100     }
101 
102     /** {@inheritDoc} */
103     @Override
104     public DerivativeStructure withValue(final double value) {
105         final DerivativeStructure ds = factory.build();
106         System.arraycopy(data, 1, ds.data, 1, data.length - 1);
107         ds.data[0] = value;
108         return ds;
109     }
110 
111     /** Get the factory that built the instance.
112      * @return factory that built the instance
113      */
114     public DSFactory getFactory() {
115         return factory;
116     }
117 
118     /** {@inheritDoc} */
119     @Override
120     public int getFreeParameters() {
121         return getFactory().getCompiler().getFreeParameters();
122     }
123 
124     /** {@inheritDoc} */
125     @Override
126     public int getOrder() {
127         return getFactory().getCompiler().getOrder();
128     }
129 
130     /** Set a derivative component.
131      * <p>
132      * This method is package-private (no modifier specified), as it is intended
133      * to be used only by Hipparchus classes since it relied on the ordering of
134      * derivatives within the class. This allows avoiding checks on the index,
135      * for performance reasons.
136      * </p>
137      * @param index index of the derivative
138      * @param value of the derivative to set
139      * @since 1.4
140      */
141     void setDerivativeComponent(final int index, final double value) {
142         data[index] = value;
143     }
144 
145     /** Get a derivative component.
146      * <p>
147      * This method is package-private (no modifier specified), as it is intended
148      * to be used only by Hipparchus classes since it relied on the ordering of
149      * derivatives within the class. This allows avoiding checks on the index,
150      * for performance reasons.
151      * </p>
152      * @param index index of the derivative
153      * @return value of the derivative
154      * @since 2.2
155      */
156     double getDerivativeComponent(final int index) {
157         return data[index];
158     }
159 
160     /** {@inheritDoc} */
161     @Override
162     public DerivativeStructure getAddendum() {
163         final double[] addendum = data.clone();
164         addendum[0] = 0;
165         return new DerivativeStructure(factory, addendum);
166     }
167 
168     /** Get the value part of the derivative structure.
169      * @return value part of the derivative structure
170      * @see #getPartialDerivative(int...)
171      */
172     @Override
173     public double getValue() {
174         return data[0];
175     }
176 
177     /** {@inheritDoc} */
178     @Override
179     public double getPartialDerivative(final int ... orders)
180         throws MathIllegalArgumentException {
181         return data[getFactory().getCompiler().getPartialDerivativeIndex(orders)];
182     }
183 
184     /** Get all partial derivatives.
185      * @return a fresh copy of partial derivatives, in an array sorted according to
186      * {@link DSCompiler#getPartialDerivativeIndex(int...)}
187      */
188     public double[] getAllDerivatives() {
189         return data.clone();
190     }
191 
192     /** {@inheritDoc}
193      * @exception MathIllegalArgumentException if number of free parameters
194      * or orders do not match
195      */
196     @Override
197     public DerivativeStructure add(final DerivativeStructure a)
198         throws MathIllegalArgumentException {
199         factory.checkCompatibility(a.factory);
200         final DerivativeStructure ds = factory.build();
201         factory.getCompiler().add(data, 0, a.data, 0, ds.data, 0);
202         return ds;
203     }
204 
205     /** {@inheritDoc}
206      * @exception MathIllegalArgumentException if number of free parameters
207      * or orders do not match
208      */
209     @Override
210     public DerivativeStructure subtract(final DerivativeStructure a)
211         throws MathIllegalArgumentException {
212         factory.checkCompatibility(a.factory);
213         final DerivativeStructure ds = factory.build();
214         factory.getCompiler().subtract(data, 0, a.data, 0, ds.data, 0);
215         return ds;
216     }
217 
218     /** {@inheritDoc}
219      */
220     @Override
221     public DerivativeStructure multiply(final double a) {
222         final DerivativeStructure ds = factory.build();
223         for (int i = 0; i < ds.data.length; ++i) {
224             ds.data[i] = data[i] * a;
225         }
226         return ds;
227     }
228 
229     /** {@inheritDoc}
230      * @exception MathIllegalArgumentException if number of free parameters
231      * or orders do not match
232      */
233     @Override
234     public DerivativeStructure multiply(final DerivativeStructure a)
235         throws MathIllegalArgumentException {
236         factory.checkCompatibility(a.factory);
237         final DerivativeStructure result = factory.build();
238         factory.getCompiler().multiply(data, 0, a.data, 0, result.data, 0);
239         return result;
240     }
241 
242     /** {@inheritDoc} */
243     @Override
244     public DerivativeStructure square() {
245         return multiply(this);
246     }
247 
248     /** {@inheritDoc}
249      */
250     @Override
251     public DerivativeStructure divide(final double a) {
252         final DerivativeStructure ds = factory.build();
253         final double inv = 1.0 / a;
254         for (int i = 0; i < ds.data.length; ++i) {
255             ds.data[i] = data[i] * inv;
256         }
257         return ds;
258     }
259 
260     /** {@inheritDoc}
261      * @exception MathIllegalArgumentException if number of free parameters
262      * or orders do not match
263      */
264     @Override
265     public DerivativeStructure divide(final DerivativeStructure a)
266         throws MathIllegalArgumentException {
267         factory.checkCompatibility(a.factory);
268         final DerivativeStructure result = factory.build();
269         factory.getCompiler().divide(data, 0, a.data, 0, result.data, 0);
270         return result;
271     }
272 
273     /** {@inheritDoc}
274      * @exception MathIllegalArgumentException if number of free parameters
275      * or orders do not match
276      */
277     @Override
278     public DerivativeStructure remainder(final DerivativeStructure a)
279         throws MathIllegalArgumentException {
280         factory.checkCompatibility(a.factory);
281         final DerivativeStructure result = factory.build();
282         factory.getCompiler().remainder(data, 0, a.data, 0, result.data, 0);
283         return result;
284     }
285 
286     /** {@inheritDoc} */
287     @Override
288     public DerivativeStructure negate() {
289         final DerivativeStructure ds = factory.build();
290         for (int i = 0; i < ds.data.length; ++i) {
291             ds.data[i] = -data[i];
292         }
293         return ds;
294     }
295 
296     /** {@inheritDoc}
297      */
298     @Override
299     public DerivativeStructure abs() {
300         if (Double.doubleToLongBits(data[0]) < 0) {
301             // we use the bits representation to also handle -0.0
302             return negate();
303         } else {
304             return this;
305         }
306     }
307 
308     /** {@inheritDoc}
309      */
310     @Override
311     public DerivativeStructure copySign(final DerivativeStructure sign) {
312         long m = Double.doubleToLongBits(data[0]);
313         long s = Double.doubleToLongBits(sign.data[0]);
314         if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
315             return this;
316         }
317         return negate(); // flip sign
318     }
319 
320     /** {@inheritDoc}
321      */
322     @Override
323     public DerivativeStructure copySign(final double sign) {
324         long m = Double.doubleToLongBits(data[0]);
325         long s = Double.doubleToLongBits(sign);
326         if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
327             return this;
328         }
329         return negate(); // flip sign
330     }
331 
332     /** {@inheritDoc}
333      */
334     @Override
335     public DerivativeStructure scalb(final int n) {
336         final DerivativeStructure ds = factory.build();
337         for (int i = 0; i < ds.data.length; ++i) {
338             ds.data[i] = FastMath.scalb(data[i], n);
339         }
340         return ds;
341     }
342 
343     /** {@inheritDoc}
344      * @exception MathIllegalArgumentException if number of free parameters
345      * or orders do not match
346      */
347     @Override
348     public DerivativeStructure hypot(final DerivativeStructure y)
349         throws MathIllegalArgumentException {
350 
351         factory.checkCompatibility(y.factory);
352 
353         if (Double.isInfinite(data[0]) || Double.isInfinite(y.data[0])) {
354             return factory.constant(Double.POSITIVE_INFINITY);
355         } else if (Double.isNaN(data[0]) || Double.isNaN(y.data[0])) {
356             return factory.constant(Double.NaN);
357         } else {
358 
359             final int expX = getExponent();
360             final int expY = y.getExponent();
361             if (expX > expY + 27) {
362                 // y is negligible with respect to x
363                 return abs();
364             } else if (expY > expX + 27) {
365                 // x is negligible with respect to y
366                 return y.abs();
367             } else {
368 
369                 // find an intermediate scale to avoid both overflow and underflow
370                 final int middleExp = (expX + expY) / 2;
371 
372                 // scale parameters without losing precision
373                 final DerivativeStructure scaledX = scalb(-middleExp);
374                 final DerivativeStructure scaledY = y.scalb(-middleExp);
375 
376                 // compute scaled hypotenuse
377                 final DerivativeStructure scaledH =
378                         scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
379 
380                 // remove scaling
381                 return scaledH.scalb(middleExp);
382 
383             }
384 
385         }
386     }
387 
388     /**
389      * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
390      * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
391      * avoiding intermediate overflow or underflow.
392      *
393      * <ul>
394      * <li> If either argument is infinite, then the result is positive infinity.</li>
395      * <li> else, if either argument is NaN then the result is NaN.</li>
396      * </ul>
397      *
398      * @param x a value
399      * @param y a value
400      * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
401      * @exception MathIllegalArgumentException if number of free parameters
402      * or orders do not match
403      */
404     public static DerivativeStructure hypot(final DerivativeStructure x, final DerivativeStructure y)
405         throws MathIllegalArgumentException {
406         return x.hypot(y);
407     }
408 
409     /** Compute composition of the instance by a univariate function.
410      * @param f array of value and derivatives of the function at
411      * the current point (i.e. [f({@link #getValue()}),
412      * f'({@link #getValue()}), f''({@link #getValue()})...]).
413      * @return f(this)
414      * @exception MathIllegalArgumentException if the number of derivatives
415      * in the array is not equal to {@link #getOrder() order} + 1
416      */
417     @Override
418     public DerivativeStructure compose(final double ... f)
419         throws MathIllegalArgumentException {
420 
421         MathUtils.checkDimension(f.length, getOrder() + 1);
422         final DerivativeStructure result = factory.build();
423         factory.getCompiler().compose(data, 0, f, result.data, 0);
424         return result;
425     }
426 
427     /** {@inheritDoc} */
428     @Override
429     public DerivativeStructure reciprocal() {
430         final DerivativeStructure result = factory.build();
431         factory.getCompiler().reciprocal(data, 0, result.data, 0);
432         return result;
433     }
434 
435     /** {@inheritDoc}
436      */
437     @Override
438     public DerivativeStructure sqrt() {
439         final DerivativeStructure result = factory.build();
440         factory.getCompiler().sqrt(data, 0, result.data, 0);
441         return result;
442     }
443 
444     /** {@inheritDoc}
445      */
446     @Override
447     public DerivativeStructure rootN(final int n) {
448         final DerivativeStructure result = factory.build();
449         factory.getCompiler().rootN(data, 0, n, result.data, 0);
450         return result;
451     }
452 
453     /** {@inheritDoc} */
454     @Override
455     public Field<DerivativeStructure> getField() {
456         return factory.getDerivativeField();
457     }
458 
459     /** Compute a<sup>x</sup> where a is a double and x a {@link DerivativeStructure}
460      * @param a number to exponentiate
461      * @param x power to apply
462      * @return a<sup>x</sup>
463      */
464     public static DerivativeStructure pow(final double a, final DerivativeStructure x) {
465         final DerivativeStructure result = x.factory.build();
466         x.factory.getCompiler().pow(a, x.data, 0, result.data, 0);
467         return result;
468     }
469 
470     /** {@inheritDoc}
471      */
472     @Override
473     public DerivativeStructure pow(final double p) {
474         final DerivativeStructure result = factory.build();
475         factory.getCompiler().pow(data, 0, p, result.data, 0);
476         return result;
477     }
478 
479     /** {@inheritDoc}
480      */
481     @Override
482     public DerivativeStructure pow(final int n) {
483         final DerivativeStructure result = factory.build();
484         factory.getCompiler().pow(data, 0, n, result.data, 0);
485         return result;
486     }
487 
488     /** {@inheritDoc}
489      * @exception MathIllegalArgumentException if number of free parameters
490      * or orders do not match
491      */
492     @Override
493     public DerivativeStructure pow(final DerivativeStructure e)
494         throws MathIllegalArgumentException {
495         factory.checkCompatibility(e.factory);
496         final DerivativeStructure result = factory.build();
497         factory.getCompiler().pow(data, 0, e.data, 0, result.data, 0);
498         return result;
499     }
500 
501     /** {@inheritDoc}
502      */
503     @Override
504     public DerivativeStructure exp() {
505         final DerivativeStructure result = factory.build();
506         factory.getCompiler().exp(data, 0, result.data, 0);
507         return result;
508     }
509 
510     /** {@inheritDoc}
511      */
512     @Override
513     public DerivativeStructure expm1() {
514         final DerivativeStructure result = factory.build();
515         factory.getCompiler().expm1(data, 0, result.data, 0);
516         return result;
517     }
518 
519     /** {@inheritDoc}
520      */
521     @Override
522     public DerivativeStructure log() {
523         final DerivativeStructure result = factory.build();
524         factory.getCompiler().log(data, 0, result.data, 0);
525         return result;
526     }
527 
528     /** {@inheritDoc}
529      */
530     @Override
531     public DerivativeStructure log1p() {
532         final DerivativeStructure result = factory.build();
533         factory.getCompiler().log1p(data, 0, result.data, 0);
534         return result;
535     }
536 
537     /** Base 10 logarithm.
538      * @return base 10 logarithm of the instance
539      */
540     @Override
541     public DerivativeStructure log10() {
542         final DerivativeStructure result = factory.build();
543         factory.getCompiler().log10(data, 0, result.data, 0);
544         return result;
545     }
546 
547     /** {@inheritDoc}
548      */
549     @Override
550     public DerivativeStructure cos() {
551         final DerivativeStructure result = factory.build();
552         factory.getCompiler().cos(data, 0, result.data, 0);
553         return result;
554     }
555 
556     /** {@inheritDoc}
557      */
558     @Override
559     public DerivativeStructure sin() {
560         final DerivativeStructure result = factory.build();
561         factory.getCompiler().sin(data, 0, result.data, 0);
562         return result;
563     }
564 
565     /** {@inheritDoc}
566      */
567     @Override
568     public FieldSinCos<DerivativeStructure> sinCos() {
569         final DerivativeStructure sin = factory.build();
570         final DerivativeStructure cos = factory.build();
571         factory.getCompiler().sinCos(data, 0, sin.data, 0, cos.data, 0);
572         return new FieldSinCos<>(sin, cos);
573     }
574 
575     /** {@inheritDoc}
576      */
577     @Override
578     public DerivativeStructure tan() {
579         final DerivativeStructure result = factory.build();
580         factory.getCompiler().tan(data, 0, result.data, 0);
581         return result;
582     }
583 
584     /** {@inheritDoc}
585      */
586     @Override
587     public DerivativeStructure acos() {
588         final DerivativeStructure result = factory.build();
589         factory.getCompiler().acos(data, 0, result.data, 0);
590         return result;
591     }
592 
593     /** {@inheritDoc}
594      */
595     @Override
596     public DerivativeStructure asin() {
597         final DerivativeStructure result = factory.build();
598         factory.getCompiler().asin(data, 0, result.data, 0);
599         return result;
600     }
601 
602     /** {@inheritDoc}
603      */
604     @Override
605     public DerivativeStructure atan() {
606         final DerivativeStructure result = factory.build();
607         factory.getCompiler().atan(data, 0, result.data, 0);
608         return result;
609     }
610 
611     /** {@inheritDoc}
612      */
613     @Override
614     public DerivativeStructure atan2(final DerivativeStructure x)
615         throws MathIllegalArgumentException {
616         factory.checkCompatibility(x.factory);
617         final DerivativeStructure result = factory.build();
618         factory.getCompiler().atan2(data, 0, x.data, 0, result.data, 0);
619         return result;
620     }
621 
622     /** Two arguments arc tangent operation.
623      * @param y first argument of the arc tangent
624      * @param x second argument of the arc tangent
625      * @return atan2(y, x)
626      * @exception MathIllegalArgumentException if number of free parameters
627      * or orders do not match
628      */
629     public static DerivativeStructure atan2(final DerivativeStructure y, final DerivativeStructure x)
630         throws MathIllegalArgumentException {
631         return y.atan2(x);
632     }
633 
634     /** {@inheritDoc}
635      */
636     @Override
637     public DerivativeStructure cosh() {
638         final DerivativeStructure result = factory.build();
639         factory.getCompiler().cosh(data, 0, result.data, 0);
640         return result;
641     }
642 
643     /** {@inheritDoc}
644      */
645     @Override
646     public DerivativeStructure sinh() {
647         final DerivativeStructure result = factory.build();
648         factory.getCompiler().sinh(data, 0, result.data, 0);
649         return result;
650     }
651 
652     /** {@inheritDoc}
653      */
654     @Override
655     public FieldSinhCosh<DerivativeStructure> sinhCosh() {
656         final DerivativeStructure sinh = factory.build();
657         final DerivativeStructure cosh = factory.build();
658         factory.getCompiler().sinhCosh(data, 0, sinh.data, 0, cosh.data, 0);
659         return new FieldSinhCosh<>(sinh, cosh);
660     }
661 
662     /** {@inheritDoc}
663      */
664     @Override
665     public DerivativeStructure tanh() {
666         final DerivativeStructure result = factory.build();
667         factory.getCompiler().tanh(data, 0, result.data, 0);
668         return result;
669     }
670 
671     /** {@inheritDoc}
672      */
673     @Override
674     public DerivativeStructure acosh() {
675         final DerivativeStructure result = factory.build();
676         factory.getCompiler().acosh(data, 0, result.data, 0);
677         return result;
678     }
679 
680     /** {@inheritDoc}
681      */
682     @Override
683     public DerivativeStructure asinh() {
684         final DerivativeStructure result = factory.build();
685         factory.getCompiler().asinh(data, 0, result.data, 0);
686         return result;
687     }
688 
689     /** {@inheritDoc}
690      */
691     @Override
692     public DerivativeStructure atanh() {
693         final DerivativeStructure result = factory.build();
694         factory.getCompiler().atanh(data, 0, result.data, 0);
695         return result;
696     }
697 
698     /** {@inheritDoc} */
699     @Override
700     public DerivativeStructure toDegrees() {
701         final DerivativeStructure ds = factory.build();
702         for (int i = 0; i < ds.data.length; ++i) {
703             ds.data[i] = FastMath.toDegrees(data[i]);
704         }
705         return ds;
706     }
707 
708     /** {@inheritDoc} */
709     @Override
710     public DerivativeStructure toRadians() {
711         final DerivativeStructure ds = factory.build();
712         for (int i = 0; i < ds.data.length; ++i) {
713             ds.data[i] = FastMath.toRadians(data[i]);
714         }
715         return ds;
716     }
717 
718     /** Integrate w.r.t. one independent variable.
719      * <p>
720      * Rigorously, if the derivatives of a function are known up to
721      * order N, the ones of its M-th integral w.r.t. a given variable
722      * (seen as a function itself) are actually known up to order N+M.
723      * However, this method still casts the output as a DerivativeStructure
724      * of order N. The integration constants are systematically set to zero.
725      * </p>
726      * @param varIndex Index of independent variable w.r.t. which integration is done.
727      * @param integrationOrder Number of times the integration operator must be applied. If non-positive, call the
728      *                         differentiation operator.
729      * @return DerivativeStructure on which integration operator has been applied a certain number of times.
730      * @since 2.2
731      */
732     public DerivativeStructure integrate(final int varIndex, final int integrationOrder) {
733 
734         // Deal first with trivial case
735         if (integrationOrder > getOrder()) {
736             return factory.constant(0.);
737         } else if (integrationOrder == 0) {
738             return factory.build(data);
739         }
740 
741         // Call 'inverse' (not rigorously) operation if necessary
742         if (integrationOrder < 0) {
743             return differentiate(varIndex, -integrationOrder);
744         }
745 
746         final double[] newData = new double[data.length];
747         final DSCompiler dsCompiler = factory.getCompiler();
748         for (int i = 0; i < newData.length; i++) {
749             if (data[i] != 0.) {
750                 final int[] orders = dsCompiler.getPartialDerivativeOrders(i);
751                 int sum = 0;
752                 for (int order : orders) {
753                     sum += order;
754                 }
755                 if (sum + integrationOrder <= getOrder()) {
756                     final int saved = orders[varIndex];
757                     orders[varIndex] += integrationOrder;
758                     final int index = dsCompiler.getPartialDerivativeIndex(orders);
759                     orders[varIndex] = saved;
760                     newData[index] = data[i];
761                 }
762             }
763         }
764 
765         return factory.build(newData);
766     }
767 
768     /** Differentiate w.r.t. one independent variable.
769      * <p>
770      * Rigorously, if the derivatives of a function are known up to
771      * order N, the ones of its M-th derivative w.r.t. a given variable
772      * (seen as a function itself) are only known up to order N-M.
773      * However, this method still casts the output as a DerivativeStructure
774      * of order N with zeroes for the higher order terms.
775      * </p>
776      * @param varIndex Index of independent variable w.r.t. which differentiation is done.
777      * @param differentiationOrder Number of times the differentiation operator must be applied. If non-positive, call
778      *                             the integration operator instead.
779      * @return DerivativeStructure on which differentiation operator has been applied a certain number of times
780      * @since 2.2
781      */
782     public DerivativeStructure differentiate(final int varIndex, final int differentiationOrder) {
783 
784         // Deal first with trivial case
785         if (differentiationOrder > getOrder()) {
786             return factory.constant(0.);
787         } else if (differentiationOrder == 0) {
788             return factory.build(data);
789         }
790 
791         // Call 'inverse' (not rigorously) operation if necessary
792         if (differentiationOrder < 0) {
793             return integrate(varIndex, -differentiationOrder);
794         }
795 
796         final double[] newData = new double[data.length];
797         final DSCompiler dsCompiler = factory.getCompiler();
798         for (int i = 0; i < newData.length; i++) {
799             if (data[i] != 0.) {
800                 final int[] orders = dsCompiler.getPartialDerivativeOrders(i);
801                 if (orders[varIndex] - differentiationOrder >= 0) {
802                     final int saved = orders[varIndex];
803                     orders[varIndex] -= differentiationOrder;
804                     final int index = dsCompiler.getPartialDerivativeIndex(orders);
805                     orders[varIndex] = saved;
806                     newData[index] = data[i];
807                 }
808             }
809         }
810 
811         return factory.build(newData);
812     }
813 
814     /** Evaluate Taylor expansion a derivative structure.
815      * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
816      * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
817      * @throws MathRuntimeException if factorials becomes too large
818      */
819     public double taylor(final double ... delta) throws MathRuntimeException {
820         return factory.getCompiler().taylor(data, 0, delta);
821     }
822 
823     /** Rebase instance with respect to low level parameter functions.
824      * <p>
825      * The instance is considered to be a function of {@link #getFreeParameters()
826      * n free parameters} up to order {@link #getOrder() o} \(f(p_0, p_1, \ldots p_{n-1})\).
827      * Its {@link #getPartialDerivative(int...) partial derivatives} are therefore
828      * \(f, \frac{\partial f}{\partial p_0}, \frac{\partial f}{\partial p_1}, \ldots
829      * \frac{\partial^2 f}{\partial p_0^2}, \frac{\partial^2 f}{\partial p_0 p_1},
830      * \ldots \frac{\partial^o f}{\partial p_{n-1}^o}\). The free parameters
831      * \(p_0, p_1, \ldots p_{n-1}\) are considered to be functions of \(m\) lower
832      * level other parameters \(q_0, q_1, \ldots q_{m-1}\).
833      * </p>
834      * \( \begin{align}
835      * p_0 &amp; = p_0(q_0, q_1, \ldots q_{m-1})\\
836      * p_1 &amp; = p_1(q_0, q_1, \ldots q_{m-1})\\
837      * p_{n-1} &amp; = p_{n-1}(q_0, q_1, \ldots q_{m-1})
838      * \end{align}\)
839      * <p>
840      * This method compute the composition of the partial derivatives of \(f\)
841      * and the partial derivatives of \(p_0, p_1, \ldots p_{n-1}\), i.e. the
842      * {@link #getPartialDerivative(int...) partial derivatives} of the value
843      * returned will be
844      * \(f, \frac{\partial f}{\partial q_0}, \frac{\partial f}{\partial q_1}, \ldots
845      * \frac{\partial^2 f}{\partial q_0^2}, \frac{\partial^2 f}{\partial q_0 q_1},
846      * \ldots \frac{\partial^o f}{\partial q_{m-1}^o}\).
847      * </p>
848      * <p>
849      * The number of parameters must match {@link #getFreeParameters()} and the
850      * derivation orders of the instance and parameters must also match.
851      * </p>
852      * @param p base parameters with respect to which partial derivatives
853      * were computed in the instance
854      * @return derivative structure with partial derivatives computed
855      * with respect to the lower level parameters used in the \(p_i\)
856      * @since 2.2
857      */
858     public DerivativeStructure rebase(final DerivativeStructure... p) {
859 
860         MathUtils.checkDimension(getFreeParameters(), p.length);
861 
862         // handle special case of no variables at all
863         if (p.length == 0) {
864             return this;
865         }
866 
867         final int pSize = p[0].getFactory().getCompiler().getSize();
868         final double[] pData = new double[p.length * pSize];
869         for (int i = 0; i < p.length; ++i) {
870             MathUtils.checkDimension(getOrder(), p[i].getOrder());
871             MathUtils.checkDimension(p[0].getFreeParameters(), p[i].getFreeParameters());
872             System.arraycopy(p[i].data, 0, pData, i * pSize, pSize);
873         }
874 
875         final DerivativeStructure result = p[0].factory.build();
876         factory.getCompiler().rebase(data, 0, p[0].factory.getCompiler(), pData, result.data, 0);
877         return result;
878 
879     }
880 
881     /** {@inheritDoc}
882      * @exception MathIllegalArgumentException if number of free parameters
883      * or orders do not match
884      */
885     @Override
886     public DerivativeStructure linearCombination(final DerivativeStructure[] a, final DerivativeStructure[] b)
887         throws MathIllegalArgumentException {
888 
889         // compute an accurate value, taking care of cancellations
890         final double[] aDouble = new double[a.length];
891         for (int i = 0; i < a.length; ++i) {
892             aDouble[i] = a[i].getValue();
893         }
894         final double[] bDouble = new double[b.length];
895         for (int i = 0; i < b.length; ++i) {
896             bDouble[i] = b[i].getValue();
897         }
898         final double accurateValue = MathArrays.linearCombination(aDouble, bDouble);
899 
900         // compute a simple value, with all partial derivatives
901         DerivativeStructure simpleValue = a[0].getField().getZero();
902         for (int i = 0; i < a.length; ++i) {
903             simpleValue = simpleValue.add(a[i].multiply(b[i]));
904         }
905 
906         // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
907         final double[] all = simpleValue.getAllDerivatives();
908         all[0] = accurateValue;
909         return factory.build(all);
910 
911     }
912 
913     /** {@inheritDoc}
914      * @exception MathIllegalArgumentException if number of free parameters
915      * or orders do not match
916      */
917     @Override
918     public DerivativeStructure linearCombination(final double[] a, final DerivativeStructure[] b)
919         throws MathIllegalArgumentException {
920 
921         // compute an accurate value, taking care of cancellations
922         final double[] bDouble = new double[b.length];
923         for (int i = 0; i < b.length; ++i) {
924             bDouble[i] = b[i].getValue();
925         }
926         final double accurateValue = MathArrays.linearCombination(a, bDouble);
927 
928         // compute a simple value, with all partial derivatives
929         DerivativeStructure simpleValue = b[0].getField().getZero();
930         for (int i = 0; i < a.length; ++i) {
931             simpleValue = simpleValue.add(b[i].multiply(a[i]));
932         }
933 
934         // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
935         final double[] all = simpleValue.getAllDerivatives();
936         all[0] = accurateValue;
937         return factory.build(all);
938 
939     }
940 
941     /** {@inheritDoc}
942      * @exception MathIllegalArgumentException if number of free parameters
943      * or orders do not match
944      */
945     @Override
946     public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
947                                                  final DerivativeStructure a2, final DerivativeStructure b2)
948         throws MathIllegalArgumentException {
949 
950         // compute an accurate value, taking care of cancellations
951         final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
952                                                                   a2.getValue(), b2.getValue());
953 
954         // compute a simple value, with all partial derivatives
955         final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2));
956 
957         // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
958         final double[] all = simpleValue.getAllDerivatives();
959         all[0] = accurateValue;
960         return factory.build(all);
961 
962     }
963 
964     /** {@inheritDoc}
965      * @exception MathIllegalArgumentException if number of free parameters
966      * or orders do not match
967      */
968     @Override
969     public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
970                                                  final double a2, final DerivativeStructure b2)
971         throws MathIllegalArgumentException {
972 
973         factory.checkCompatibility(b1.factory);
974         factory.checkCompatibility(b2.factory);
975 
976         final DerivativeStructure ds = factory.build();
977         factory.getCompiler().linearCombination(a1, b1.data, 0,
978                                                 a2, b2.data, 0,
979                                                 ds.data, 0);
980 
981         return ds;
982 
983     }
984 
985     /** {@inheritDoc}
986      * @exception MathIllegalArgumentException if number of free parameters
987      * or orders do not match
988      */
989     @Override
990     public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
991                                                  final DerivativeStructure a2, final DerivativeStructure b2,
992                                                  final DerivativeStructure a3, final DerivativeStructure b3)
993         throws MathIllegalArgumentException {
994 
995         // compute an accurate value, taking care of cancellations
996         final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
997                                                                   a2.getValue(), b2.getValue(),
998                                                                   a3.getValue(), b3.getValue());
999 
1000         // compute a simple value, with all partial derivatives
1001         final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
1002 
1003         // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1004         final double[] all = simpleValue.getAllDerivatives();
1005         all[0] = accurateValue;
1006         return factory.build(all);
1007 
1008     }
1009 
1010     /** {@inheritDoc}
1011      * @exception MathIllegalArgumentException if number of free parameters
1012      * or orders do not match
1013      */
1014     @Override
1015     public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1016                                                  final double a2, final DerivativeStructure b2,
1017                                                  final double a3, final DerivativeStructure b3)
1018         throws MathIllegalArgumentException {
1019 
1020         factory.checkCompatibility(b1.factory);
1021         factory.checkCompatibility(b2.factory);
1022         factory.checkCompatibility(b3.factory);
1023 
1024         final DerivativeStructure ds = factory.build();
1025         factory.getCompiler().linearCombination(a1, b1.data, 0,
1026                                                 a2, b2.data, 0,
1027                                                 a3, b3.data, 0,
1028                                                 ds.data, 0);
1029 
1030         return ds;
1031 
1032     }
1033 
1034     /** {@inheritDoc}
1035      * @exception MathIllegalArgumentException if number of free parameters
1036      * or orders do not match
1037      */
1038     @Override
1039     public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
1040                                                  final DerivativeStructure a2, final DerivativeStructure b2,
1041                                                  final DerivativeStructure a3, final DerivativeStructure b3,
1042                                                  final DerivativeStructure a4, final DerivativeStructure b4)
1043         throws MathIllegalArgumentException {
1044 
1045         // compute an accurate value, taking care of cancellations
1046         final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
1047                                                                   a2.getValue(), b2.getValue(),
1048                                                                   a3.getValue(), b3.getValue(),
1049                                                                   a4.getValue(), b4.getValue());
1050 
1051         // compute a simple value, with all partial derivatives
1052         final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
1053 
1054         // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1055         final double[] all = simpleValue.getAllDerivatives();
1056         all[0] = accurateValue;
1057         return factory.build(all);
1058 
1059     }
1060 
1061     /** {@inheritDoc}
1062      * @exception MathIllegalArgumentException if number of free parameters
1063      * or orders do not match
1064      */
1065     @Override
1066     public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1067                                                  final double a2, final DerivativeStructure b2,
1068                                                  final double a3, final DerivativeStructure b3,
1069                                                  final double a4, final DerivativeStructure b4)
1070         throws MathIllegalArgumentException {
1071 
1072         factory.checkCompatibility(b1.factory);
1073         factory.checkCompatibility(b2.factory);
1074         factory.checkCompatibility(b3.factory);
1075         factory.checkCompatibility(b4.factory);
1076 
1077         final DerivativeStructure ds = factory.build();
1078         factory.getCompiler().linearCombination(a1, b1.data, 0,
1079                                                 a2, b2.data, 0,
1080                                                 a3, b3.data, 0,
1081                                                 a4, b4.data, 0,
1082                                                 ds.data, 0);
1083 
1084         return ds;
1085 
1086     }
1087 
1088     /** {@inheritDoc}
1089      */
1090     @Override
1091     public DerivativeStructure getPi() {
1092         return factory.getDerivativeField().getPi();
1093     }
1094 
1095     /**
1096      * Test for the equality of two derivative structures.
1097      * <p>
1098      * Derivative structures are considered equal if they have the same number
1099      * of free parameters, the same derivation order, and the same derivatives.
1100      * </p>
1101      * @param other Object to test for equality to this
1102      * @return true if two derivative structures are equal
1103      */
1104     @Override
1105     public boolean equals(Object other) {
1106 
1107         if (this == other) {
1108             return true;
1109         }
1110 
1111         if (other instanceof DerivativeStructure) {
1112             final DerivativeStructure rhs = (DerivativeStructure)other;
1113             return (getFreeParameters() == rhs.getFreeParameters()) &&
1114                    (getOrder() == rhs.getOrder()) &&
1115                    MathArrays.equals(data, rhs.data);
1116         }
1117 
1118         return false;
1119 
1120     }
1121 
1122     /**
1123      * Get a hashCode for the derivative structure.
1124      * @return a hash code value for this object
1125      */
1126     @Override
1127     public int hashCode() {
1128         return 227 + 229 * getFreeParameters() + 233 * getOrder() + 239 * MathUtils.hash(data);
1129     }
1130 
1131 }