001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.analysis.differentiation;
018
019 import java.io.Serializable;
020
021 import org.apache.commons.math3.RealFieldElement;
022 import org.apache.commons.math3.Field;
023 import org.apache.commons.math3.FieldElement;
024 import org.apache.commons.math3.exception.DimensionMismatchException;
025 import org.apache.commons.math3.exception.MathArithmeticException;
026 import org.apache.commons.math3.exception.NumberIsTooLargeException;
027 import org.apache.commons.math3.util.FastMath;
028 import org.apache.commons.math3.util.MathArrays;
029 import 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://www.math.american.edu/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 static
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 1462423 2013-03-29 07:25:18Z luc $
060 * @since 3.1
061 */
062 public 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 /** {@inheritDoc}
236 * @since 3.2
237 */
238 public double getReal() {
239 return data[0];
240 }
241
242 /** Get the value part of the derivative structure.
243 * @return value part of the derivative structure
244 * @see #getPartialDerivative(int...)
245 */
246 public double getValue() {
247 return data[0];
248 }
249
250 /** Get a partial derivative.
251 * @param orders derivation orders with respect to each variable (if all orders are 0,
252 * the value is returned)
253 * @return partial derivative
254 * @see #getValue()
255 * @exception DimensionMismatchException if the numbers of variables does not
256 * match the instance
257 * @exception NumberIsTooLargeException if sum of derivation orders is larger
258 * than the instance limits
259 */
260 public double getPartialDerivative(final int ... orders)
261 throws DimensionMismatchException, NumberIsTooLargeException {
262 return data[compiler.getPartialDerivativeIndex(orders)];
263 }
264
265 /** Get all partial derivatives.
266 * @return a fresh copy of partial derivatives, in an array sorted according to
267 * {@link DSCompiler#getPartialDerivativeIndex(int...)}
268 */
269 public double[] getAllDerivatives() {
270 return data.clone();
271 }
272
273 /** {@inheritDoc}
274 * @since 3.2
275 */
276 public DerivativeStructure add(final double a) {
277 final DerivativeStructure ds = new DerivativeStructure(this);
278 ds.data[0] += a;
279 return ds;
280 }
281
282 /** {@inheritDoc}
283 * @exception DimensionMismatchException if number of free parameters
284 * or orders do not match
285 */
286 public DerivativeStructure add(final DerivativeStructure a)
287 throws DimensionMismatchException {
288 compiler.checkCompatibility(a.compiler);
289 final DerivativeStructure ds = new DerivativeStructure(this);
290 compiler.add(data, 0, a.data, 0, ds.data, 0);
291 return ds;
292 }
293
294 /** {@inheritDoc}
295 * @since 3.2
296 */
297 public DerivativeStructure subtract(final double a) {
298 return add(-a);
299 }
300
301 /** {@inheritDoc}
302 * @exception DimensionMismatchException if number of free parameters
303 * or orders do not match
304 */
305 public DerivativeStructure subtract(final DerivativeStructure a)
306 throws DimensionMismatchException {
307 compiler.checkCompatibility(a.compiler);
308 final DerivativeStructure ds = new DerivativeStructure(this);
309 compiler.subtract(data, 0, a.data, 0, ds.data, 0);
310 return ds;
311 }
312
313 /** {@inheritDoc} */
314 public DerivativeStructure multiply(final int n) {
315 return multiply((double) n);
316 }
317
318 /** {@inheritDoc}
319 * @since 3.2
320 */
321 public DerivativeStructure multiply(final double a) {
322 final DerivativeStructure ds = new DerivativeStructure(this);
323 for (int i = 0; i < ds.data.length; ++i) {
324 ds.data[i] *= a;
325 }
326 return ds;
327 }
328
329 /** {@inheritDoc}
330 * @exception DimensionMismatchException if number of free parameters
331 * or orders do not match
332 */
333 public DerivativeStructure multiply(final DerivativeStructure a)
334 throws DimensionMismatchException {
335 compiler.checkCompatibility(a.compiler);
336 final DerivativeStructure result = new DerivativeStructure(compiler);
337 compiler.multiply(data, 0, a.data, 0, result.data, 0);
338 return result;
339 }
340
341 /** {@inheritDoc}
342 * @since 3.2
343 */
344 public DerivativeStructure divide(final double a) {
345 final DerivativeStructure ds = new DerivativeStructure(this);
346 for (int i = 0; i < ds.data.length; ++i) {
347 ds.data[i] /= a;
348 }
349 return ds;
350 }
351
352 /** {@inheritDoc}
353 * @exception DimensionMismatchException if number of free parameters
354 * or orders do not match
355 */
356 public DerivativeStructure divide(final DerivativeStructure a)
357 throws DimensionMismatchException {
358 compiler.checkCompatibility(a.compiler);
359 final DerivativeStructure result = new DerivativeStructure(compiler);
360 compiler.divide(data, 0, a.data, 0, result.data, 0);
361 return result;
362 }
363
364 /** {@inheritDoc} */
365 public DerivativeStructure remainder(final double a) {
366 final DerivativeStructure ds = new DerivativeStructure(this);
367 ds.data[0] = FastMath.IEEEremainder(ds.data[0], a);
368 return ds;
369 }
370
371 /** {@inheritDoc}
372 * @exception DimensionMismatchException if number of free parameters
373 * or orders do not match
374 * @since 3.2
375 */
376 public DerivativeStructure remainder(final DerivativeStructure a)
377 throws DimensionMismatchException {
378 compiler.checkCompatibility(a.compiler);
379 final DerivativeStructure result = new DerivativeStructure(compiler);
380 compiler.remainder(data, 0, a.data, 0, result.data, 0);
381 return result;
382 }
383
384 /** {@inheritDoc} */
385 public DerivativeStructure negate() {
386 final DerivativeStructure ds = new DerivativeStructure(compiler);
387 for (int i = 0; i < ds.data.length; ++i) {
388 ds.data[i] = -data[i];
389 }
390 return ds;
391 }
392
393 /** {@inheritDoc}
394 * @since 3.2
395 */
396 public DerivativeStructure abs() {
397 if (Double.doubleToLongBits(data[0]) < 0) {
398 // we use the bits representation to also handle -0.0
399 return negate();
400 } else {
401 return this;
402 }
403 }
404
405 /** {@inheritDoc}
406 * @since 3.2
407 */
408 public DerivativeStructure ceil() {
409 return new DerivativeStructure(compiler.getFreeParameters(),
410 compiler.getOrder(),
411 FastMath.ceil(data[0]));
412 }
413
414 /** {@inheritDoc}
415 * @since 3.2
416 */
417 public DerivativeStructure floor() {
418 return new DerivativeStructure(compiler.getFreeParameters(),
419 compiler.getOrder(),
420 FastMath.floor(data[0]));
421 }
422
423 /** {@inheritDoc}
424 * @since 3.2
425 */
426 public DerivativeStructure rint() {
427 return new DerivativeStructure(compiler.getFreeParameters(),
428 compiler.getOrder(),
429 FastMath.rint(data[0]));
430 }
431
432 /** {@inheritDoc} */
433 public long round() {
434 return FastMath.round(data[0]);
435 }
436
437 /** {@inheritDoc}
438 * @since 3.2
439 */
440 public DerivativeStructure signum() {
441 return new DerivativeStructure(compiler.getFreeParameters(),
442 compiler.getOrder(),
443 FastMath.signum(data[0]));
444 }
445
446 /** {@inheritDoc}
447 * @since 3.2
448 */
449 public DerivativeStructure copySign(final DerivativeStructure sign){
450 long m = Double.doubleToLongBits(data[0]);
451 long s = Double.doubleToLongBits(sign.data[0]);
452 if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
453 return this;
454 }
455 return negate(); // flip sign
456 }
457
458 /** {@inheritDoc}
459 * @since 3.2
460 */
461 public DerivativeStructure copySign(final double sign) {
462 long m = Double.doubleToLongBits(data[0]);
463 long s = Double.doubleToLongBits(sign);
464 if ((m >= 0 && s >= 0) || (m < 0 && s < 0)) { // Sign is currently OK
465 return this;
466 }
467 return negate(); // flip sign
468 }
469
470 /**
471 * Return the exponent of the instance value, removing the bias.
472 * <p>
473 * For double numbers of the form 2<sup>x</sup>, the unbiased
474 * exponent is exactly x.
475 * </p>
476 * @return exponent for instance in IEEE754 representation, without bias
477 */
478 public int getExponent() {
479 return FastMath.getExponent(data[0]);
480 }
481
482 /** {@inheritDoc}
483 * @since 3.2
484 */
485 public DerivativeStructure scalb(final int n) {
486 final DerivativeStructure ds = new DerivativeStructure(compiler);
487 for (int i = 0; i < ds.data.length; ++i) {
488 ds.data[i] = FastMath.scalb(data[i], n);
489 }
490 return ds;
491 }
492
493 /** {@inheritDoc}
494 * @exception DimensionMismatchException if number of free parameters
495 * or orders do not match
496 * @since 3.2
497 */
498 public DerivativeStructure hypot(final DerivativeStructure y)
499 throws DimensionMismatchException {
500
501 compiler.checkCompatibility(y.compiler);
502
503 if (Double.isInfinite(data[0]) || Double.isInfinite(y.data[0])) {
504 return new DerivativeStructure(compiler.getFreeParameters(),
505 compiler.getFreeParameters(),
506 Double.POSITIVE_INFINITY);
507 } else if (Double.isNaN(data[0]) || Double.isNaN(y.data[0])) {
508 return new DerivativeStructure(compiler.getFreeParameters(),
509 compiler.getFreeParameters(),
510 Double.NaN);
511 } else {
512
513 final int expX = getExponent();
514 final int expY = y.getExponent();
515 if (expX > expY + 27) {
516 // y is neglectible with respect to x
517 return abs();
518 } else if (expY > expX + 27) {
519 // x is neglectible with respect to y
520 return y.abs();
521 } else {
522
523 // find an intermediate scale to avoid both overflow and underflow
524 final int middleExp = (expX + expY) / 2;
525
526 // scale parameters without losing precision
527 final DerivativeStructure scaledX = scalb(-middleExp);
528 final DerivativeStructure scaledY = y.scalb(-middleExp);
529
530 // compute scaled hypotenuse
531 final DerivativeStructure scaledH =
532 scaledX.multiply(scaledX).add(scaledY.multiply(scaledY)).sqrt();
533
534 // remove scaling
535 return scaledH.scalb(middleExp);
536
537 }
538
539 }
540 }
541
542 /**
543 * Returns the hypotenuse of a triangle with sides {@code x} and {@code y}
544 * - sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)<br/>
545 * avoiding intermediate overflow or underflow.
546 *
547 * <ul>
548 * <li> If either argument is infinite, then the result is positive infinity.</li>
549 * <li> else, if either argument is NaN then the result is NaN.</li>
550 * </ul>
551 *
552 * @param x a value
553 * @param y a value
554 * @return sqrt(<i>x</i><sup>2</sup> +<i>y</i><sup>2</sup>)
555 * @exception DimensionMismatchException if number of free parameters
556 * or orders do not match
557 * @since 3.2
558 */
559 public static DerivativeStructure hypot(final DerivativeStructure x, final DerivativeStructure y)
560 throws DimensionMismatchException {
561 return x.hypot(y);
562 }
563
564 /** Compute composition of the instance by a univariate function.
565 * @param f array of value and derivatives of the function at
566 * the current point (i.e. [f({@link #getValue()}),
567 * f'({@link #getValue()}), f''({@link #getValue()})...]).
568 * @return f(this)
569 * @exception DimensionMismatchException if the number of derivatives
570 * in the array is not equal to {@link #getOrder() order} + 1
571 */
572 public DerivativeStructure compose(final double ... f)
573 throws DimensionMismatchException {
574 if (f.length != getOrder() + 1) {
575 throw new DimensionMismatchException(f.length, getOrder() + 1);
576 }
577 final DerivativeStructure result = new DerivativeStructure(compiler);
578 compiler.compose(data, 0, f, result.data, 0);
579 return result;
580 }
581
582 /** {@inheritDoc} */
583 public DerivativeStructure reciprocal() {
584 final DerivativeStructure result = new DerivativeStructure(compiler);
585 compiler.pow(data, 0, -1, result.data, 0);
586 return result;
587 }
588
589 /** {@inheritDoc}
590 * @since 3.2
591 */
592 public DerivativeStructure sqrt() {
593 return rootN(2);
594 }
595
596 /** {@inheritDoc}
597 * @since 3.2
598 */
599 public DerivativeStructure cbrt() {
600 return rootN(3);
601 }
602
603 /** {@inheritDoc}
604 * @since 3.2
605 */
606 public DerivativeStructure rootN(final int n) {
607 final DerivativeStructure result = new DerivativeStructure(compiler);
608 compiler.rootN(data, 0, n, result.data, 0);
609 return result;
610 }
611
612 /** {@inheritDoc} */
613 public Field<DerivativeStructure> getField() {
614 return new Field<DerivativeStructure>() {
615
616 /** {@inheritDoc} */
617 public DerivativeStructure getZero() {
618 return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 0.0);
619 }
620
621 /** {@inheritDoc} */
622 public DerivativeStructure getOne() {
623 return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 1.0);
624 }
625
626 /** {@inheritDoc} */
627 public Class<? extends FieldElement<DerivativeStructure>> getRuntimeClass() {
628 return DerivativeStructure.class;
629 }
630
631 };
632 }
633
634 /** {@inheritDoc}
635 * @since 3.2
636 */
637 public DerivativeStructure pow(final double p) {
638 final DerivativeStructure result = new DerivativeStructure(compiler);
639 compiler.pow(data, 0, p, result.data, 0);
640 return result;
641 }
642
643 /** {@inheritDoc}
644 * @since 3.2
645 */
646 public DerivativeStructure pow(final int n) {
647 final DerivativeStructure result = new DerivativeStructure(compiler);
648 compiler.pow(data, 0, n, result.data, 0);
649 return result;
650 }
651
652 /** {@inheritDoc}
653 * @exception DimensionMismatchException if number of free parameters
654 * or orders do not match
655 * @since 3.2
656 */
657 public DerivativeStructure pow(final DerivativeStructure e)
658 throws DimensionMismatchException {
659 compiler.checkCompatibility(e.compiler);
660 final DerivativeStructure result = new DerivativeStructure(compiler);
661 compiler.pow(data, 0, e.data, 0, result.data, 0);
662 return result;
663 }
664
665 /** {@inheritDoc}
666 * @since 3.2
667 */
668 public DerivativeStructure exp() {
669 final DerivativeStructure result = new DerivativeStructure(compiler);
670 compiler.exp(data, 0, result.data, 0);
671 return result;
672 }
673
674 /** {@inheritDoc}
675 * @since 3.2
676 */
677 public DerivativeStructure expm1() {
678 final DerivativeStructure result = new DerivativeStructure(compiler);
679 compiler.expm1(data, 0, result.data, 0);
680 return result;
681 }
682
683 /** {@inheritDoc}
684 * @since 3.2
685 */
686 public DerivativeStructure log() {
687 final DerivativeStructure result = new DerivativeStructure(compiler);
688 compiler.log(data, 0, result.data, 0);
689 return result;
690 }
691
692 /** {@inheritDoc}
693 * @since 3.2
694 */
695 public DerivativeStructure log1p() {
696 final DerivativeStructure result = new DerivativeStructure(compiler);
697 compiler.log1p(data, 0, result.data, 0);
698 return result;
699 }
700
701 /** Base 10 logarithm.
702 * @return base 10 logarithm of the instance
703 */
704 public DerivativeStructure log10() {
705 final DerivativeStructure result = new DerivativeStructure(compiler);
706 compiler.log10(data, 0, result.data, 0);
707 return result;
708 }
709
710 /** {@inheritDoc}
711 * @since 3.2
712 */
713 public DerivativeStructure cos() {
714 final DerivativeStructure result = new DerivativeStructure(compiler);
715 compiler.cos(data, 0, result.data, 0);
716 return result;
717 }
718
719 /** {@inheritDoc}
720 * @since 3.2
721 */
722 public DerivativeStructure sin() {
723 final DerivativeStructure result = new DerivativeStructure(compiler);
724 compiler.sin(data, 0, result.data, 0);
725 return result;
726 }
727
728 /** {@inheritDoc}
729 * @since 3.2
730 */
731 public DerivativeStructure tan() {
732 final DerivativeStructure result = new DerivativeStructure(compiler);
733 compiler.tan(data, 0, result.data, 0);
734 return result;
735 }
736
737 /** {@inheritDoc}
738 * @since 3.2
739 */
740 public DerivativeStructure acos() {
741 final DerivativeStructure result = new DerivativeStructure(compiler);
742 compiler.acos(data, 0, result.data, 0);
743 return result;
744 }
745
746 /** {@inheritDoc}
747 * @since 3.2
748 */
749 public DerivativeStructure asin() {
750 final DerivativeStructure result = new DerivativeStructure(compiler);
751 compiler.asin(data, 0, result.data, 0);
752 return result;
753 }
754
755 /** {@inheritDoc}
756 * @since 3.2
757 */
758 public DerivativeStructure atan() {
759 final DerivativeStructure result = new DerivativeStructure(compiler);
760 compiler.atan(data, 0, result.data, 0);
761 return result;
762 }
763
764 /** {@inheritDoc}
765 * @since 3.2
766 */
767 public DerivativeStructure atan2(final DerivativeStructure x)
768 throws DimensionMismatchException {
769 compiler.checkCompatibility(x.compiler);
770 final DerivativeStructure result = new DerivativeStructure(compiler);
771 compiler.atan2(data, 0, x.data, 0, result.data, 0);
772 return result;
773 }
774
775 /** Two arguments arc tangent operation.
776 * @param y first argument of the arc tangent
777 * @param x second argument of the arc tangent
778 * @return atan2(y, x)
779 * @exception DimensionMismatchException if number of free parameters
780 * or orders do not match
781 * @since 3.2
782 */
783 public static DerivativeStructure atan2(final DerivativeStructure y, final DerivativeStructure x)
784 throws DimensionMismatchException {
785 return y.atan2(x);
786 }
787
788 /** {@inheritDoc}
789 * @since 3.2
790 */
791 public DerivativeStructure cosh() {
792 final DerivativeStructure result = new DerivativeStructure(compiler);
793 compiler.cosh(data, 0, result.data, 0);
794 return result;
795 }
796
797 /** {@inheritDoc}
798 * @since 3.2
799 */
800 public DerivativeStructure sinh() {
801 final DerivativeStructure result = new DerivativeStructure(compiler);
802 compiler.sinh(data, 0, result.data, 0);
803 return result;
804 }
805
806 /** {@inheritDoc}
807 * @since 3.2
808 */
809 public DerivativeStructure tanh() {
810 final DerivativeStructure result = new DerivativeStructure(compiler);
811 compiler.tanh(data, 0, result.data, 0);
812 return result;
813 }
814
815 /** {@inheritDoc}
816 * @since 3.2
817 */
818 public DerivativeStructure acosh() {
819 final DerivativeStructure result = new DerivativeStructure(compiler);
820 compiler.acosh(data, 0, result.data, 0);
821 return result;
822 }
823
824 /** {@inheritDoc}
825 * @since 3.2
826 */
827 public DerivativeStructure asinh() {
828 final DerivativeStructure result = new DerivativeStructure(compiler);
829 compiler.asinh(data, 0, result.data, 0);
830 return result;
831 }
832
833 /** {@inheritDoc}
834 * @since 3.2
835 */
836 public DerivativeStructure atanh() {
837 final DerivativeStructure result = new DerivativeStructure(compiler);
838 compiler.atanh(data, 0, result.data, 0);
839 return result;
840 }
841
842 /** Convert radians to degrees, with error of less than 0.5 ULP
843 * @return instance converted into degrees
844 */
845 public DerivativeStructure toDegrees() {
846 final DerivativeStructure ds = new DerivativeStructure(compiler);
847 for (int i = 0; i < ds.data.length; ++i) {
848 ds.data[i] = FastMath.toDegrees(data[i]);
849 }
850 return ds;
851 }
852
853 /** Convert degrees to radians, with error of less than 0.5 ULP
854 * @return instance converted into radians
855 */
856 public DerivativeStructure toRadians() {
857 final DerivativeStructure ds = new DerivativeStructure(compiler);
858 for (int i = 0; i < ds.data.length; ++i) {
859 ds.data[i] = FastMath.toRadians(data[i]);
860 }
861 return ds;
862 }
863
864 /** Evaluate Taylor expansion a derivative structure.
865 * @param delta parameters offsets (Δx, Δy, ...)
866 * @return value of the Taylor expansion at x + Δx, y + Δy, ...
867 * @throws MathArithmeticException if factorials becomes too large
868 */
869 public double taylor(final double ... delta) throws MathArithmeticException {
870 return compiler.taylor(data, 0, delta);
871 }
872
873 /** {@inheritDoc}
874 * @exception DimensionMismatchException if number of free parameters
875 * or orders do not match
876 * @since 3.2
877 */
878 public DerivativeStructure linearCombination(final DerivativeStructure[] a, final DerivativeStructure[] b)
879 throws DimensionMismatchException {
880
881 // compute an accurate value, taking care of cancellations
882 final double[] aDouble = new double[a.length];
883 for (int i = 0; i < a.length; ++i) {
884 aDouble[i] = a[i].getValue();
885 }
886 final double[] bDouble = new double[b.length];
887 for (int i = 0; i < b.length; ++i) {
888 bDouble[i] = b[i].getValue();
889 }
890 final double accurateValue = MathArrays.linearCombination(aDouble, bDouble);
891
892 // compute a simple value, with all partial derivatives
893 DerivativeStructure simpleValue = a[0].getField().getZero();
894 for (int i = 0; i < a.length; ++i) {
895 simpleValue = simpleValue.add(a[i].multiply(b[i]));
896 }
897
898 // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
899 final double[] all = simpleValue.getAllDerivatives();
900 all[0] = accurateValue;
901 return new DerivativeStructure(simpleValue.getFreeParameters(), simpleValue.getOrder(), all);
902
903 }
904
905 /** {@inheritDoc}
906 * @exception DimensionMismatchException if number of free parameters
907 * or orders do not match
908 * @since 3.2
909 */
910 public DerivativeStructure linearCombination(final double[] a, final DerivativeStructure[] b)
911 throws DimensionMismatchException {
912
913 // compute an accurate value, taking care of cancellations
914 final double[] bDouble = new double[b.length];
915 for (int i = 0; i < b.length; ++i) {
916 bDouble[i] = b[i].getValue();
917 }
918 final double accurateValue = MathArrays.linearCombination(a, bDouble);
919
920 // compute a simple value, with all partial derivatives
921 DerivativeStructure simpleValue = b[0].getField().getZero();
922 for (int i = 0; i < a.length; ++i) {
923 simpleValue = simpleValue.add(b[i].multiply(a[i]));
924 }
925
926 // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
927 final double[] all = simpleValue.getAllDerivatives();
928 all[0] = accurateValue;
929 return new DerivativeStructure(simpleValue.getFreeParameters(), simpleValue.getOrder(), all);
930
931 }
932
933 /** {@inheritDoc}
934 * @exception DimensionMismatchException if number of free parameters
935 * or orders do not match
936 * @since 3.2
937 */
938 public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
939 final DerivativeStructure a2, final DerivativeStructure b2)
940 throws DimensionMismatchException {
941
942 // compute an accurate value, taking care of cancellations
943 final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
944 a2.getValue(), b2.getValue());
945
946 // compute a simple value, with all partial derivatives
947 final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2));
948
949 // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
950 final double[] all = simpleValue.getAllDerivatives();
951 all[0] = accurateValue;
952 return new DerivativeStructure(getFreeParameters(), getOrder(), all);
953
954 }
955
956 /** {@inheritDoc}
957 * @exception DimensionMismatchException if number of free parameters
958 * or orders do not match
959 * @since 3.2
960 */
961 public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
962 final double a2, final DerivativeStructure b2)
963 throws DimensionMismatchException {
964
965 // compute an accurate value, taking care of cancellations
966 final double accurateValue = MathArrays.linearCombination(a1, b1.getValue(),
967 a2, b2.getValue());
968
969 // compute a simple value, with all partial derivatives
970 final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2));
971
972 // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
973 final double[] all = simpleValue.getAllDerivatives();
974 all[0] = accurateValue;
975 return new DerivativeStructure(getFreeParameters(), getOrder(), all);
976
977 }
978
979 /** {@inheritDoc}
980 * @exception DimensionMismatchException if number of free parameters
981 * or orders do not match
982 * @since 3.2
983 */
984 public DerivativeStructure linearCombination(final DerivativeStructure a1, final DerivativeStructure b1,
985 final DerivativeStructure a2, final DerivativeStructure b2,
986 final DerivativeStructure a3, final DerivativeStructure b3)
987 throws DimensionMismatchException {
988
989 // compute an accurate value, taking care of cancellations
990 final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
991 a2.getValue(), b2.getValue(),
992 a3.getValue(), b3.getValue());
993
994 // compute a simple value, with all partial derivatives
995 final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3));
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 double a1, final DerivativeStructure b1,
1010 final double a2, final DerivativeStructure b2,
1011 final double a3, final DerivativeStructure b3)
1012 throws DimensionMismatchException {
1013
1014 // compute an accurate value, taking care of cancellations
1015 final double accurateValue = MathArrays.linearCombination(a1, b1.getValue(),
1016 a2, b2.getValue(),
1017 a3, b3.getValue());
1018
1019 // compute a simple value, with all partial derivatives
1020 final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3));
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 DerivativeStructure a1, final DerivativeStructure b1,
1035 final DerivativeStructure a2, final DerivativeStructure b2,
1036 final DerivativeStructure a3, final DerivativeStructure b3,
1037 final DerivativeStructure a4, final DerivativeStructure b4)
1038 throws DimensionMismatchException {
1039
1040 // compute an accurate value, taking care of cancellations
1041 final double accurateValue = MathArrays.linearCombination(a1.getValue(), b1.getValue(),
1042 a2.getValue(), b2.getValue(),
1043 a3.getValue(), b3.getValue(),
1044 a4.getValue(), b4.getValue());
1045
1046 // compute a simple value, with all partial derivatives
1047 final DerivativeStructure simpleValue = a1.multiply(b1).add(a2.multiply(b2)).add(a3.multiply(b3)).add(a4.multiply(b4));
1048
1049 // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1050 final double[] all = simpleValue.getAllDerivatives();
1051 all[0] = accurateValue;
1052 return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1053
1054 }
1055
1056 /** {@inheritDoc}
1057 * @exception DimensionMismatchException if number of free parameters
1058 * or orders do not match
1059 * @since 3.2
1060 */
1061 public DerivativeStructure linearCombination(final double a1, final DerivativeStructure b1,
1062 final double a2, final DerivativeStructure b2,
1063 final double a3, final DerivativeStructure b3,
1064 final double a4, final DerivativeStructure b4)
1065 throws DimensionMismatchException {
1066
1067 // compute an accurate value, taking care of cancellations
1068 final double accurateValue = MathArrays.linearCombination(a1, b1.getValue(),
1069 a2, b2.getValue(),
1070 a3, b3.getValue(),
1071 a4, b4.getValue());
1072
1073 // compute a simple value, with all partial derivatives
1074 final DerivativeStructure simpleValue = b1.multiply(a1).add(b2.multiply(a2)).add(b3.multiply(a3)).add(b4.multiply(a4));
1075
1076 // create a result with accurate value and all derivatives (not necessarily as accurate as the value)
1077 final double[] all = simpleValue.getAllDerivatives();
1078 all[0] = accurateValue;
1079 return new DerivativeStructure(getFreeParameters(), getOrder(), all);
1080
1081 }
1082
1083 /**
1084 * Test for the equality of two derivative structures.
1085 * <p>
1086 * Derivative structures are considered equal if they have the same number
1087 * of free parameters, the same derivation order, and the same derivatives.
1088 * </p>
1089 * @param other Object to test for equality to this
1090 * @return true if two derivative structures are equal
1091 * @since 3.2
1092 */
1093 @Override
1094 public boolean equals(Object other) {
1095
1096 if (this == other) {
1097 return true;
1098 }
1099
1100 if (other instanceof DerivativeStructure) {
1101 final DerivativeStructure rhs = (DerivativeStructure)other;
1102 return (getFreeParameters() == rhs.getFreeParameters()) &&
1103 (getOrder() == rhs.getOrder()) &&
1104 MathArrays.equals(data, rhs.data);
1105 }
1106
1107 return false;
1108
1109 }
1110
1111 /**
1112 * Get a hashCode for the derivative structure.
1113 * @return a hash code value for this object
1114 * @since 3.2
1115 */
1116 @Override
1117 public int hashCode() {
1118 return 227 + 229 * getFreeParameters() + 233 * getOrder() + 239 * MathUtils.hash(data);
1119 }
1120
1121 /**
1122 * Replace the instance with a data transfer object for serialization.
1123 * @return data transfer object that will be serialized
1124 */
1125 private Object writeReplace() {
1126 return new DataTransferObject(compiler.getFreeParameters(), compiler.getOrder(), data);
1127 }
1128
1129 /** Internal class used only for serialization. */
1130 private static class DataTransferObject implements Serializable {
1131
1132 /** Serializable UID. */
1133 private static final long serialVersionUID = 20120730L;
1134
1135 /** Number of variables.
1136 * @serial
1137 */
1138 private final int variables;
1139
1140 /** Derivation order.
1141 * @serial
1142 */
1143 private final int order;
1144
1145 /** Partial derivatives.
1146 * @serial
1147 */
1148 private final double[] data;
1149
1150 /** Simple constructor.
1151 * @param variables number of variables
1152 * @param order derivation order
1153 * @param data partial derivatives
1154 */
1155 public DataTransferObject(final int variables, final int order, final double[] data) {
1156 this.variables = variables;
1157 this.order = order;
1158 this.data = data;
1159 }
1160
1161 /** Replace the deserialized data transfer object with a {@link DerivativeStructure}.
1162 * @return replacement {@link DerivativeStructure}
1163 */
1164 private Object readResolve() {
1165 return new DerivativeStructure(variables, order, data);
1166 }
1167
1168 }
1169
1170 }