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 */
017package org.apache.commons.math3.analysis.differentiation;
018
019import java.io.Serializable;
020
021import org.apache.commons.math3.RealFieldElement;
022import org.apache.commons.math3.Field;
023import org.apache.commons.math3.FieldElement;
024import org.apache.commons.math3.exception.DimensionMismatchException;
025import org.apache.commons.math3.exception.MathArithmeticException;
026import org.apache.commons.math3.exception.NumberIsTooLargeException;
027import org.apache.commons.math3.util.FastMath;
028import org.apache.commons.math3.util.MathArrays;
029import org.apache.commons.math3.util.MathUtils;
030
031/** Class representing both the value and the differentials of a function.
032 * <p>This class is the workhorse of the differentiation package.</p>
033 * <p>This class is an implementation of the extension to Rall's
034 * numbers described in Dan Kalman's paper <a
035 * href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
036 * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
037 * no. 3, June 2002.</p>. Rall's numbers are an extension to the real numbers used
038 * throughout mathematical expressions; they hold the derivative together with the
039 * value of a function. Dan Kalman's derivative structures hold all partial derivatives
040 * up to any specified order, with respect to any number of free parameters. Rall's
041 * numbers therefore can be seen as derivative structures for order one derivative and
042 * one free parameter, and real numbers can be seen as derivative structures with zero
043 * order derivative and no free parameters.</p>
044 * <p>{@link DerivativeStructure} instances can be used directly thanks to
045 * the arithmetic operators to the mathematical functions provided as
046 * methods by this class (+, -, *, /, %, sin, cos ...).</p>
047 * <p>Implementing complex expressions by hand using these classes is
048 * a tedious and error-prone task but has the advantage of having no limitation
049 * on the derivation order despite no requiring users to compute the derivatives by
050 * themselves. Implementing complex expression can also be done by developing computation
051 * code using standard primitive double values and to use {@link
052 * UnivariateFunctionDifferentiator differentiators} to create the {@link
053 * DerivativeStructure}-based instances. This method is simpler but may be limited in
054 * the accuracy and derivation orders and may be computationally intensive (this is
055 * typically the case for {@link FiniteDifferencesDifferentiator finite differences
056 * differentiator}.</p>
057 * <p>Instances of this class are guaranteed to be immutable.</p>
058 * @see DSCompiler
059 * @since 3.1
060 */
061public class DerivativeStructure implements RealFieldElement<DerivativeStructure>, Serializable {
062
063    /** Serializable UID. */
064    private static final long serialVersionUID = 20120730L;
065
066    /** Compiler for the current dimensions. */
067    private transient DSCompiler compiler;
068
069    /** Combined array holding all values. */
070    private final double[] data;
071
072    /** Build an instance with all values and derivatives set to 0.
073     * @param compiler compiler to use for computation
074     */
075    private DerivativeStructure(final DSCompiler compiler) {
076        this.compiler = compiler;
077        this.data     = new double[compiler.getSize()];
078    }
079
080    /** Build an instance with all values and derivatives set to 0.
081     * @param parameters number of free parameters
082     * @param order derivation order
083     * @throws NumberIsTooLargeException if order is too large
084     */
085    public DerivativeStructure(final int parameters, final int order)
086        throws NumberIsTooLargeException {
087        this(DSCompiler.getCompiler(parameters, order));
088    }
089
090    /** Build an instance representing a constant value.
091     * @param parameters number of free parameters
092     * @param order derivation order
093     * @param value value of the constant
094     * @throws NumberIsTooLargeException if order is too large
095     * @see #DerivativeStructure(int, int, int, double)
096     */
097    public DerivativeStructure(final int parameters, final int order, final double value)
098        throws NumberIsTooLargeException {
099        this(parameters, order);
100        this.data[0] = value;
101    }
102
103    /** Build an instance representing a variable.
104     * <p>Instances built using this constructor are considered
105     * to be the free variables with respect to which differentials
106     * are computed. As such, their differential with respect to
107     * themselves is +1.</p>
108     * @param parameters number of free parameters
109     * @param order derivation order
110     * @param index index of the variable (from 0 to {@code parameters - 1})
111     * @param value value of the variable
112     * @exception NumberIsTooLargeException if {@code index >= parameters}.
113     * @see #DerivativeStructure(int, int, double)
114     */
115    public DerivativeStructure(final int parameters, final int order,
116                               final int index, final double value)
117        throws NumberIsTooLargeException {
118        this(parameters, order, value);
119
120        if (index >= parameters) {
121            throw new NumberIsTooLargeException(index, parameters, false);
122        }
123
124        if (order > 0) {
125            // the derivative of the variable with respect to itself is 1.
126            data[DSCompiler.getCompiler(index, order).getSize()] = 1.0;
127        }
128
129    }
130
131    /** Linear combination constructor.
132     * The derivative structure built will be a1 * ds1 + a2 * ds2
133     * @param a1 first scale factor
134     * @param ds1 first base (unscaled) derivative structure
135     * @param a2 second scale factor
136     * @param ds2 second base (unscaled) derivative structure
137     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
138     */
139    public DerivativeStructure(final double a1, final DerivativeStructure ds1,
140                               final double a2, final DerivativeStructure ds2)
141        throws DimensionMismatchException {
142        this(ds1.compiler);
143        compiler.checkCompatibility(ds2.compiler);
144        compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0, data, 0);
145    }
146
147    /** Linear combination constructor.
148     * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3
149     * @param a1 first scale factor
150     * @param ds1 first base (unscaled) derivative structure
151     * @param a2 second scale factor
152     * @param ds2 second base (unscaled) derivative structure
153     * @param a3 third scale factor
154     * @param ds3 third base (unscaled) derivative structure
155     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
156     */
157    public DerivativeStructure(final double a1, final DerivativeStructure ds1,
158                               final double a2, final DerivativeStructure ds2,
159                               final double a3, final DerivativeStructure ds3)
160        throws DimensionMismatchException {
161        this(ds1.compiler);
162        compiler.checkCompatibility(ds2.compiler);
163        compiler.checkCompatibility(ds3.compiler);
164        compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0, a3, ds3.data, 0, data, 0);
165    }
166
167    /** Linear combination constructor.
168     * The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
169     * @param a1 first scale factor
170     * @param ds1 first base (unscaled) derivative structure
171     * @param a2 second scale factor
172     * @param ds2 second base (unscaled) derivative structure
173     * @param a3 third scale factor
174     * @param ds3 third base (unscaled) derivative structure
175     * @param a4 fourth scale factor
176     * @param ds4 fourth base (unscaled) derivative structure
177     * @exception DimensionMismatchException if number of free parameters or orders are inconsistent
178     */
179    public DerivativeStructure(final double a1, final DerivativeStructure ds1,
180                               final double a2, final DerivativeStructure ds2,
181                               final double a3, final DerivativeStructure ds3,
182                               final double a4, final DerivativeStructure ds4)
183        throws DimensionMismatchException {
184        this(ds1.compiler);
185        compiler.checkCompatibility(ds2.compiler);
186        compiler.checkCompatibility(ds3.compiler);
187        compiler.checkCompatibility(ds4.compiler);
188        compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0,
189                                   a3, ds3.data, 0, a4, ds4.data, 0,
190                                   data, 0);
191    }
192
193    /** Build an instance from all its derivatives.
194     * @param parameters number of free parameters
195     * @param order derivation order
196     * @param derivatives derivatives sorted according to
197     * {@link DSCompiler#getPartialDerivativeIndex(int...)}
198     * @exception DimensionMismatchException if derivatives array does not match the
199     * {@link DSCompiler#getSize() size} expected by the compiler
200     * @throws NumberIsTooLargeException if order is too large
201     * @see #getAllDerivatives()
202     */
203    public DerivativeStructure(final int parameters, final int order, final double ... derivatives)
204        throws DimensionMismatchException, NumberIsTooLargeException {
205        this(parameters, order);
206        if (derivatives.length != data.length) {
207            throw new DimensionMismatchException(derivatives.length, data.length);
208        }
209        System.arraycopy(derivatives, 0, data, 0, data.length);
210    }
211
212    /** Copy constructor.
213     * @param ds instance to copy
214     */
215    private DerivativeStructure(final DerivativeStructure ds) {
216        this.compiler = ds.compiler;
217        this.data     = ds.data.clone();
218    }
219
220    /** Get the number of free parameters.
221     * @return number of free parameters
222     */
223    public int getFreeParameters() {
224        return compiler.getFreeParameters();
225    }
226
227    /** Get the derivation order.
228     * @return derivation order
229     */
230    public int getOrder() {
231        return compiler.getOrder();
232    }
233
234    /** Create a constant compatible with instance order and number of parameters.
235     * <p>
236     * This method is a convenience factory method, it simply calls
237     * {@code new DerivativeStructure(getFreeParameters(), getOrder(), c)}
238     * </p>
239     * @param c value of the constant
240     * @return a constant compatible with instance order and number of parameters
241     * @see #DerivativeStructure(int, int, double)
242     * @since 3.3
243     */
244    public DerivativeStructure createConstant(final double c) {
245        return new DerivativeStructure(getFreeParameters(), getOrder(), c);
246    }
247
248    /** {@inheritDoc}
249     * @since 3.2
250     */
251    public double getReal() {
252        return data[0];
253    }
254
255    /** Get the value part of the derivative structure.
256     * @return value part of the derivative structure
257     * @see #getPartialDerivative(int...)
258     */
259    public double getValue() {
260        return data[0];
261    }
262
263    /** Get a partial derivative.
264     * @param orders derivation orders with respect to each variable (if all orders are 0,
265     * the value is returned)
266     * @return partial derivative
267     * @see #getValue()
268     * @exception DimensionMismatchException if the numbers of variables does not
269     * match the instance
270     * @exception NumberIsTooLargeException if sum of derivation orders is larger
271     * than the instance limits
272     */
273    public double getPartialDerivative(final int ... orders)
274        throws DimensionMismatchException, NumberIsTooLargeException {
275        return data[compiler.getPartialDerivativeIndex(orders)];
276    }
277
278    /** Get all partial derivatives.
279     * @return a fresh copy of partial derivatives, in an array sorted according to
280     * {@link DSCompiler#getPartialDerivativeIndex(int...)}
281     */
282    public double[] getAllDerivatives() {
283        return data.clone();
284    }
285
286    /** {@inheritDoc}
287     * @since 3.2
288     */
289    public DerivativeStructure add(final double a) {
290        final DerivativeStructure ds = new DerivativeStructure(this);
291        ds.data[0] += a;
292        return ds;
293    }
294
295    /** {@inheritDoc}
296     * @exception DimensionMismatchException if number of free parameters
297     * or orders do not match
298     */
299    public DerivativeStructure add(final DerivativeStructure a)
300        throws DimensionMismatchException {
301        compiler.checkCompatibility(a.compiler);
302        final DerivativeStructure ds = new DerivativeStructure(this);
303        compiler.add(data, 0, a.data, 0, ds.data, 0);
304        return ds;
305    }
306
307    /** {@inheritDoc}
308     * @since 3.2
309     */
310    public DerivativeStructure subtract(final double a) {
311        return add(-a);
312    }
313
314    /** {@inheritDoc}
315     * @exception DimensionMismatchException if number of free parameters
316     * or orders do not match
317     */
318    public DerivativeStructure subtract(final DerivativeStructure a)
319        throws DimensionMismatchException {
320        compiler.checkCompatibility(a.compiler);
321        final DerivativeStructure ds = new DerivativeStructure(this);
322        compiler.subtract(data, 0, a.data, 0, ds.data, 0);
323        return ds;
324    }
325
326    /** {@inheritDoc} */
327    public DerivativeStructure multiply(final int n) {
328        return multiply((double) n);
329    }
330
331    /** {@inheritDoc}
332     * @since 3.2
333     */
334    public DerivativeStructure multiply(final double a) {
335        final DerivativeStructure ds = new DerivativeStructure(this);
336        for (int i = 0; i < ds.data.length; ++i) {
337            ds.data[i] *= a;
338        }
339        return ds;
340    }
341
342    /** {@inheritDoc}
343     * @exception DimensionMismatchException if number of free parameters
344     * or orders do not match
345     */
346    public DerivativeStructure multiply(final DerivativeStructure a)
347        throws DimensionMismatchException {
348        compiler.checkCompatibility(a.compiler);
349        final DerivativeStructure result = new DerivativeStructure(compiler);
350        compiler.multiply(data, 0, a.data, 0, result.data, 0);
351        return result;
352    }
353
354    /** {@inheritDoc}
355     * @since 3.2
356     */
357    public DerivativeStructure divide(final double a) {
358        final DerivativeStructure ds = new DerivativeStructure(this);
359        for (int i = 0; i < ds.data.length; ++i) {
360            ds.data[i] /= a;
361        }
362        return ds;
363    }
364
365    /** {@inheritDoc}
366     * @exception DimensionMismatchException if number of free parameters
367     * or orders do not match
368     */
369    public DerivativeStructure divide(final DerivativeStructure a)
370        throws DimensionMismatchException {
371        compiler.checkCompatibility(a.compiler);
372        final DerivativeStructure result = new DerivativeStructure(compiler);
373        compiler.divide(data, 0, a.data, 0, result.data, 0);
374        return result;
375    }
376
377    /** {@inheritDoc} */
378    public DerivativeStructure remainder(final double a) {
379        final DerivativeStructure ds = new DerivativeStructure(this);
380        ds.data[0] = FastMath.IEEEremainder(ds.data[0], a);
381        return ds;
382    }
383
384    /** {@inheritDoc}
385     * @exception DimensionMismatchException if number of free parameters
386     * or orders do not match
387     * @since 3.2
388     */
389    public DerivativeStructure remainder(final DerivativeStructure a)
390        throws DimensionMismatchException {
391        compiler.checkCompatibility(a.compiler);
392        final DerivativeStructure result = new DerivativeStructure(compiler);
393        compiler.remainder(data, 0, a.data, 0, result.data, 0);
394        return result;
395    }
396
397    /** {@inheritDoc} */
398    public DerivativeStructure negate() {
399        final DerivativeStructure ds = new DerivativeStructure(compiler);
400        for (int i = 0; i < ds.data.length; ++i) {
401            ds.data[i] = -data[i];
402        }
403        return ds;
404    }
405
406    /** {@inheritDoc}
407     * @since 3.2
408     */
409    public DerivativeStructure abs() {
410        if (Double.doubleToLongBits(data[0]) < 0) {
411            // we use the bits representation to also handle -0.0
412            return negate();
413        } else {
414            return this;
415        }
416    }
417
418    /** {@inheritDoc}
419     * @since 3.2
420     */
421    public DerivativeStructure ceil() {
422        return new DerivativeStructure(compiler.getFreeParameters(),
423                                       compiler.getOrder(),
424                                       FastMath.ceil(data[0]));
425    }
426
427    /** {@inheritDoc}
428     * @since 3.2
429     */
430    public DerivativeStructure floor() {
431        return new DerivativeStructure(compiler.getFreeParameters(),
432                                       compiler.getOrder(),
433                                       FastMath.floor(data[0]));
434    }
435
436    /** {@inheritDoc}
437     * @since 3.2
438     */
439    public DerivativeStructure rint() {
440        return new DerivativeStructure(compiler.getFreeParameters(),
441                                       compiler.getOrder(),
442                                       FastMath.rint(data[0]));
443    }
444
445    /** {@inheritDoc} */
446    public long round() {
447        return FastMath.round(data[0]);
448    }
449
450    /** {@inheritDoc}
451     * @since 3.2
452     */
453    public DerivativeStructure signum() {
454        return new DerivativeStructure(compiler.getFreeParameters(),
455                                       compiler.getOrder(),
456                                       FastMath.signum(data[0]));
457    }
458
459    /** {@inheritDoc}
460     * @since 3.2
461     */
462    public DerivativeStructure copySign(final DerivativeStructure sign){
463        long m = Double.doubleToLongBits(data[0]);
464        long s = Double.doubleToLongBits(sign.data[0]);
465        if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
466            return this;
467        }
468        return negate(); // flip sign
469    }
470
471    /** {@inheritDoc}
472     * @since 3.2
473     */
474    public DerivativeStructure copySign(final double sign) {
475        long m = Double.doubleToLongBits(data[0]);
476        long s = Double.doubleToLongBits(sign);
477        if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
478            return this;
479        }
480        return negate(); // flip sign
481    }
482
483    /**
484     * Return the exponent of the instance value, removing the bias.
485     * <p>
486     * For double numbers of the form 2<sup>x</sup>, the unbiased
487     * exponent is exactly x.
488     * </p>
489     * @return exponent for instance in IEEE754 representation, without bias
490     */
491    public int getExponent() {
492        return FastMath.getExponent(data[0]);
493    }
494
495    /** {@inheritDoc}
496     * @since 3.2
497     */
498    public DerivativeStructure scalb(final int n) {
499        final DerivativeStructure ds = new DerivativeStructure(compiler);
500        for (int i = 0; i < ds.data.length; ++i) {
501            ds.data[i] = FastMath.scalb(data[i], n);
502        }
503        return ds;
504    }
505
506    /** {@inheritDoc}
507     * @exception DimensionMismatchException if number of free parameters
508     * or orders do not match
509     * @since 3.2
510     */
511    public DerivativeStructure hypot(final DerivativeStructure y)
512        throws DimensionMismatchException {
513
514        compiler.checkCompatibility(y.compiler);
515
516        if (Double.isInfinite(data[0]) || Double.isInfinite(y.data[0])) {
517            return new DerivativeStructure(compiler.getFreeParameters(),
518                                           compiler.getFreeParameters(),
519                                           Double.POSITIVE_INFINITY);
520        } else if (Double.isNaN(data[0]) || Double.isNaN(y.data[0])) {
521            return new DerivativeStructure(compiler.getFreeParameters(),
522                                           compiler.getFreeParameters(),
523                                           Double.NaN);
524        } else {
525
526            final int expX = getExponent();
527            final int expY = y.getExponent();
528            if (expX > expY + 27) {
529                // y is neglectible with respect to x
530                return abs();
531            } else if (expY > expX + 27) {
532                // x is neglectible with respect to y
533                return y.abs();
534            } else {
535
536                // find an intermediate scale to avoid both overflow and underflow
537                final int middleExp = (expX + expY) / 2;
538
539                // scale parameters without losing precision
540                final DerivativeStructure scaledX = scalb(-middleExp);
541                final DerivativeStructure scaledY = y.scalb(-middleExp);
542
543                // compute scaled hypotenuse
544                final DerivativeStructure scaledH =
545                        scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
546
547                // remove scaling
548                return scaledH.scalb(middleExp);
549
550            }
551
552        }
553    }
554
555    /**
556     * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
557     * - sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)<br/>
558     * avoiding intermediate overflow or underflow.
559     *
560     * <ul>
561     * <li> If either argument is infinite, then the result is positive infinity.</li>
562     * <li> else, if either argument is NaN then the result is NaN.</li>
563     * </ul>
564     *
565     * @param x a value
566     * @param y a value
567     * @return sqrt(<i>x</i><sup>2</sup>&nbsp;+<i>y</i><sup>2</sup>)
568     * @exception DimensionMismatchException if number of free parameters
569     * or orders do not match
570     * @since 3.2
571     */
572    public static DerivativeStructure hypot(final DerivativeStructure x, final DerivativeStructure y)
573        throws DimensionMismatchException {
574        return x.hypot(y);
575    }
576
577    /** Compute composition of the instance by a univariate function.
578     * @param f array of value and derivatives of the function at
579     * the current point (i.e. [f({@link #getValue()}),
580     * f'({@link #getValue()}), f''({@link #getValue()})...]).
581     * @return f(this)
582     * @exception DimensionMismatchException if the number of derivatives
583     * in the array is not equal to {@link #getOrder() order} + 1
584     */
585    public DerivativeStructure compose(final double ... f)
586        throws DimensionMismatchException {
587        if (f.length != getOrder() + 1) {
588            throw new DimensionMismatchException(f.length, getOrder() + 1);
589        }
590        final DerivativeStructure result = new DerivativeStructure(compiler);
591        compiler.compose(data, 0, f, result.data, 0);
592        return result;
593    }
594
595    /** {@inheritDoc} */
596    public DerivativeStructure reciprocal() {
597        final DerivativeStructure result = new DerivativeStructure(compiler);
598        compiler.pow(data, 0, -1, result.data, 0);
599        return result;
600    }
601
602    /** {@inheritDoc}
603     * @since 3.2
604     */
605    public DerivativeStructure sqrt() {
606        return rootN(2);
607    }
608
609    /** {@inheritDoc}
610     * @since 3.2
611     */
612    public DerivativeStructure cbrt() {
613        return rootN(3);
614    }
615
616    /** {@inheritDoc}
617     * @since 3.2
618     */
619    public DerivativeStructure rootN(final int n) {
620        final DerivativeStructure result = new DerivativeStructure(compiler);
621        compiler.rootN(data, 0, n, result.data, 0);
622        return result;
623    }
624
625    /** {@inheritDoc} */
626    public Field<DerivativeStructure> getField() {
627        return new Field<DerivativeStructure>() {
628
629            /** {@inheritDoc} */
630            public DerivativeStructure getZero() {
631                return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 0.0);
632            }
633
634            /** {@inheritDoc} */
635            public DerivativeStructure getOne() {
636                return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 1.0);
637            }
638
639            /** {@inheritDoc} */
640            public Class<? extends FieldElement<DerivativeStructure>> getRuntimeClass() {
641                return DerivativeStructure.class;
642            }
643
644        };
645    }
646
647    /** Compute a<sup>x</sup> where a is a double and x a {@link DerivativeStructure}
648     * @param a number to exponentiate
649     * @param x power to apply
650     * @return a<sup>x</sup>
651     * @since 3.3
652     */
653    public static DerivativeStructure pow(final double a, final DerivativeStructure x) {
654        final DerivativeStructure result = new DerivativeStructure(x.compiler);
655        x.compiler.pow(a, x.data, 0, result.data, 0);
656        return result;
657    }
658
659    /** {@inheritDoc}
660     * @since 3.2
661     */
662    public DerivativeStructure pow(final double p) {
663        final DerivativeStructure result = new DerivativeStructure(compiler);
664        compiler.pow(data, 0, p, result.data, 0);
665        return result;
666    }
667
668    /** {@inheritDoc}
669     * @since 3.2
670     */
671    public DerivativeStructure pow(final int n) {
672        final DerivativeStructure result = new DerivativeStructure(compiler);
673        compiler.pow(data, 0, n, result.data, 0);
674        return result;
675    }
676
677    /** {@inheritDoc}
678     * @exception DimensionMismatchException if number of free parameters
679     * or orders do not match
680     * @since 3.2
681     */
682    public DerivativeStructure pow(final DerivativeStructure e)
683        throws DimensionMismatchException {
684        compiler.checkCompatibility(e.compiler);
685        final DerivativeStructure result = new DerivativeStructure(compiler);
686        compiler.pow(data, 0, e.data, 0, result.data, 0);
687        return result;
688    }
689
690    /** {@inheritDoc}
691     * @since 3.2
692     */
693    public DerivativeStructure exp() {
694        final DerivativeStructure result = new DerivativeStructure(compiler);
695        compiler.exp(data, 0, result.data, 0);
696        return result;
697    }
698
699    /** {@inheritDoc}
700     * @since 3.2
701     */
702    public DerivativeStructure expm1() {
703        final DerivativeStructure result = new DerivativeStructure(compiler);
704        compiler.expm1(data, 0, result.data, 0);
705        return result;
706    }
707
708    /** {@inheritDoc}
709     * @since 3.2
710     */
711    public DerivativeStructure log() {
712        final DerivativeStructure result = new DerivativeStructure(compiler);
713        compiler.log(data, 0, result.data, 0);
714        return result;
715    }
716
717    /** {@inheritDoc}
718     * @since 3.2
719     */
720    public DerivativeStructure log1p() {
721        final DerivativeStructure result = new DerivativeStructure(compiler);
722        compiler.log1p(data, 0, result.data, 0);
723        return result;
724    }
725
726    /** Base 10 logarithm.
727     * @return base 10 logarithm of the instance
728     */
729    public DerivativeStructure log10() {
730        final DerivativeStructure result = new DerivativeStructure(compiler);
731        compiler.log10(data, 0, result.data, 0);
732        return result;
733    }
734
735    /** {@inheritDoc}
736     * @since 3.2
737     */
738    public DerivativeStructure cos() {
739        final DerivativeStructure result = new DerivativeStructure(compiler);
740        compiler.cos(data, 0, result.data, 0);
741        return result;
742    }
743
744    /** {@inheritDoc}
745     * @since 3.2
746     */
747    public DerivativeStructure sin() {
748        final DerivativeStructure result = new DerivativeStructure(compiler);
749        compiler.sin(data, 0, result.data, 0);
750        return result;
751    }
752
753    /** {@inheritDoc}
754     * @since 3.2
755     */
756    public DerivativeStructure tan() {
757        final DerivativeStructure result = new DerivativeStructure(compiler);
758        compiler.tan(data, 0, result.data, 0);
759        return result;
760    }
761
762    /** {@inheritDoc}
763     * @since 3.2
764     */
765    public DerivativeStructure acos() {
766        final DerivativeStructure result = new DerivativeStructure(compiler);
767        compiler.acos(data, 0, result.data, 0);
768        return result;
769    }
770
771    /** {@inheritDoc}
772     * @since 3.2
773     */
774    public DerivativeStructure asin() {
775        final DerivativeStructure result = new DerivativeStructure(compiler);
776        compiler.asin(data, 0, result.data, 0);
777        return result;
778    }
779
780    /** {@inheritDoc}
781     * @since 3.2
782     */
783    public DerivativeStructure atan() {
784        final DerivativeStructure result = new DerivativeStructure(compiler);
785        compiler.atan(data, 0, result.data, 0);
786        return result;
787    }
788
789    /** {@inheritDoc}
790     * @since 3.2
791     */
792    public DerivativeStructure atan2(final DerivativeStructure x)
793        throws DimensionMismatchException {
794        compiler.checkCompatibility(x.compiler);
795        final DerivativeStructure result = new DerivativeStructure(compiler);
796        compiler.atan2(data, 0, x.data, 0, result.data, 0);
797        return result;
798    }
799
800    /** Two arguments arc tangent operation.
801     * @param y first argument of the arc tangent
802     * @param x second argument of the arc tangent
803     * @return atan2(y, x)
804     * @exception DimensionMismatchException if number of free parameters
805     * or orders do not match
806     * @since 3.2
807     */
808    public static DerivativeStructure atan2(final DerivativeStructure y, final DerivativeStructure x)
809        throws DimensionMismatchException {
810        return y.atan2(x);
811    }
812
813    /** {@inheritDoc}
814     * @since 3.2
815     */
816    public DerivativeStructure cosh() {
817        final DerivativeStructure result = new DerivativeStructure(compiler);
818        compiler.cosh(data, 0, result.data, 0);
819        return result;
820    }
821
822    /** {@inheritDoc}
823     * @since 3.2
824     */
825    public DerivativeStructure sinh() {
826        final DerivativeStructure result = new DerivativeStructure(compiler);
827        compiler.sinh(data, 0, result.data, 0);
828        return result;
829    }
830
831    /** {@inheritDoc}
832     * @since 3.2
833     */
834    public DerivativeStructure tanh() {
835        final DerivativeStructure result = new DerivativeStructure(compiler);
836        compiler.tanh(data, 0, result.data, 0);
837        return result;
838    }
839
840    /** {@inheritDoc}
841     * @since 3.2
842     */
843    public DerivativeStructure acosh() {
844        final DerivativeStructure result = new DerivativeStructure(compiler);
845        compiler.acosh(data, 0, result.data, 0);
846        return result;
847    }
848
849    /** {@inheritDoc}
850     * @since 3.2
851     */
852    public DerivativeStructure asinh() {
853        final DerivativeStructure result = new DerivativeStructure(compiler);
854        compiler.asinh(data, 0, result.data, 0);
855        return result;
856    }
857
858    /** {@inheritDoc}
859     * @since 3.2
860     */
861    public DerivativeStructure atanh() {
862        final DerivativeStructure result = new DerivativeStructure(compiler);
863        compiler.atanh(data, 0, result.data, 0);
864        return result;
865    }
866
867    /** Convert radians to degrees, with error of less than 0.5 ULP
868     *  @return instance converted into degrees
869     */
870    public DerivativeStructure toDegrees() {
871        final DerivativeStructure ds = new DerivativeStructure(compiler);
872        for (int i = 0; i < ds.data.length; ++i) {
873            ds.data[i] = FastMath.toDegrees(data[i]);
874        }
875        return ds;
876    }
877
878    /** Convert degrees to radians, with error of less than 0.5 ULP
879     *  @return instance converted into radians
880     */
881    public DerivativeStructure toRadians() {
882        final DerivativeStructure ds = new DerivativeStructure(compiler);
883        for (int i = 0; i < ds.data.length; ++i) {
884            ds.data[i] = FastMath.toRadians(data[i]);
885        }
886        return ds;
887    }
888
889    /** Evaluate Taylor expansion a derivative structure.
890     * @param delta parameters offsets (&Delta;x, &Delta;y, ...)
891     * @return value of the Taylor expansion at x + &Delta;x, y + &Delta;y, ...
892     * @throws MathArithmeticException if factorials becomes too large
893     */
894    public double taylor(final double ... delta) throws MathArithmeticException {
895        return compiler.taylor(data, 0, delta);
896    }
897
898    /** {@inheritDoc}
899     * @exception DimensionMismatchException if number of free parameters
900     * or orders do not match
901     * @since 3.2
902     */
903    public DerivativeStructure linearCombination(final DerivativeStructure[] a, final DerivativeStructure[] b)
904        throws DimensionMismatchException {
905
906        // compute an accurate value, taking care of cancellations
907        final double[] aDouble = new double[a.length];
908        for (int i = 0; i < a.length; ++i) {
909            aDouble[i] = a[i].getValue();
910        }
911        final double[] bDouble = new double[b.length];
912        for (int i = 0; i < b.length; ++i) {
913            bDouble[i] = b[i].getValue();
914        }
915        final double accurateValue = MathArrays.linearCombination(aDouble, bDouble);
916
917        // compute a simple value, with all partial derivatives
918        DerivativeStructure simpleValue = a[0].getField().getZero();
919        for (int i = 0; i < a.length; ++i) {
920            simpleValue = simpleValue.add(a[i].multiply(b[i]));
921        }
922
923        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
924        final double[] all = simpleValue.getAllDerivatives();
925        all[0] = accurateValue;
926        return new DerivativeStructure(simpleValue.getFreeParameters(), simpleValue.getOrder(), all);
927
928    }
929
930    /** {@inheritDoc}
931     * @exception DimensionMismatchException if number of free parameters
932     * or orders do not match
933     * @since 3.2
934     */
935    public DerivativeStructure linearCombination(final double[] a, final DerivativeStructure[] b)
936        throws DimensionMismatchException {
937
938        // compute an accurate value, taking care of cancellations
939        final double[] bDouble = new double[b.length];
940        for (int i = 0; i < b.length; ++i) {
941            bDouble[i] = b[i].getValue();
942        }
943        final double accurateValue = MathArrays.linearCombination(a, bDouble);
944
945        // compute a simple value, with all partial derivatives
946        DerivativeStructure simpleValue = b[0].getField().getZero();
947        for (int i = 0; i < a.length; ++i) {
948            simpleValue = simpleValue.add(b[i].multiply(a[i]));
949        }
950
951        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
952        final double[] all = simpleValue.getAllDerivatives();
953        all[0] = accurateValue;
954        return new DerivativeStructure(simpleValue.getFreeParameters(), simpleValue.getOrder(), all);
955
956    }
957
958    /** {@inheritDoc}
959     * @exception DimensionMismatchException if number of free parameters
960     * or orders do not match
961     * @since 3.2
962     */
963    public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
964                                                 final DerivativeStructure a2, final DerivativeStructure b2)
965        throws DimensionMismatchException {
966
967        // compute an accurate value, taking care of cancellations
968        final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
969                                                                  a2.getValue(), b2.getValue());
970
971        // compute a simple value, with all partial derivatives
972        final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2));
973
974        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
975        final double[] all = simpleValue.getAllDerivatives();
976        all[0] = accurateValue;
977        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
978
979    }
980
981    /** {@inheritDoc}
982     * @exception DimensionMismatchException if number of free parameters
983     * or orders do not match
984     * @since 3.2
985     */
986    public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
987                                                 final double a2, final DerivativeStructure b2)
988        throws DimensionMismatchException {
989
990        // compute an accurate value, taking care of cancellations
991        final double accurateValue = MathArrays.linearCombination(a1, b1.getValue(),
992                                                                  a2, b2.getValue());
993
994        // compute a simple value, with all partial derivatives
995        final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2));
996
997        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
998        final double[] all = simpleValue.getAllDerivatives();
999        all[0] = accurateValue;
1000        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1001
1002    }
1003
1004    /** {@inheritDoc}
1005     * @exception DimensionMismatchException if number of free parameters
1006     * or orders do not match
1007     * @since 3.2
1008     */
1009    public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
1010                                                 final DerivativeStructure a2, final DerivativeStructure b2,
1011                                                 final DerivativeStructure a3, final DerivativeStructure b3)
1012        throws DimensionMismatchException {
1013
1014        // compute an accurate value, taking care of cancellations
1015        final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
1016                                                                  a2.getValue(), b2.getValue(),
1017                                                                  a3.getValue(), b3.getValue());
1018
1019        // compute a simple value, with all partial derivatives
1020        final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
1021
1022        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1023        final double[] all = simpleValue.getAllDerivatives();
1024        all[0] = accurateValue;
1025        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1026
1027    }
1028
1029    /** {@inheritDoc}
1030     * @exception DimensionMismatchException if number of free parameters
1031     * or orders do not match
1032     * @since 3.2
1033     */
1034    public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1035                                                 final double a2, final DerivativeStructure b2,
1036                                                 final double a3, final DerivativeStructure b3)
1037        throws DimensionMismatchException {
1038
1039        // compute an accurate value, taking care of cancellations
1040        final double accurateValue = MathArrays.linearCombination(a1, b1.getValue(),
1041                                                                  a2, b2.getValue(),
1042                                                                  a3, b3.getValue());
1043
1044        // compute a simple value, with all partial derivatives
1045        final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3));
1046
1047        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1048        final double[] all = simpleValue.getAllDerivatives();
1049        all[0] = accurateValue;
1050        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1051
1052    }
1053
1054    /** {@inheritDoc}
1055     * @exception DimensionMismatchException if number of free parameters
1056     * or orders do not match
1057     * @since 3.2
1058     */
1059    public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
1060                                                 final DerivativeStructure a2, final DerivativeStructure b2,
1061                                                 final DerivativeStructure a3, final DerivativeStructure b3,
1062                                                 final DerivativeStructure a4, final DerivativeStructure b4)
1063        throws DimensionMismatchException {
1064
1065        // compute an accurate value, taking care of cancellations
1066        final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
1067                                                                  a2.getValue(), b2.getValue(),
1068                                                                  a3.getValue(), b3.getValue(),
1069                                                                  a4.getValue(), b4.getValue());
1070
1071        // compute a simple value, with all partial derivatives
1072        final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
1073
1074        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1075        final double[] all = simpleValue.getAllDerivatives();
1076        all[0] = accurateValue;
1077        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1078
1079    }
1080
1081    /** {@inheritDoc}
1082     * @exception DimensionMismatchException if number of free parameters
1083     * or orders do not match
1084     * @since 3.2
1085     */
1086    public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1087                                                 final double a2, final DerivativeStructure b2,
1088                                                 final double a3, final DerivativeStructure b3,
1089                                                 final double a4, final DerivativeStructure b4)
1090        throws DimensionMismatchException {
1091
1092        // compute an accurate value, taking care of cancellations
1093        final double accurateValue = MathArrays.linearCombination(a1, b1.getValue(),
1094                                                                  a2, b2.getValue(),
1095                                                                  a3, b3.getValue(),
1096                                                                  a4, b4.getValue());
1097
1098        // compute a simple value, with all partial derivatives
1099        final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4));
1100
1101        // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1102        final double[] all = simpleValue.getAllDerivatives();
1103        all[0] = accurateValue;
1104        return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1105
1106    }
1107
1108    /**
1109     * Test for the equality of two derivative structures.
1110     * <p>
1111     * Derivative structures are considered equal if they have the same number
1112     * of free parameters, the same derivation order, and the same derivatives.
1113     * </p>
1114     * @param other Object to test for equality to this
1115     * @return true if two derivative structures are equal
1116     * @since 3.2
1117     */
1118    @Override
1119    public boolean equals(Object other) {
1120
1121        if (this == other) {
1122            return true;
1123        }
1124
1125        if (other instanceof DerivativeStructure) {
1126            final DerivativeStructure rhs = (DerivativeStructure)other;
1127            return (getFreeParameters() == rhs.getFreeParameters()) &&
1128                   (getOrder() == rhs.getOrder()) &&
1129                   MathArrays.equals(data, rhs.data);
1130        }
1131
1132        return false;
1133
1134    }
1135
1136    /**
1137     * Get a hashCode for the derivative structure.
1138     * @return a hash code value for this object
1139     * @since 3.2
1140     */
1141    @Override
1142    public int hashCode() {
1143        return 227 + 229 * getFreeParameters() + 233 * getOrder() + 239 * MathUtils.hash(data);
1144    }
1145
1146    /**
1147     * Replace the instance with a data transfer object for serialization.
1148     * @return data transfer object that will be serialized
1149     */
1150    private Object writeReplace() {
1151        return new DataTransferObject(compiler.getFreeParameters(), compiler.getOrder(), data);
1152    }
1153
1154    /** Internal class used only for serialization. */
1155    private static class DataTransferObject implements Serializable {
1156
1157        /** Serializable UID. */
1158        private static final long serialVersionUID = 20120730L;
1159
1160        /** Number of variables.
1161         * @serial
1162         */
1163        private final int variables;
1164
1165        /** Derivation order.
1166         * @serial
1167         */
1168        private final int order;
1169
1170        /** Partial derivatives.
1171         * @serial
1172         */
1173        private final double[] data;
1174
1175        /** Simple constructor.
1176         * @param variables number of variables
1177         * @param order derivation order
1178         * @param data partial derivatives
1179         */
1180        public DataTransferObject(final int variables, final int order, final double[] data) {
1181            this.variables = variables;
1182            this.order     = order;
1183            this.data      = data;
1184        }
1185
1186        /** Replace the deserialized data transfer object with a {@link DerivativeStructure}.
1187         * @return replacement {@link DerivativeStructure}
1188         */
1189        private Object readResolve() {
1190            return new DerivativeStructure(variables, order, data);
1191        }
1192
1193    }
1194
1195}