You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2015/05/17 16:35:40 UTC
[1/3] [math] Added a fast implementation of IEEEremainder in FastMath.
Repository: commons-math
Updated Branches:
refs/heads/master 83c61da2c -> 0c0455fd6
Added a fast implementation of IEEEremainder in FastMath.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/15a24dc0
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/15a24dc0
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/15a24dc0
Branch: refs/heads/master
Commit: 15a24dc0fc9067fabe00857f8546aee4a5234a78
Parents: 83c61da
Author: Luc Maisonobe <lu...@apache.org>
Authored: Sat May 16 14:22:46 2015 +0200
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Sat May 16 14:24:13 2015 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 3 ++
.../org/apache/commons/math4/util/FastMath.java | 17 +++++++--
.../userguide/FastMathTestPerformance.java | 40 ++++++++++++++++++--
3 files changed, 54 insertions(+), 6 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/15a24dc0/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 6afa831..1b448a8 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -54,6 +54,9 @@ If the output is not quite correct, check for invisible trailing spaces!
</release>
<release version="4.0" date="XXXX-XX-XX" description="">
+ <action dev="luc" type="add" >
+ Added a fast implementation of IEEEremainder in FastMath.
+ </action>
<action dev="luc" type="fix" issue="MATH-1222" due-to="Benedikt Ritter">
Use Double.isNaN rather than x != x in FastMath.
</action>
http://git-wip-us.apache.org/repos/asf/commons-math/blob/15a24dc0/src/main/java/org/apache/commons/math4/util/FastMath.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/util/FastMath.java b/src/main/java/org/apache/commons/math4/util/FastMath.java
index cd56192..3e96f69 100644
--- a/src/main/java/org/apache/commons/math4/util/FastMath.java
+++ b/src/main/java/org/apache/commons/math4/util/FastMath.java
@@ -3628,13 +3628,24 @@ public class FastMath {
* <li>If the dividend is finite and the divisor is an infinity, the result equals the dividend.</li>
* <li>If the dividend is a zero and the divisor is finite, the result equals the dividend.</li>
* </ul>
- * <p><b>Note:</b> this implementation currently delegates to {@link StrictMath#IEEEremainder}
* @param dividend the number to be divided
* @param divisor the number by which to divide
* @return the remainder, rounded
*/
- public static double IEEEremainder(double dividend, double divisor) {
- return StrictMath.IEEEremainder(dividend, divisor); // TODO provide our own implementation
+ public static double IEEEremainder(final double dividend, final double divisor) {
+ if (getExponent(dividend) == 1024 || getExponent(divisor) == 1024 || divisor == 0.0) {
+ // we are in one of the special cases
+ if (Double.isInfinite(divisor) && !Double.isInfinite(dividend)) {
+ return dividend;
+ } else {
+ return Double.NaN;
+ }
+ } else {
+ // we are in the general case
+ final double n = FastMath.rint(dividend / divisor);
+ final double remainder = Double.isInfinite(n) ? 0.0 : dividend - divisor * n;
+ return (remainder == 0) ? FastMath.copySign(remainder, dividend) : remainder;
+ }
}
/** Convert a long to interger, detecting overflows
http://git-wip-us.apache.org/repos/asf/commons-math/blob/15a24dc0/src/userguide/java/org/apache/commons/math4/userguide/FastMathTestPerformance.java
----------------------------------------------------------------------
diff --git a/src/userguide/java/org/apache/commons/math4/userguide/FastMathTestPerformance.java b/src/userguide/java/org/apache/commons/math4/userguide/FastMathTestPerformance.java
index 8af6fca..697bd60 100644
--- a/src/userguide/java/org/apache/commons/math4/userguide/FastMathTestPerformance.java
+++ b/src/userguide/java/org/apache/commons/math4/userguide/FastMathTestPerformance.java
@@ -28,9 +28,9 @@ public class FastMathTestPerformance {
private static final double F1 = 1d / RUNS;
// Header format
- private static final String FMT_HDR = "%-5s %13s %13s %13s Runs=%d Java %s (%s) %s (%s)";
+ private static final String FMT_HDR = "%-13s %13s %13s %13s Runs=%d Java %s (%s) %s (%s)";
// Detail format
- private static final String FMT_DTL = "%-5s %6d %6.1f %6d %6.4f %6d %6.4f";
+ private static final String FMT_DTL = "%-13s %6d %6.1f %6d %6.4f %6d %6.4f";
public static void main(String[] args) {
System.out.println(String.format(FMT_HDR,
@@ -60,6 +60,7 @@ public class FastMathTestPerformance {
testSqrt();
testTan();
testTanh();
+ testIEEEremainder();
}
@@ -429,7 +430,40 @@ public class FastMathTestPerformance {
report("hypot",strictTime,fastTime,mathTime);
assertTrue(!Double.isNaN(x));
}
-
+
+ private static void testIEEEremainder() {
+ double x = 0;
+ long time = System.nanoTime();
+ int max = (int) FastMath.floor(FastMath.sqrt(RUNS));
+ for (int i = 0; i < max; i++) {
+ for (int j = 0; j < max; j++) {
+ x += StrictMath.IEEEremainder((i - max/2) * (100.0 / max), (j + 1) * (100.0 / max));
+ }
+ }
+ long strictTime = System.nanoTime() - time;
+
+ x = 0;
+ time = System.nanoTime();
+ for (int i = 0; i < max; i++) {
+ for (int j = 0; j < max; j++) {
+ x += FastMath.IEEEremainder((i - max/2) * (100.0 / max), (j + 1) * (100.0 / max));
+ }
+ }
+ long fastTime = System.nanoTime() - time;
+
+ x = 0;
+ time = System.nanoTime();
+ for (int i = 0; i < max; i++) {
+ for (int j = 0; j < max; j++) {
+ x += Math.IEEEremainder((i - max/2) * (100.0 / max), (j + 1) * (100.0 / max));
+ }
+ }
+ long mathTime = System.nanoTime() - time;
+
+ report("IEEEremainder",strictTime,fastTime,mathTime);
+ assertTrue(!Double.isNaN(x));
+ }
+
private static void testCbrt() {
double x = 0;
long time = System.nanoTime();
[3/3] [math] Attempt to re-implement the pow function.
Posted by lu...@apache.org.
Attempt to re-implement the pow function.
The attempts are motivated by platform-specific failures, which seem to
be platform-specific, and probably due to JIT optimization bugs.
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/0c0455fd
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/0c0455fd
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/0c0455fd
Branch: refs/heads/master
Commit: 0c0455fd66eed4113bf53829cfdb501901cc70f9
Parents: 9b6a649
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue May 12 15:08:07 2015 +0200
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Sun May 17 16:31:12 2015 +0200
----------------------------------------------------------------------
src/changes/changes.xml | 7 +
.../org/apache/commons/math4/util/FastMath.java | 422 +++++++++++--------
.../distribution/GammaDistributionTest.java | 2 +-
.../scalar/noderiv/BOBYQAOptimizerTest.java | 2 +-
.../apache/commons/math4/util/FastMathTest.java | 54 ++-
5 files changed, 302 insertions(+), 185 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0c0455fd/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 1b448a8..66f0175 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -55,6 +55,13 @@ If the output is not quite correct, check for invisible trailing spaces!
<release version="4.0" date="XXXX-XX-XX" description="">
<action dev="luc" type="add" >
+ Reimplemented pow(double, double) in FastMath, for better accuracy in
+ integral power cases and trying to fix erroneous JIT optimization again.
+ </action>
+ <action dev="luc" type="add" >
+ Added a pow(double, long) method in FastMath.
+ </action>
+ <action dev="luc" type="add" >
Added a fast implementation of IEEEremainder in FastMath.
</action>
<action dev="luc" type="fix" issue="MATH-1222" due-to="Benedikt Ritter">
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0c0455fd/src/main/java/org/apache/commons/math4/util/FastMath.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/util/FastMath.java b/src/main/java/org/apache/commons/math4/util/FastMath.java
index 3e96f69..a4a9a1b 100644
--- a/src/main/java/org/apache/commons/math4/util/FastMath.java
+++ b/src/main/java/org/apache/commons/math4/util/FastMath.java
@@ -315,10 +315,17 @@ public class FastMath {
/** Mask used to clear the non-sign part of a long. */
private static final long MASK_NON_SIGN_LONG = 0x7fffffffffffffffl;
+ /** Mask used to extract exponent from double bits. */
+ private static final long MASK_DOUBLE_EXPONENT = 0x7ff0000000000000L;
+
+ /** Mask used to extract mantissa from double bits. */
+ private static final long MASK_DOUBLE_MANTISSA = 0x000fffffffffffffL;
+
+ /** Mask used to add implicit high order bit for normalized double. */
+ private static final long IMPLICIT_HIGH_BIT = 0x0010000000000000L;
+
/** 2^52 - double numbers this large must be integral (no fraction) or NaN or Infinite */
private static final double TWO_POWER_52 = 4503599627370496.0;
- /** 2^53 - double numbers this large must be even. */
- private static final double TWO_POWER_53 = 2 * TWO_POWER_52;
/** Constant: {@value}. */
private static final double F_1_3 = 1d / 3d;
@@ -1458,144 +1465,141 @@ public class FastMath {
* @return double
*/
public static double pow(final double x, final double y) {
- final double lns[] = new double[2];
- if (y == 0.0) {
+ if (y == 0) {
+ // y = -0 or y = +0
return 1.0;
- } else if (Double.isNaN(x)) {
- return x;
- } else if (Double.isNaN(y)) {
- return y;
- } else if (x == 0) {
- long bits = Double.doubleToRawLongBits(x);
- if ((bits & 0x8000000000000000L) != 0) {
- // -zero
- long yi = (long) y;
-
- if (y < 0 && y == yi && (yi & 1) == 1) {
- return Double.NEGATIVE_INFINITY;
- }
+ } else {
- if (y > 0 && y == yi && (yi & 1) == 1) {
- return -0.0;
+ final long yBits = Double.doubleToRawLongBits(y);
+ final int yRawExp = (int) ((yBits & MASK_DOUBLE_EXPONENT) >> 52);
+ final long yRawMantissa = yBits & MASK_DOUBLE_MANTISSA;
+ final long xBits = Double.doubleToRawLongBits(x);
+ final int xRawExp = (int) ((xBits & MASK_DOUBLE_EXPONENT) >> 52);
+ final long xRawMantissa = xBits & MASK_DOUBLE_MANTISSA;
+
+ if (yRawExp > 1085) {
+ // y is either a very large integral value that does not fit in a long or it is a special number
+
+ if ((yRawExp == 2047 && yRawMantissa != 0) ||
+ (xRawExp == 2047 && xRawMantissa != 0)) {
+ // NaN
+ return Double.NaN;
+ } else if (xRawExp == 1023 && xRawMantissa == 0) {
+ // x = -1.0 or x = +1.0
+ if (yRawExp == 2047) {
+ // y is infinite
+ return Double.NaN;
+ } else {
+ // y is a large even integer
+ return 1.0;
+ }
+ } else {
+ // the absolute value of x is either greater or smaller than 1.0
+
+ // if yRawExp == 2047 and mantissa is 0, y = -infinity or y = +infinity
+ // if 1085 < yRawExp < 2047, y is simply a large number, however, due to limited
+ // accuracy, at this magnitude it behaves just like infinity with regards to x
+ if ((y > 0) ^ (xRawExp < 1023)) {
+ // either y = +infinity (or large engouh) and abs(x) > 1.0
+ // or y = -infinity (or large engouh) and abs(x) < 1.0
+ return Double.POSITIVE_INFINITY;
+ } else {
+ // either y = +infinity (or large engouh) and abs(x) < 1.0
+ // or y = -infinity (or large engouh) and abs(x) > 1.0
+ return +0.0;
+ }
}
- }
- if (y < 0) {
- return Double.POSITIVE_INFINITY;
- }
- if (y > 0) {
- return 0.0;
- }
-
- return Double.NaN;
- } else if (x == Double.POSITIVE_INFINITY) {
- if (y < 0.0) {
- return 0.0;
- } else {
- return Double.POSITIVE_INFINITY;
- }
- } else if (y == Double.POSITIVE_INFINITY) {
- if (x * x == 1.0) {
- return Double.NaN;
- }
-
- if (x * x > 1.0) {
- return Double.POSITIVE_INFINITY;
} else {
- return 0.0;
- }
- } else if (x == Double.NEGATIVE_INFINITY) {
- if (y < 0) {
- long yi = (long) y;
- if (y == yi && (yi & 1) == 1) {
- return -0.0;
+ // y is a regular non-zero number
+
+ if (yRawExp >= 1023) {
+ // y may be an integral value, which should be handled specifically
+ final long yFullMantissa = IMPLICIT_HIGH_BIT | yRawMantissa;
+ if (yRawExp < 1075) {
+ // normal number with negative shift that may have a fractional part
+ final long integralMask = (-1L) << (1075 - yRawExp);
+ if ((yFullMantissa & integralMask) == yFullMantissa) {
+ // all fractional bits are 0, the number is really integral
+ final long l = yFullMantissa >> (1075 - yRawExp);
+ return FastMath.pow(x, (y < 0) ? -l : l);
+ }
+ } else {
+ // normal number with positive shift, always an integral value
+ // we know it fits in a primitive long because yRawExp > 1085 has been handled above
+ final long l = yFullMantissa << (yRawExp - 1075);
+ return FastMath.pow(x, (y < 0) ? -l : l);
+ }
}
- return 0.0;
- }
+ // y is a non-integral value
+
+ if (x == 0) {
+ // x = -0 or x = +0
+ // the integer powers have already been handled above
+ return y < 0 ? Double.POSITIVE_INFINITY : +0.0;
+ } else if (xRawExp == 2047) {
+ if (xRawMantissa == 0) {
+ // x = -infinity or x = +infinity
+ return (y < 0) ? +0.0 : Double.POSITIVE_INFINITY;
+ } else {
+ // NaN
+ return Double.NaN;
+ }
+ } else if (x < 0) {
+ // the integer powers have already been handled above
+ return Double.NaN;
+ } else {
- if (y > 0) {
- long yi = (long) y;
- if (y == yi && (yi & 1) == 1) {
- return Double.NEGATIVE_INFINITY;
- }
+ // this is the general case, for regular fractional numbers x and y
- return Double.POSITIVE_INFINITY;
- }
- } else if (y == Double.NEGATIVE_INFINITY) {
- if (x * x == 1.0) {
- return Double.NaN;
- }
+ // Split y into ya and yb such that y = ya+yb
+ final double tmp = y * HEX_40000000;
+ final double ya = (y + tmp) - tmp;
+ final double yb = y - ya;
- if (x * x < 1.0) {
- return Double.POSITIVE_INFINITY;
- } else {
- return 0.0;
- }
- } else if (x < 0) { // Handle special case x<0
- // y is an even integer in this case
- if (y >= TWO_POWER_53 || y <= -TWO_POWER_53) {
- return pow(-x, y);
- }
+ /* Compute ln(x) */
+ final double lns[] = new double[2];
+ final double lores = log(x, lns);
+ if (Double.isInfinite(lores)) { // don't allow this to be converted to NaN
+ return lores;
+ }
- if (y == (long) y) {
- // If y is an integer
- return ((long)y & 1) == 0 ? pow(-x, y) : -pow(-x, y);
- } else {
- return Double.NaN;
- }
- }
+ double lna = lns[0];
+ double lnb = lns[1];
- /* Split y into ya and yb such that y = ya+yb */
- double ya;
- double yb;
- if (y < 8e298 && y > -8e298) {
- double tmp1 = y * HEX_40000000;
- ya = y + tmp1 - tmp1;
- yb = y - ya;
- } else {
- double tmp1 = y * 9.31322574615478515625E-10;
- double tmp2 = tmp1 * 9.31322574615478515625E-10;
- ya = (tmp1 + tmp2 - tmp1) * HEX_40000000 * HEX_40000000;
- yb = y - ya;
- }
+ /* resplit lns */
+ final double tmp1 = lna * HEX_40000000;
+ final double tmp2 = (lna + tmp1) - tmp1;
+ lnb += lna - tmp2;
+ lna = tmp2;
- /* Compute ln(x) */
- final double lores = log(x, lns);
- if (Double.isInfinite(lores)){ // don't allow this to be converted to NaN
- return lores;
- }
+ // y*ln(x) = (aa+ab)
+ final double aa = lna * ya;
+ final double ab = lna * yb + lnb * ya + lnb * yb;
- double lna = lns[0];
- double lnb = lns[1];
+ lna = aa+ab;
+ lnb = -(lna - aa - ab);
- /* resplit lns */
- double tmp1 = lna * HEX_40000000;
- double tmp2 = lna + tmp1 - tmp1;
- lnb += lna - tmp2;
- lna = tmp2;
+ double z = 1.0 / 120.0;
+ z = z * lnb + (1.0 / 24.0);
+ z = z * lnb + (1.0 / 6.0);
+ z = z * lnb + 0.5;
+ z = z * lnb + 1.0;
+ z *= lnb;
- // y*ln(x) = (aa+ab)
- final double aa = lna * ya;
- final double ab = lna * yb + lnb * ya + lnb * yb;
+ final double result = exp(lna, z, null);
+ //result = result + result * z;
+ return result;
- lna = aa+ab;
- lnb = -(lna - aa - ab);
+ }
+ }
- double z = 1.0 / 120.0;
- z = z * lnb + (1.0 / 24.0);
- z = z * lnb + (1.0 / 6.0);
- z = z * lnb + 0.5;
- z = z * lnb + 1.0;
- z *= lnb;
+ }
- final double result = exp(lna, z, null);
- //result = result + result * z;
- return result;
}
-
/**
* Raise a double to an int power.
*
@@ -1605,68 +1609,150 @@ public class FastMath {
* @since 3.1
*/
public static double pow(double d, int e) {
+ return pow(d, (long) e);
+ }
+
+ /**
+ * Raise a double to a long power.
+ *
+ * @param d Number to raise.
+ * @param e Exponent.
+ * @return d<sup>e</sup>
+ * @since 4.0
+ */
+ public static double pow(double d, long e) {
if (e == 0) {
return 1.0;
- } else if (e < 0) {
- e = -e;
- d = 1.0 / d;
- }
-
- // split d as one 26 bits number and one 27 bits number
- // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
- final double d1High = Double.longBitsToDouble(Double.doubleToRawLongBits(d) & ((-1L) << 27));
- final double d1Low = d - d1High;
-
- // prepare result
- double resultHigh = 1;
- double resultLow = 0;
-
- // d^(2p)
- double d2p = d;
- double d2pHigh = d1High;
- double d2pLow = d1Low;
-
- while (e != 0) {
-
- if ((e & 0x1) != 0) {
- // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
- // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
- final double tmpHigh = resultHigh * d2p;
- final double rHH = Double.longBitsToDouble(Double.doubleToRawLongBits(resultHigh) & ((-1L) << 27));
- final double rHL = resultHigh - rHH;
- final double tmpLow = rHL * d2pLow - (((tmpHigh - rHH * d2pHigh) - rHL * d2pHigh) - rHH * d2pLow);
- resultHigh = tmpHigh;
- resultLow = resultLow * d2p + tmpLow;
- }
+ } else if (e > 0) {
+ return new Split(d).pow(e).full;
+ } else {
+ return new Split(d).reciprocal().pow(-e).full;
+ }
+ }
+
+ /** Class operator on double numbers split into one 26 bits number and one 27 bits number. */
+ private static class Split {
+
+ /** Split version of NaN. */
+ public static final Split NAN = new Split(Double.NaN, 0);
+
+ /** Split version of positive infinity. */
+ public static final Split POSITIVE_INFINITY = new Split(Double.POSITIVE_INFINITY, 0);
+
+ /** Split version of negative infinity. */
+ public static final Split NEGATIVE_INFINITY = new Split(Double.NEGATIVE_INFINITY, 0);
+
+ /** Full number. */
+ private final double full;
- // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
+ /** High order bits. */
+ private final double high;
+
+ /** Low order bits. */
+ private final double low;
+
+ /** Simple constructor.
+ * @param x number to split
+ */
+ public Split(final double x) {
+ full = x;
+ high = Double.longBitsToDouble(Double.doubleToRawLongBits(x) & ((-1L) << 27));
+ low = x - high;
+ }
+
+ /** Simple constructor.
+ * @param high high order bits
+ * @param low low order bits
+ */
+ public Split(final double high, final double low) {
+ this(high + low, high, low);
+ }
+
+ /** Simple constructor.
+ * @param full full number
+ * @param high high order bits
+ * @param low low order bits
+ */
+ public Split(final double full, final double high, final double low) {
+ this.full = full;
+ this.high = high;
+ this.low = low;
+ }
+
+ /** Multiply the instance by another one.
+ * @param b other instance to multiply by
+ * @return product
+ */
+ public Split multiply(final Split b) {
+ // beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
+ final Split mulBasic = new Split(full * b.full);
+ final double mulError = low * b.low - (((mulBasic.full - high * b.high) - low * b.high) - high * b.low);
+ return new Split(mulBasic.high, mulBasic.low + mulError);
+ }
+
+ /** Compute the reciprocal of the instance.
+ * @return reciprocal of the instance
+ */
+ public Split reciprocal() {
+
+ final double approximateInv = 1.0 / full;
+ final Split splitInv = new Split(approximateInv);
+
+ // if 1.0/d were computed perfectly, remultiplying it by d should give 1.0
+ // we want to estimate the error so we can fix the low order bits of approximateInvLow
// beware the following expressions must NOT be simplified, they rely on floating point arithmetic properties
- final double tmpHigh = d2pHigh * d2p;
- final double cD2pH = Double.longBitsToDouble(Double.doubleToRawLongBits(d2pHigh) & ((-1L) << 27));
- final double d2pHH = cD2pH - (cD2pH - d2pHigh);
- final double d2pHL = d2pHigh - d2pHH;
- final double tmpLow = d2pHL * d2pLow - (((tmpHigh - d2pHH * d2pHigh) - d2pHL * d2pHigh) - d2pHH * d2pLow);
- d2pHigh = Double.longBitsToDouble(Double.doubleToRawLongBits(tmpHigh) & ((-1L) << 27));
- d2pLow = d2pLow * d2p + tmpLow + (tmpHigh - d2pHigh);
- d2p = d2pHigh + d2pLow;
+ final Split product = multiply(splitInv);
+ final double error = (product.high - 1) + product.low;
- e >>= 1;
+ // better accuracy estimate of reciprocal
+ return Double.isNaN(error) ? splitInv : new Split(splitInv.high, splitInv.low - error / full);
}
- final double result = resultHigh + resultLow;
+ /** Computes this^e.
+ * @param e exponent (beware, here it MUST be > 0)
+ * @return d^e, split in high and low bits
+ * @since 4.0
+ */
+ private Split pow(final long e) {
- if (Double.isNaN(result)) {
- if (Double.isNaN(d)) {
- return Double.NaN;
+ // prepare result
+ Split result = new Split(1);
+
+ // d^(2p)
+ Split d2p = new Split(full, high, low);
+
+ for (long p = e; p != 0; p >>= 1) {
+
+ if ((p & 0x1) != 0) {
+ // accurate multiplication result = result * d^(2p) using Veltkamp TwoProduct algorithm
+ result = result.multiply(d2p);
+ }
+
+ // accurate squaring d^(2(p+1)) = d^(2p) * d^(2p) using Veltkamp TwoProduct algorithm
+ d2p = d2p.multiply(d2p);
+
+ }
+
+ if (Double.isNaN(result.full)) {
+ if (Double.isNaN(full)) {
+ return Split.NAN;
+ } else {
+ // some intermediate numbers exceeded capacity,
+ // and the low order bits became NaN (because infinity - infinity = NaN)
+ if (FastMath.abs(full) < 1) {
+ return new Split(FastMath.copySign(0.0, full), 0.0);
+ } else if (full < 0 && (e & 0x1) == 1) {
+ return Split.NEGATIVE_INFINITY;
+ } else {
+ return Split.POSITIVE_INFINITY;
+ }
+ }
} else {
- // some intermediate numbers exceeded capacity,
- // and the low order bits became NaN (because infinity - infinity = NaN)
- return (d < 0 && (e & 0x1) == 1) ? Double.NEGATIVE_INFINITY : Double.POSITIVE_INFINITY;
+ return result;
}
- } else {
- return result;
+
}
}
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0c0455fd/src/test/java/org/apache/commons/math4/distribution/GammaDistributionTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/distribution/GammaDistributionTest.java b/src/test/java/org/apache/commons/math4/distribution/GammaDistributionTest.java
index 6945e3f..437ed5a 100644
--- a/src/test/java/org/apache/commons/math4/distribution/GammaDistributionTest.java
+++ b/src/test/java/org/apache/commons/math4/distribution/GammaDistributionTest.java
@@ -348,7 +348,7 @@ public class GammaDistributionTest extends RealDistributionAbstractTest {
@Test
public void testMath753Shape142() throws IOException {
- doTestMath753(142.0, 0.5, 1.5, 40.0, 40.0, "gamma-distribution-shape-142.csv");
+ doTestMath753(142.0, 3.3, 1.6, 40.0, 40.0, "gamma-distribution-shape-142.csv");
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0c0455fd/src/test/java/org/apache/commons/math4/optim/nonlinear/scalar/noderiv/BOBYQAOptimizerTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/optim/nonlinear/scalar/noderiv/BOBYQAOptimizerTest.java b/src/test/java/org/apache/commons/math4/optim/nonlinear/scalar/noderiv/BOBYQAOptimizerTest.java
index 2bb5cb6..5064f55 100644
--- a/src/test/java/org/apache/commons/math4/optim/nonlinear/scalar/noderiv/BOBYQAOptimizerTest.java
+++ b/src/test/java/org/apache/commons/math4/optim/nonlinear/scalar/noderiv/BOBYQAOptimizerTest.java
@@ -189,7 +189,7 @@ public class BOBYQAOptimizerTest {
new PointValuePair(point(DIM/2,0.0),0.0);
doTest(new DiffPow(), startPoint, boundaries,
GoalType.MINIMIZE,
- 1e-8, 1e-1, 12000, expected);
+ 1e-8, 1e-1, 21000, expected);
}
@Test
http://git-wip-us.apache.org/repos/asf/commons-math/blob/0c0455fd/src/test/java/org/apache/commons/math4/util/FastMathTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math4/util/FastMathTest.java b/src/test/java/org/apache/commons/math4/util/FastMathTest.java
index 1a2bca4..b95deba 100644
--- a/src/test/java/org/apache/commons/math4/util/FastMathTest.java
+++ b/src/test/java/org/apache/commons/math4/util/FastMathTest.java
@@ -321,9 +321,9 @@ public class FastMathTest {
@Test
public void testLogSpecialCases() {
- Assert.assertTrue("Log of zero should be -Inf", Double.isInfinite(FastMath.log(0.0)));
+ Assert.assertEquals("Log of zero should be -Inf", Double.NEGATIVE_INFINITY, FastMath.log(0.0), 1.0);
- Assert.assertTrue("Log of -zero should be -Inf", Double.isInfinite(FastMath.log(-0.0)));
+ Assert.assertEquals("Log of -zero should be -Inf", Double.NEGATIVE_INFINITY, FastMath.log(-0.0), 1.0);
Assert.assertTrue("Log of NaN should be NaN", Double.isNaN(FastMath.log(Double.NaN)));
@@ -331,8 +331,9 @@ public class FastMathTest {
Assert.assertEquals("Log of Double.MIN_VALUE should be -744.4400719213812", -744.4400719213812, FastMath.log(Double.MIN_VALUE), Precision.EPSILON);
- Assert.assertTrue("Log of infinity should be infinity", Double.isInfinite(FastMath.log(Double.POSITIVE_INFINITY)));
+ Assert.assertEquals("Log of infinity should be infinity", Double.POSITIVE_INFINITY, FastMath.log(Double.POSITIVE_INFINITY), 1.0);
}
+
@Test
public void testExpSpecialCases() {
@@ -343,7 +344,7 @@ public class FastMathTest {
Assert.assertTrue("exp of NaN should be NaN", Double.isNaN(FastMath.exp(Double.NaN)));
- Assert.assertTrue("exp of infinity should be infinity", Double.isInfinite(FastMath.exp(Double.POSITIVE_INFINITY)));
+ Assert.assertEquals("exp of infinity should be infinity", Double.POSITIVE_INFINITY, FastMath.exp(Double.POSITIVE_INFINITY), 1.0);
Assert.assertEquals("exp of -infinity should be 0.0", 0.0, FastMath.exp(Double.NEGATIVE_INFINITY), Precision.EPSILON);
@@ -365,9 +366,9 @@ public class FastMathTest {
Assert.assertTrue("pow(NaN, PI) should be NaN", Double.isNaN(FastMath.pow(Double.NaN, Math.PI)));
- Assert.assertTrue("pow(2.0, Infinity) should be Infinity", Double.isInfinite(FastMath.pow(2.0, Double.POSITIVE_INFINITY)));
+ Assert.assertEquals("pow(2.0, Infinity) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(2.0, Double.POSITIVE_INFINITY), 1.0);
- Assert.assertTrue("pow(0.5, -Infinity) should be Infinity", Double.isInfinite(FastMath.pow(0.5, Double.NEGATIVE_INFINITY)));
+ Assert.assertEquals("pow(0.5, -Infinity) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(0.5, Double.NEGATIVE_INFINITY), 1.0);
Assert.assertEquals("pow(0.5, Infinity) should be 0.0", 0.0, FastMath.pow(0.5, Double.POSITIVE_INFINITY), Precision.EPSILON);
@@ -377,23 +378,25 @@ public class FastMathTest {
Assert.assertEquals("pow(Infinity, -0.5) should be 0.0", 0.0, FastMath.pow(Double.POSITIVE_INFINITY, -0.5), Precision.EPSILON);
- Assert.assertTrue("pow(0.0, -0.5) should be Inf", Double.isInfinite(FastMath.pow(0.0, -0.5)));
+ Assert.assertEquals("pow(0.0, -0.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(0.0, -0.5), 1.0);
- Assert.assertTrue("pow(Inf, 0.5) should be Inf", Double.isInfinite(FastMath.pow(Double.POSITIVE_INFINITY, 0.5)));
+ Assert.assertEquals("pow(Inf, 0.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(Double.POSITIVE_INFINITY, 0.5), 1.0);
- Assert.assertTrue("pow(-0.0, -3.0) should be -Inf", Double.isInfinite(FastMath.pow(-0.0, -3.0)));
+ Assert.assertEquals("pow(-0.0, -3.0) should be -Inf", Double.NEGATIVE_INFINITY, FastMath.pow(-0.0, -3.0), 1.0);
Assert.assertEquals("pow(-0.0, Infinity) should be 0.0", 0.0, FastMath.pow(-0.0, Double.POSITIVE_INFINITY), Precision.EPSILON);
Assert.assertTrue("pow(-0.0, NaN) should be NaN", Double.isNaN(FastMath.pow(-0.0, Double.NaN)));
- Assert.assertTrue("pow(-0.0, -tiny) should be Infinity", Double.isInfinite(FastMath.pow(-0.0, -Double.MIN_VALUE)));
+ Assert.assertEquals("pow(-0.0, -tiny) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(-0.0, -Double.MIN_VALUE), 1.0);
+
+ Assert.assertEquals("pow(-0.0, -huge) should be Infinity", Double.POSITIVE_INFINITY, FastMath.pow(-0.0, -Double.MAX_VALUE), 1.0);
- Assert.assertTrue("pow(-Inf, -3.0) should be -Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 3.0)));
+ Assert.assertEquals("pow(-Inf, -3.0) should be -Inf", Double.NEGATIVE_INFINITY, FastMath.pow(Double.NEGATIVE_INFINITY, 3.0), 1.0);
- Assert.assertTrue("pow(-0.0, -3.5) should be Inf", Double.isInfinite(FastMath.pow(-0.0, -3.5)));
+ Assert.assertEquals("pow(-0.0, -3.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(-0.0, -3.5), 1.0);
- Assert.assertTrue("pow(Inf, 3.5) should be Inf", Double.isInfinite(FastMath.pow(Double.POSITIVE_INFINITY, 3.5)));
+ Assert.assertEquals("pow(Inf, 3.5) should be Inf", Double.POSITIVE_INFINITY, FastMath.pow(Double.POSITIVE_INFINITY, 3.5), 1.0);
Assert.assertEquals("pow(-2.0, 3.0) should be -8.0", -8.0, FastMath.pow(-2.0, 3.0), Precision.EPSILON);
@@ -409,6 +412,16 @@ public class FastMathTest {
Assert.assertTrue("pow(-huge, huge) should be +Inf", Double.isInfinite(FastMath.pow(-Double.MAX_VALUE, Double.MAX_VALUE)));
+ Assert.assertTrue("pow(NaN, -Infinity) should be NaN", Double.isNaN(FastMath.pow(Double.NaN, Double.NEGATIVE_INFINITY)));
+
+ Assert.assertEquals("pow(NaN, 0.0) should be 1.0", 1.0, FastMath.pow(Double.NaN, 0.0), Precision.EPSILON);
+
+ Assert.assertEquals("pow(-Infinity, -Infinity) should be 0.0", 0.0, FastMath.pow(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY), Precision.EPSILON);
+
+ Assert.assertEquals("pow(-huge, -huge) should be 0.0", 0.0, FastMath.pow(-Double.MAX_VALUE, -Double.MAX_VALUE), Precision.EPSILON);
+
+ Assert.assertEquals("pow(-huge, huge) should be +Inf", Double.POSITIVE_INFINITY, FastMath.pow(-Double.MAX_VALUE, Double.MAX_VALUE), 1.0);
+
// Added tests for a 100% coverage
Assert.assertTrue("pow(+Inf, NaN) should be NaN", Double.isNaN(FastMath.pow(Double.POSITIVE_INFINITY, Double.NaN)));
@@ -421,15 +434,26 @@ public class FastMathTest {
Assert.assertEquals("pow(-Inf, -2.0) should be 0.0", 0.0, FastMath.pow(Double.NEGATIVE_INFINITY, -2.0), Precision.EPSILON);
- Assert.assertTrue("pow(-Inf, 1.0) should be -Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 1.0)));
+ Assert.assertEquals("pow(-Inf, 1.0) should be -Inf", Double.NEGATIVE_INFINITY, FastMath.pow(Double.NEGATIVE_INFINITY, 1.0), 1.0);
- Assert.assertTrue("pow(-Inf, 2.0) should be +Inf", Double.isInfinite(FastMath.pow(Double.NEGATIVE_INFINITY, 2.0)));
+ Assert.assertEquals("pow(-Inf, 2.0) should be +Inf", Double.POSITIVE_INFINITY, FastMath.pow(Double.NEGATIVE_INFINITY, 2.0), 1.0);
Assert.assertTrue("pow(1.0, -Inf) should be NaN", Double.isNaN(FastMath.pow(1.0, Double.NEGATIVE_INFINITY)));
}
@Test
+ public void testPowLargeIntegralDouble() {
+ double y = FastMath.scalb(1.0, 65);
+ Assert.assertEquals(Double.POSITIVE_INFINITY, FastMath.pow(FastMath.nextUp(1.0), y), 1.0);
+ Assert.assertEquals(1.0, FastMath.pow(1.0, y), 1.0);
+ Assert.assertEquals(0.0, FastMath.pow(FastMath.nextDown(1.0), y), 1.0);
+ Assert.assertEquals(0.0, FastMath.pow(FastMath.nextUp(-1.0), y), 1.0);
+ Assert.assertEquals(1.0, FastMath.pow(-1.0, y), 1.0);
+ Assert.assertEquals(Double.POSITIVE_INFINITY, FastMath.pow(FastMath.nextDown(-1.0), y), 1.0);
+ }
+
+ @Test
public void testAtan2SpecialCases() {
Assert.assertTrue("atan2(NaN, 0.0) should be NaN", Double.isNaN(FastMath.atan2(Double.NaN, 0.0)));
[2/3] [math] updated ODE userguide documentation.
Posted by lu...@apache.org.
updated ODE userguide documentation.
JIRA: MATH-1225
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/9b6a649f
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/9b6a649f
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/9b6a649f
Branch: refs/heads/master
Commit: 9b6a649f9f9a35e28d8e95b69dc4261fdc2d03d4
Parents: 15a24dc
Author: Luc Maisonobe <lu...@apache.org>
Authored: Sat May 16 18:31:18 2015 +0200
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Sun May 17 15:07:31 2015 +0200
----------------------------------------------------------------------
src/site/xdoc/userguide/ode.xml | 243 ++++++++++++++++++++---------------
1 file changed, 142 insertions(+), 101 deletions(-)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/9b6a649f/src/site/xdoc/userguide/ode.xml
----------------------------------------------------------------------
diff --git a/src/site/xdoc/userguide/ode.xml b/src/site/xdoc/userguide/ode.xml
index 37b266c..d3b2eea 100644
--- a/src/site/xdoc/userguide/ode.xml
+++ b/src/site/xdoc/userguide/ode.xml
@@ -62,6 +62,13 @@
problem, integration range, and step size or error control settings.
</p>
<p>
+ All integrators support expanding the main ODE with one or more secondary ODE to manage
+ additional state that will be integrated together with the main state. This can be used
+ for example to integrate variational equations and compute not only the main state but also
+ its partial derivatives with respect to either the initial state or some parameters, these
+ derivatives being handled be secondary ODE (see below for an example).
+ </p>
+ <p>
The user should describe his problem in his own classes which should implement the
<a href="../apidocs/org/apache/commons/math4/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
interface. Then he should pass it to the integrator he prefers among all the classes that implement
@@ -159,10 +166,11 @@ integrator.addStepHandler(stepHandler);
<p>
Some integrators (the simple ones) use fixed steps that are set at creation time. The more efficient
integrators use variable steps that are handled internally in order to control the integration error
- with respect to a specified accuracy (these integrators extend the
+ of the main state with respect to a specified accuracy (these integrators extend the
<a href="../apidocs/org/apache/commons/math4/ode/AdaptiveStepsizeIntegrator.html">AdaptiveStepsizeIntegrator</a>
- abstract class). In this case, the step handler which is called after each successful step shows up
- the variable stepsize. The <a href="../apidocs/org/apache/commons/math4/ode/sampling/StepNormalizer.html">StepNormalizer</a>
+ abstract class). The secondary equations are explicitly ignored for step size control, in order to get reproducible
+ results regardless of the secondary equations being integrated or not. The step handler which is called after each
+ successful step shows up the variable stepsize. The <a href="../apidocs/org/apache/commons/math4/ode/sampling/StepNormalizer.html">StepNormalizer</a>
class can be used to convert the variable stepsize into a fixed stepsize that can be handled by classes
implementing the <a href="../apidocs/org/apache/commons/math4/ode/sampling/FixedStepHandler.html">FixedStepHandler</a>
interface. Adaptive stepsize integrators can automatically compute the initial stepsize by themselves,
@@ -282,111 +290,94 @@ public int eventOccurred(double t, double[] y, boolean increasing) {
</subsection>
<subsection name="13.5 Derivatives" href="derivatives">
<p>
- If in addition to state y(t) the user needs to compute the sensitivity of the state to
- the initial state or some parameter of the ODE, he will use the classes and interfaces
- from the <a
- href="../apidocs/org/apache/commons/math4/ode/jacobians/package-summary.html">org.apache.commons.ode.jacobians</a>
- package instead of the top level ode package. These classes compute the jacobians
- dy(t)/dy<sub>0</sub> and dy(t)/dp where y<sub>0</sub> is the initial state
- and p is some ODE parameter.
+ If in addition to state y(t) the user needs to compute the sensitivity of the final state with respect to
+ the initial state (dy/dy<sub>0</sub>) or the sensitivity of the final state with respect to some parameters
+ of the ODE (dy/dp<sub>k</sub>), he needs to register the variational equations as a set of secondary equations
+ appended to the main state before the integration starts. Then the integration will propagate the compound
+ state composed of both the main state and its partial derivatives. At the end of the integration, the Jacobian
+ matrices are extracted from the integrated secondary state. The <a
+ href="../apidocs/org/apache/commons/math4/ode/JacobianMatrices.html">JacobianMatrices</a> class can do most of
+ this as long as the local derivatives are provided to it. It will set up the variational equations, register
+ them as secondary equations into the ODE, and it will set up the initial values and retrieve the intermediate
+ and finale values as Jacobian matrices.
</p>
<p>
- The classes and interfaces in this package mimic the behavior of the classes and
- interfaces of the top level ode package, only adding parameters arrays for the jacobians.
- The behavior of these classes is to create a compound state vector z containing both
- the state y(t) and its derivatives dy(t)/dy<sub>0</sub> and dy(t)/dp and
- to set up an extended problem by adding the equations for the jacobians automatically.
- These extended state and problems are then provided to a classical underlying integrator
- chosen by user.
+ If for example the original state dimension is 6 and there are 3 parameters, the compound state will be a 60
+ elements array. The first 6 elements will be the original state, the next 36 elements will represent the 6x6
+ Jacobian matrix of the final state with respect to the initial state, and the remaining 18 elements will
+ represent the 6x3 Jacobian matrix of the final state with respect to the 3 parameters. The <a
+ href="../apidocs/org/apache/commons/math4/ode/JacobianMatrices.html">JacobianMatrices</a> class does the mapping
+ between the 60 elements compound state and the Jacobian matrices and sets up the correcsponding secondary equations.
</p>
<p>
- This behavior imply there will be a top level integrator knowing about state and jacobians
- and a low level integrator knowing only about compound state (which may be big). If the user
- wants to deal with the top level only, he will use the specialized step handler and event
- handler classes registered at top level. He can also register classical step handlers and
- event handlers, but in this case will see the big compound state. This state is guaranteed
- to contain the original state in the first elements, followed by the jacobian with respect
- to initial state (in row order), followed by the jacobian with respect to parameters (in
- row order). If for example the original state dimension is 6 and there are 3 parameters,
- the compound state will be a 60 elements array. The first 6 elements will be the original
- state, the next 36 elements will be the jacobian with respect to initial state, and the
- remaining 18 will be the jacobian with respect to parameters. Dealing with low level
- step handlers and event handlers is cumbersome if one really needs the jacobians in these
- methods, but it also prevents many data being copied back and forth between state and
- jacobians on one side and compound state on the other side.
+ As the variational equations are considered to be secondary equations here, variable step integrators ignore
+ them for step size control: they rely only on the main state. This feature is a design choice. The rationale is
+ to get exactly the same steps, regardless of the Jacobians being computed or not, hence ensuring reproducible
+ results in both cases.
</p>
<p>
- In order to compute dy(t)/dy<sub>0</sub> and dy(t/dp for any t, the algorithm
- needs not only the ODE function f such that y'=f(t,y) but also its local jacobians
- df(t, y, p)/dy and df(t, y, p)/dp.
+ What remains of user responsibility is to provide the local Jacobians df(t, y, p)/dy and df(t, y, p)/dp<sub>k</sub>
+ corresponding the the main ODE y'=f(t, y, p). The main ODE is as usual provided by the user as a class implementing
+ the <a href="../apidocs/org/apache/commons/math4/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
+ interface or a sub-interface.
</p>
<p>
- If the function f is too complex, the user can simply rely on internal differentiation
- using finite differences to compute these local jacobians. So rather than the <a
+ If the ODE is simple enough that the user can implement df(t, y, p)/dy directly, then instead of providing an
+ implementation of the <a
+ href="../apidocs/org/apache/commons/math4/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
+ interface only, the user should rather provide an implementation of the <a
+ href="../apidocs/org/apache/commons/math4/ode/MainStateJacobianProvider.html">MainStateJacobianProvider</a> interface,
+ which extends the previous interface by adding a method to compute df(t, y, p)/dy. The user class is used as a
+ constructor parameter of the <a href="../apidocs/org/apache/commons/math4/ode/JacobianMatrices.html">JacobianMatrices</a>
+ class. If the ODE is too complex or the user simply does not bother implementing df(t, y, p)/dy directly, then
+ the ODE can still be implemented using the simple <a
href="../apidocs/org/apache/commons/math4/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
- interface he will implement the <a
- href="../apidocs/org/apache/commons/math4/ode/jacobians/ParameterizedODE.html">ParameterizedODE</a>
- interface. Considering again our example where only ω is considered a parameter, we get:
+ interface and given as such to another constructor of the <a
+ href="../apidocs/org/apache/commons/math4/ode/JacobianMatrices.html">JacobianMatrices</a> class, but in this case an array
+ hy must also be provided that will contain the step size to use form each component of the main state vector y, and
+ the Jacobian f(t, y, p)/dy will be computed internally using finite differences. This will of course trigger more evaluations
+ of the ODE at each step and will suffer from finite differences errors, but it is much simpler to implement from a user
+ point of view.
</p>
- <source>
-public class BasicCircleODE implements ParameterizedODE {
-
- private double[] c;
- private double omega;
-
- public BasicCircleODE(double[] c, double omega) {
- this.c = c;
- this.omega = omega;
- }
-
- public int getDimension() {
- return 2;
- }
-
- public void computeDerivatives(double t, double[] y, double[] yDot) {
- yDot[0] = omega * (c[1] - y[1]);
- yDot[1] = omega * (y[0] - c[0]);
- }
-
- public int getParametersDimension() {
- // we are only interested in the omega parameter
- return 1;
- }
-
- public void setParameter(int i, double value) {
- omega = value;
- }
-
-}
- </source>
<p>
- This ODE is provided to the specialized integrator with two arrays specifying the
- step sizes to use for finite differences (one array for derivation with respect to
- state y, one array for derivation with respect to parameters p):
+ The parameters are identified by a name (a simple user defined string), which are also specified at <a
+ href="../apidocs/org/apache/commons/math4/ode/JacobianMatrices.html">JacobianMatrices</a> class construction. If the ODE
+ is simple enough that the user can implement df(t, y, p)/dp<sub>k</sub> directly for some of the parameters p<sub>k</sub>,
+ then he can provide one or more classes implementing the <a
+ href="../apidocs/org/apache/commons/math4/ode/ParameterJacobianProvider.html">ParameterJacobianProvider</a> interface by
+ calling the JacobianMatrices.addParameterJacobianProvide method. The parameters are handled one at a time, but all the calls to
+ ParameterJacobianProvider.computeParameterJacobian will be grouped in one sequence after the call to MainStateJacobianProvider.computeMainStateJacobian
+ This feature can be used when all the derivatives share a lot of costly computation. In this case, the user
+ is advised to compute all the needed derivatives at once during the call to computeMainStateJacobian, including the
+ partial derivatives with respect to the parameters and to store the derivatives temporary. Then when the next calls to
+ computeParameterJacobian will be triggerred, it will be sufficient to return the already computed derivatives. With this
+ architecture, many computation can be saved. This of course implies that the classes implementing both interfaces know
+ each other (they can even be the same class if desired, but it is not required). If the ODE is too complex or the user
+ simply does not bother implementing df(t, y, p)/dp<sub>k</sub> directly for some k, then
+ the JacobianMatrices.setParameterStep method should be called so finite differences are used to compute the derivatives
+ for this parameter. It is possible to have some parameters for which derivatives are provided by a direct implementation
+ while other parameters are computed using finite differences during the same integration.
</p>
- <source>
-double[] hY = new double[] { 0.001, 0.001 };
-double[] hP = new double[] { 1.0e-6 };
-FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode, hY, hP);
-integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
- </source>
<p>
- If the function f is simple, the user can simply provide the local jacobians
- by himself. So rather than the <a
- href="../apidocs/org/apache/commons/math4/ode/FirstOrderDifferentialEquations.html">FirstOrderDifferentialEquations</a>
- interface he will implement the <a
- href="../apidocs/org/apache/commons/math4/ode/jacobians/ODEWithJacobians.html">ODEWithJacobians</a>
- interface. Considering again our example where only ω is considered a parameter, we get:
+ The following example corresponds to a simple case where all derivatives can be computed analytically. The state is
+ a 2D point travelling along a circle. There are three parameters : the two coordinates of the center and the
+ angular velocity.
</p>
<source>
-public class EnhancedCircleODE implements ODEWithJacobians {
+public class CircleODE implements MainStateJacobianProvider, ParameterJacobianProvider {
+
+ public static final String CENTER_X = "cx";
+ public static final String CENTER_Y = "cy";
+ public static final String OMEGA = "omega";
private double[] c;
private double omega;
+ private double[][] savedDfDp;
- public EnhancedCircleODE(double[] c, double omega) {
+ public CircleODE(double[] c, double omega) {
this.c = c;
this.omega = omega;
+ this.savedDfDp = new double[2][3];
}
public int getDimension() {
@@ -394,39 +385,89 @@ public class EnhancedCircleODE implements ODEWithJacobians {
}
public void computeDerivatives(double t, double[] y, double[] yDot) {
+ // the state is a 2D point, the ODE therefore corresponds to the velocity
yDot[0] = omega * (c[1] - y[1]);
yDot[1] = omega * (y[0] - c[0]);
}
- public int getParametersDimension() {
- // we are only interested in the omega parameter
- return 1;
+ public Collection<String> getParametersNames() {
+ return Arrays.asList(CENTER_X, CENTER_Y, OMEGA);
}
- public void computeJacobians(double t, double[] y, double[] yDot, double[][] dFdY, double[][] dFdP) {
+ public boolean isSupported(String name) {
+ return CENTER_X.equals(name) || CENTER_Y.equals(name) || OMEGA.equals(name);
+ }
+
+ public void computeMainStateJacobian(double t, double[] y, double[] yDot, double[][] dFdY) {
+ // compute the Jacobian of the main state
dFdY[0][0] = 0;
dFdY[0][1] = -omega;
dFdY[1][0] = omega;
dFdY[1][1] = 0;
- dFdP[0][0] = 0;
- dFdP[0][1] = omega;
- dFdP[0][2] = c[1] - y[1];
- dFdP[1][0] = -omega;
- dFdP[1][1] = 0;
- dFdP[1][2] = y[0] - c[0];
-
+ // precompute the derivatives with respect to the parameters,
+ // they will be provided back when computeParameterJacobian are called later on
+ savedDfDp[0][0] = 0;
+ savedDfDp[0][1] = omega;
+ savedDfDp[0][2] = c[1] - y[1];
+ savedDfDp[1][0] = -omega;
+ savedDfDp[1][1] = 0;
+ savedDfDp[1][2] = y[0] - c[0];
+
}
+ public void computeParameterJacobian(double t, double[] y, double[] yDot,
+ String paramName, double[] dFdP) {
+ // we simply return the derivatives precomputed earlier
+ if (CENTER_X.equals(paramName)) {
+ dFdP[0] = savedDfDp[0][0];
+ dFdP[1] = savedDfDp[1][0];
+ } else if (CENTER_Y.equals(paramName)) {
+ dFdP[0] = savedDfDp[0][1];
+ dFdP[1] = savedDfDp[1][1];
+ } else {
+ dFdP[0] = savedDfDp[0][2];
+ dFdP[1] = savedDfDp[1][2];
+ }
+ }
+
}
</source>
<p>
- This ODE is provided to the specialized integrator as is:
+ This ODE is integrated as follows:
</p>
<source>
-FirstOrderIntegratorWithJacobians integrator = new FirstOrderIntegratorWithJacobians(dp853, ode);
-integrator.integrate(t0, y0, dy0dp, t, y, dydy0, dydp);
+ CircleODE circle = new CircleODE(new double[] {1.0, 1.0 }, 0.1);
+
+ // here, we could select only a subset of the parameters, or use another order
+ JacobianMatrices jm = new JacobianMatrices(circle, CircleODE.CENTER_X, CircleODE.CENTER_Y, CircleODE.OMEGA);
+ jm.addParameterJacobianProvider(circle);
+
+ ExpandableStatefulODE efode = new ExpandableStatefulODE(circle);
+ efode.setTime(0);
+ double[] y = { 1.0, 0.0 };
+ efode.setPrimaryState(y);
+
+ // create the variational equations and append them to the main equations, as secondary equations
+ jm.registerVariationalEquations(efode);
+
+ // integrate the compound state, with both main and additional equations
+ DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0e3, 1.0e-10, 1.0e-12);
+ integrator.setMaxEvaluations(5000);
+ integrator.integrate(efode, 20.0);
+
+ // retrieve the Jacobian of the final state with respect to initial state
+ double[][] dYdY0 = new double[2][2];
+ jm.getCurrentMainSetJacobian(dYdY0);
+
+ // retrieve the Jacobians of the final state with resepct to the various parameters
+ double[] dYdCx = new double[2];
+ double[] dYdCy = new double[2];
+ double[] dYdOm = new double[2];
+ jm.getCurrentParameterJacobian(CircleODE.CENTER_X, dYdCx);
+ jm.getCurrentParameterJacobian(CircleODE.CENTER_Y, dYdCy);
+ jm.getCurrentParameterJacobian(CircleODE.OMEGA, dYdOm);
</source>
</subsection>
</section>