1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math3.analysis.differentiation; 18 19 import java.io.Serializable; 20 21 import org.apache.commons.math3.Field; 22 import org.apache.commons.math3.FieldElement; 23 import org.apache.commons.math3.RealFieldElement; 24 import org.apache.commons.math3.exception.DimensionMismatchException; 25 import org.apache.commons.math3.exception.MathArithmeticException; 26 import org.apache.commons.math3.exception.NumberIsTooLargeException; 27 import org.apache.commons.math3.util.FastMath; 28 import org.apache.commons.math3.util.MathArrays; 29 import org.apache.commons.math3.util.MathUtils; 30 31 /** Class representing both the value and the differentials of a function. 32 * <p>This class is the workhorse of the differentiation package.</p> 33 * <p>This class is an implementation of the extension to Rall's 34 * numbers described in Dan Kalman's paper <a 35 * href="http://www1.american.edu/cas/mathstat/People/kalman/pdffiles/mmgautodiff.pdf">Doubly 36 * Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75, 37 * no. 3, June 2002. Rall's numbers are an extension to the real numbers used 38 * throughout mathematical expressions; they hold the derivative together with the 39 * value of a function. Dan Kalman's derivative structures hold all partial derivatives 40 * up to any specified order, with respect to any number of free parameters. Rall's 41 * numbers therefore can be seen as derivative structures for order one derivative and 42 * one free parameter, and real numbers can be seen as derivative structures with zero 43 * order derivative and no free parameters.</p> 44 * <p>{@link DerivativeStructure} instances can be used directly thanks to 45 * the arithmetic operators to the mathematical functions provided as 46 * methods by this class (+, -, *, /, %, sin, cos ...).</p> 47 * <p>Implementing complex expressions by hand using these classes is 48 * a tedious and error-prone task but has the advantage of having no limitation 49 * on the derivation order despite no requiring users to compute the derivatives by 50 * themselves. Implementing complex expression can also be done by developing computation 51 * code using standard primitive double values and to use {@link 52 * UnivariateFunctionDifferentiator differentiators} to create the {@link 53 * DerivativeStructure}-based instances. This method is simpler but may be limited in 54 * the accuracy and derivation orders and may be computationally intensive (this is 55 * typically the case for {@link FiniteDifferencesDifferentiator finite differences 56 * differentiator}.</p> 57 * <p>Instances of this class are guaranteed to be immutable.</p> 58 * @see DSCompiler 59 * @since 3.1 60 */ 61 public class DerivativeStructure implements RealFieldElement<DerivativeStructure>, Serializable { 62 63 /** Serializable UID. */ 64 private static final long serialVersionUID = 20120730L; 65 66 /** Compiler for the current dimensions. */ 67 private transient DSCompiler compiler; 68 69 /** Combined array holding all values. */ 70 private final double[] data; 71 72 /** Build an instance with all values and derivatives set to 0. 73 * @param compiler compiler to use for computation 74 */ DerivativeStructure(final DSCompiler compiler)75 private DerivativeStructure(final DSCompiler compiler) { 76 this.compiler = compiler; 77 this.data = new double[compiler.getSize()]; 78 } 79 80 /** Build an instance with all values and derivatives set to 0. 81 * @param parameters number of free parameters 82 * @param order derivation order 83 * @throws NumberIsTooLargeException if order is too large 84 */ DerivativeStructure(final int parameters, final int order)85 public DerivativeStructure(final int parameters, final int order) 86 throws NumberIsTooLargeException { 87 this(DSCompiler.getCompiler(parameters, order)); 88 } 89 90 /** Build an instance representing a constant value. 91 * @param parameters number of free parameters 92 * @param order derivation order 93 * @param value value of the constant 94 * @throws NumberIsTooLargeException if order is too large 95 * @see #DerivativeStructure(int, int, int, double) 96 */ DerivativeStructure(final int parameters, final int order, final double value)97 public DerivativeStructure(final int parameters, final int order, final double value) 98 throws NumberIsTooLargeException { 99 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 */ DerivativeStructure(final int parameters, final int order, final int index, final double value)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 */ DerivativeStructure(final double a1, final DerivativeStructure ds1, final double a2, final DerivativeStructure ds2)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 */ DerivativeStructure(final double a1, final DerivativeStructure ds1, final double a2, final DerivativeStructure ds2, final double a3, final DerivativeStructure ds3)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 */ DerivativeStructure(final double a1, final DerivativeStructure ds1, final double a2, final DerivativeStructure ds2, final double a3, final DerivativeStructure ds3, final double a4, final DerivativeStructure ds4)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 */ DerivativeStructure(final int parameters, final int order, final double ... derivatives)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 */ DerivativeStructure(final DerivativeStructure ds)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 */ getFreeParameters()223 public int getFreeParameters() { 224 return compiler.getFreeParameters(); 225 } 226 227 /** Get the derivation order. 228 * @return derivation order 229 */ getOrder()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 */ createConstant(final double c)244 public DerivativeStructure createConstant(final double c) { 245 return new DerivativeStructure(getFreeParameters(), getOrder(), c); 246 } 247 248 /** {@inheritDoc} 249 * @since 3.2 250 */ getReal()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 */ getValue()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 */ getPartialDerivative(final int ... orders)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 */ getAllDerivatives()282 public double[] getAllDerivatives() { 283 return data.clone(); 284 } 285 286 /** {@inheritDoc} 287 * @since 3.2 288 */ add(final double a)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 */ add(final DerivativeStructure a)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 */ subtract(final double a)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 */ subtract(final DerivativeStructure a)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} */ multiply(final int n)327 public DerivativeStructure multiply(final int n) { 328 return multiply((double) n); 329 } 330 331 /** {@inheritDoc} 332 * @since 3.2 333 */ multiply(final double a)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 */ multiply(final DerivativeStructure a)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 */ divide(final double a)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 */ divide(final DerivativeStructure a)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} */ remainder(final double a)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 */ remainder(final DerivativeStructure a)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} */ negate()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 */ abs()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 */ ceil()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 */ floor()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 */ rint()439 public DerivativeStructure rint() { 440 return new DerivativeStructure(compiler.getFreeParameters(), 441 compiler.getOrder(), 442 FastMath.rint(data[0])); 443 } 444 445 /** {@inheritDoc} */ round()446 public long round() { 447 return FastMath.round(data[0]); 448 } 449 450 /** {@inheritDoc} 451 * @since 3.2 452 */ signum()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 */ copySign(final DerivativeStructure sign)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 */ copySign(final double sign)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 */ getExponent()491 public int getExponent() { 492 return FastMath.getExponent(data[0]); 493 } 494 495 /** {@inheritDoc} 496 * @since 3.2 497 */ scalb(final int n)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 */ hypot(final DerivativeStructure y)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> +<i>y</i><sup>2</sup>) 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> +<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 */ hypot(final DerivativeStructure x, final DerivativeStructure y)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 */ compose(final double ... f)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} */ reciprocal()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 */ sqrt()605 public DerivativeStructure sqrt() { 606 return rootN(2); 607 } 608 609 /** {@inheritDoc} 610 * @since 3.2 611 */ cbrt()612 public DerivativeStructure cbrt() { 613 return rootN(3); 614 } 615 616 /** {@inheritDoc} 617 * @since 3.2 618 */ rootN(final int n)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} */ getField()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 (Δx, Δy, ...) 891 * @return value of the Taylor expansion at x + Δx, y + Δ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 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 } 1196