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