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