Apache
/
commons-math
mirror of https://github.com/apache/commons-math.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Added a new package dealing with differentials. 

The package is intended to deals with one or more free parameters and
derivation order 1 or higher.

The core elements are based on Dan Kalman paper "Recursive Multivariate
Automatic Differentiation", Mathematics Magazine, vol. 75, no. 3, June
2002. For efficiency, the recursive structure is compiled as simple
loops once for each pair (number of free parameters, derivation order).

This is work in progress, there are still some features missing even in
the most basic blocks (typically the asin, acos, atan, atant2 and taylor
methods in DSCompiler). There are also still no high level
differentiator implementation.

git-svn-id: https://svn.apache.org/repos/asf/commons/proper/math/trunk@1370951 13f79535-47bb-0310-9956-ffa450edef68
This commit is contained in:
Luc Maisonobe 2012-08-08 20:33:43 +00:00
parent 13f4032d6c
commit 87b597c622
9 changed files with 3005 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/changes/changes.xml
View File

@ -52,6 +52,10 @@ If the output is not quite correct, check for invisible trailing spaces!
<body>
<release version="3.1" date="TBD" description="
">
<action dev="luc" type="add" >
Added a new package dealing with differentials, for one or more free
parameters and derivation order 1 or higher.
</action>
<action dev="tn" type="add" issue="MATH-777" due-to="Reid Hochstedler">
Added additional crossover policies: "CycleCrossover", "NPointCrossover",
"OrderedCrossover" and "UniformCrossover".

1246
src/main/java/org/apache/commons/math3/analysis/differentiation/DSCompiler.java Normal file
View File

File diff suppressed because it is too large Load Diff

576
src/main/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructure.java Normal file
View File

@ -0,0 +1,576 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import java.io.Serializable;
import org.apache.commons.math3.Field;
import org.apache.commons.math3.FieldElement;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
/** Class representing both the value and the differentials of a function.
* <p>This class is the workhorse of the differentiation package.</p>
* <p>This class is an implementation of the extension to Rall's
* numbers described in Dan Kalman's paper <a
* href="http://www.math.american.edu/People/kalman/pdffiles/mmgautodiff.pdf">Doubly
* Recursive Multivariate Automatic Differentiation</a>, Mathematics Magazine, vol. 75,
* no. 3, June 2002.</p>. Rall's numbers are an extension to the real numbers used
* throughout mathematical expressions; they hold the derivative together with the
* value of a function. Dan Kalman's derivative structures holds all partial derivatives
* up to any specified order, with respect to any number of free variables. Rall's
* number therefore can be seen as derivative structures for order one derivative and
* one free variable, and real numbers can be seen as derivative structures with zero
* order derivative and no free variables.</p>
* <p>{@link DerivativeStructure} instances can be used directly thanks to
* the arithmetic operators to the mathematical functions provided as static
* methods by this class (+, -, *, /, %, sin, cos ...).</p>
* <p>Implementing complex expressions by hand using these classes is
* however a complex and error-prone task, so the classical use is
* simply to develop computation code using standard primitive double
* values and to use {@link UnivariateDifferentiator differentiators} to create
* the {@link DerivativeStructure}-based instances.</p>
* <p>Instances of this class are guaranteed to be immutable.</p>
* @see DSCompiler
* @version $Id$
* @since 3.1
*/
public class DerivativeStructure implements FieldElement<DerivativeStructure>, Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20120730L;
/** Compiler for the current dimensions. */
private transient DSCompiler compiler;
/** Combined array holding all values. */
private final double[] data;
/** Build an instance with all values and derivatives set to 0.
* @param compiler compiler to use for computation
*/
private DerivativeStructure(final DSCompiler compiler) {
this.compiler = compiler;
this.data = new double[compiler.getSize()];
}
/** Build an instance with all values and derivatives set to 0.
* @param variables number of variables
* @param order derivation order
*/
public DerivativeStructure(final int variables, final int order) {
this(DSCompiler.getCompiler(variables, order));
}
/** Build an instance representing a constant value.
* @param variables number of variables
* @param order derivation order
* @param value value of the constant
* @see #DerivativeStructure(int, int, int, double)
*/
public DerivativeStructure(final int variables, final int order, final double value) {
this(variables, order);
this.data[0] = value;
}
/** Build an instance representing a variable.
* <p>Instances built using this constructor are considered
* to be the free variables with respect to which differentials
* are computed. As such, their differential with respect to
* themselves is +1.</p>
* @param variables number of variables
* @param order derivation order
* @param index index of the variable (from 0 to {@code variables - 1})
* @param value value of the variable
* @exception NumberIsTooLargeException if index is equal to variables or larger
* @see #DerivativeStructure(int, int, double)
*/
public DerivativeStructure(final int variables, final int order,
final int index, final double value)
throws NumberIsTooLargeException {
this(variables, order, value);
if (index >= variables) {
throw new NumberIsTooLargeException(index, variables, false);
}
if (order > 0) {
// the derivative of the variable with respect to itself is 1.0
data[DSCompiler.getCompiler(index, order).getSize()] = 1.0;
}
}
/** Linear combination constructor.
* The derivative structure built will be a1 * ds1 + a2 * ds2
* @param a1 first scale factor
* @param ds1 first base (unscaled) derivative structure
* @param a2 second scale factor
* @param ds2 second base (unscaled) derivative structure
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure(final double a1, final DerivativeStructure ds1,
final double a2, final DerivativeStructure ds2)
throws DimensionMismatchException {
this(ds1.compiler);
compiler.checkCompatibility(ds2.compiler);
compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0, data, 0);
}
/** Linear combination constructor.
* The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3
* @param a1 first scale factor
* @param ds1 first base (unscaled) derivative structure
* @param a2 second scale factor
* @param ds2 second base (unscaled) derivative structure
* @param a3 third scale factor
* @param ds3 third base (unscaled) derivative structure
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure(final double a1, final DerivativeStructure ds1,
final double a2, final DerivativeStructure ds2,
final double a3, final DerivativeStructure ds3)
throws DimensionMismatchException {
this(ds1.compiler);
compiler.checkCompatibility(ds2.compiler);
compiler.checkCompatibility(ds3.compiler);
compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0, a3, ds3.data, 0, data, 0);
}
/** Linear combination constructor.
* The derivative structure built will be a1 * ds1 + a2 * ds2 + a3 * ds3 + a4 * ds4
* @param a1 first scale factor
* @param ds1 first base (unscaled) derivative structure
* @param a2 second scale factor
* @param ds2 second base (unscaled) derivative structure
* @param a3 third scale factor
* @param ds3 third base (unscaled) derivative structure
* @param a4 fourth scale factor
* @param ds4 fourth base (unscaled) derivative structure
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public 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)
throws DimensionMismatchException {
this(ds1.compiler);
compiler.checkCompatibility(ds2.compiler);
compiler.checkCompatibility(ds3.compiler);
compiler.checkCompatibility(ds4.compiler);
compiler.linearCombination(a1, ds1.data, 0, a2, ds2.data, 0,
a3, ds3.data, 0, a4, ds4.data, 0,
data, 0);
}
/** Copy constructor.
* @param instance to copy
*/
private DerivativeStructure(final DerivativeStructure ds) {
this.compiler = ds.compiler;
this.data = ds.data.clone();
}
/** Get the number of free parameters.
* @return number of free parameters
*/
public int getFreeParameters() {
return compiler.getFreeParameters();
}
/** Get the derivation order.
* @return derivation order
*/
public int getOrder() {
return compiler.getOrder();
}
/** Get the value part of the derivative structure.
* @return value part of the derivative structure
* @see #getPartialDerivative(int...)
*/
public double getValue() {
return data[0];
}
/** Get a partial derivative.
* @param orders derivation orders with respect to each variable (if all orders are 0,
* the value is returned)
* @return partial derivative
* @see #getValue()
* @exception DimensionMismatchException if the numbers of variables does not
* match the instance
* @exception NumberIsTooLargeException if sum of derivation orders is larger
* than the instance limits
*/
public double getPartialDerivative(final int ... orders)
throws DimensionMismatchException, NumberIsTooLargeException {
return data[compiler.getPartialDerivativeIndex(orders)];
}
/** '+' operator.
* @param a right hand side parameter of the operator
* @return this+a
*/
public DerivativeStructure add(final double a) {
final DerivativeStructure ds = new DerivativeStructure(this);
ds.data[0] += a;
return ds;
}
/** '+' operator.
* @param a right hand side parameter of the operator
* @return this+a
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure add(final DerivativeStructure a)
throws DimensionMismatchException {
compiler.checkCompatibility(a.compiler);
final DerivativeStructure ds = new DerivativeStructure(this);
compiler.add(data, 0, a.data, 0, ds.data, 0);
return ds;
}
/** '-' operator.
* @param a right hand side parameter of the operator
* @return this-a
*/
public DerivativeStructure subtract(final double a) {
return add(-a);
}
/** '-' operator.
* @param a right hand side parameter of the operator
* @return this-a
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure subtract(final DerivativeStructure a)
throws DimensionMismatchException {
compiler.checkCompatibility(a.compiler);
final DerivativeStructure ds = new DerivativeStructure(this);
compiler.subtract(data, 0, a.data, 0, ds.data, 0);
return ds;
}
/** {@inheritDoc} */
public DerivativeStructure multiply(final int n) {
return multiply((double) n);
}
/** '&times;' operator.
* @param a right hand side parameter of the operator
* @return this&times;a
*/
public DerivativeStructure multiply(final double a) {
final DerivativeStructure ds = new DerivativeStructure(this);
for (int i = 0; i < ds.data.length; ++i) {
ds.data[i] *= a;
}
return ds;
}
/** '&times;' operator.
* @param a right hand side parameter of the operator
* @return this&times;a
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure multiply(final DerivativeStructure a)
throws DimensionMismatchException {
compiler.checkCompatibility(a.compiler);
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.multiply(data, 0, a.data, 0, result.data, 0);
return result;
}
/** '&divides;' operator.
* @param a right hand side parameter of the operator
* @return this&divides;a
*/
public DerivativeStructure divide(final double a) {
final DerivativeStructure ds = new DerivativeStructure(this);
for (int i = 0; i < ds.data.length; ++i) {
ds.data[i] /= a;
}
return ds;
}
/** '&divides;' operator.
* @param a right hand side parameter of the operator
* @return this&divides;a
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure divide(final DerivativeStructure a)
throws DimensionMismatchException {
compiler.checkCompatibility(a.compiler);
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.divide(data, 0, a.data, 0, result.data, 0);
return result;
}
/** '%' operator.
* @param a right hand side parameter of the operator
* @return this%a
*/
public DerivativeStructure remainder(final double a) {
final DerivativeStructure ds = new DerivativeStructure(this);
ds.data[0] = ds.data[0] % a;
return ds;
}
/** '%' operator.
* @param a right hand side parameter of the operator
* @return this%a
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public DerivativeStructure remainder(final DerivativeStructure a)
throws DimensionMismatchException {
compiler.checkCompatibility(a.compiler);
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.remainder(data, 0, a.data, 0, result.data, 0);
return result;
}
/** unary '-' operator.
* @return -this
*/
public DerivativeStructure negate() {
final DerivativeStructure ds = new DerivativeStructure(compiler);
for (int i = 0; i < ds.data.length; ++i) {
ds.data[i] = -data[i];
}
return ds;
}
/** {@inheritDoc} */
public DerivativeStructure reciprocal() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.pow(data, 0, -1, result.data, 0);
return result;
}
/** Square root.
* @return square root of the instance
*/
public DerivativeStructure sqrt() {
return rootN(2);
}
/** Cubic root.
* @return cubic root of the instance
*/
public DerivativeStructure cbrt() {
return rootN(3);
}
/** N<sup>th</sup> root.
* @param n order of the root
* @return n<sup>th</sup> root of the instance
*/
public DerivativeStructure rootN(final int n) {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.rootN(data, 0, n, result.data, 0);
return result;
}
/** {@inheritDoc} */
public Field<DerivativeStructure> getField() {
return new Field<DerivativeStructure>() {
/** {@inheritDoc} */
public DerivativeStructure getZero() {
return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 0.0);
}
/** {@inheritDoc} */
public DerivativeStructure getOne() {
return new DerivativeStructure(compiler.getFreeParameters(), compiler.getOrder(), 1.0);
}
/** {@inheritDoc} */
public Class<? extends FieldElement<DerivativeStructure>> getRuntimeClass() {
return DerivativeStructure.class;
}
};
}
/** Power operation.
* @param p power to apply
* @return this<sup>p</sup>
*/
public DerivativeStructure pow(final double p) {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.pow(data, 0, p, result.data, 0);
return result;
}
/** Integer power operation.
* @param n power to apply
* @return this<sup>n</sup>
*/
public DerivativeStructure pow(final int n) {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.pow(data, 0, n, result.data, 0);
return result;
}
/** Exponential.
* @return exponential of the instance
*/
public DerivativeStructure exp() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.exp(data, 0, result.data, 0);
return result;
}
/** Natural logarithm.
* @return logarithm of the instance
*/
public DerivativeStructure log() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.log(data, 0, result.data, 0);
return result;
}
/** Cosine operation.
* @return cos(this)
*/
public DerivativeStructure cos() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.cos(data, 0, result.data, 0);
return result;
}
/** Sine operation.
* @return sin(this)
*/
public DerivativeStructure sin() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.sin(data, 0, result.data, 0);
return result;
}
/** Tangent operation.
* @return tan(this)
*/
public DerivativeStructure tan() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.tan(data, 0, result.data, 0);
return result;
}
/** Arc cosine operation.
* @return acos(this)
*/
public DerivativeStructure acos() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.acos(data, 0, result.data, 0);
return result;
}
/** Arc sine operation.
* @return asin(this)
*/
public DerivativeStructure asin() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.asin(data, 0, result.data, 0);
return result;
}
/** Arc tangent operation.
* @return tan(this)
*/
public DerivativeStructure atan() {
final DerivativeStructure result = new DerivativeStructure(compiler);
compiler.atan(data, 0, result.data, 0);
return result;
}
/** Two arguments arc tangent operation.
* @param y first argument of the arc tangent
* @param x second argument of the arc tangent
* @return atan2(y, x)
* @exception DimensionMismatchException if number of free parameters or orders are inconsistent
*/
public static DerivativeStructure atan2(final DerivativeStructure y, final DerivativeStructure x)
throws DimensionMismatchException {
y.compiler.checkCompatibility(x.compiler);
final DerivativeStructure result = new DerivativeStructure(y.compiler);
y.compiler.atan2(y.data, 0, x.data, 0, result.data, 0);
return result;
}
/** Evaluate Taylor expansion a derivative structure.
* @param offsets parameters offsets (dx, dy, ...)
* @return value of the Taylor expansion at x+dx, y.dy, ...
*/
public double taylor(final double ... offsets) {
return compiler.taylor(data, 0, offsets);
}
/**
* Replace the instance with a data transfer object for serialization.
* @return data transfer object that will be serialized
*/
private Object writeReplace() {
return new DataTransferObject(compiler.getFreeParameters(), compiler.getOrder(), data);
}
/** Internal class used only for serialization. */
private static class DataTransferObject implements Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20120730L;
/** Number of variables.
* @Serial
*/
private final int variables;
/** Derivation order.
* @Serial
*/
private final int order;
/** Partial derivatives.
* @Serial
*/
private final double[] data;
/** Simple constructor.
* @param variables number of variables
* @param order derivation order
* @param data partial derivatives
*/
public DataTransferObject(final int variables, final int order, final double[] data) {
this.variables = variables;
this.order = order;
this.data = data;
}
/** Replace the deserialized data transfer object with a {@link DerivativeStructure}.
* @return replacement {@link DerivativeStructure}
*/
private Object readResolve() {
final DerivativeStructure ds = new DerivativeStructure(variables, order);
System.arraycopy(data, 0, ds.data, 0, ds.data.length);
return ds;
}
}
}

61
src/main/java/org/apache/commons/math3/analysis/differentiation/UnivariateDifferential.java Normal file
View File

@ -0,0 +1,61 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import org.apache.commons.math3.analysis.UnivariateFunction;
/** Interface for univariate functions derivatives.
* <p>This interface represents a simple function which computes
* both the value and the first derivative of a mathematical function.
* The derivative is computed with respect to the input variable.</p>
* @see UnivariateDifferentiable
* @see UnivariateDifferentiator
* @since 3.1
* @version $Id$
*/
public interface UnivariateDifferential {
/** Get the primitive function associated with this differential.
* <p>Each {@link UnivariateDifferential} instance is tightly bound
* to an {@link UnivariateDifferentiable} instance. If the state of
* the primitive instance changes in any way that affects the
* differential computation, this binding allows this change to
* be immediately seen by the derivative instance, there is no need
* to differentiate the primitive again. The existing instance is aware
* of the primitive changes.</p>
* <p>In other words in the following code snippet, the three values
* f1, f2 and f3 should be equal (at least at machine tolerance level)</p>
* <pre>
* UnivariateDifferential derivative = differentiator.differentiate(derivable);
* derivable.someFunctionThatMutatesHeavilyTheInstance();
* double f1 = derivable.f(t);
* double f2 = derivative.getPrimitive().f(t);
* double f3 = derivative.f(new DerivativeStructure(variables, order, index, t)).getValue();
* </pre>
* @return primitive function bound to this derivative
*/
UnivariateFunction getPrimitive();
/** Simple mathematical function.
* <p>{@link UnivariateDifferential} classes compute both the
* value and the first derivative of the function.</p>
* @param t function input value
* @return function result
*/
DerivativeStructure f(DerivativeStructure t);
}

34
src/main/java/org/apache/commons/math3/analysis/differentiation/UnivariateDifferentiator.java Normal file
View File

@ -0,0 +1,34 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import org.apache.commons.math3.analysis.UnivariateFunction;
/** Interface defining the function differentiation operation.
* @version $Id$
* @since 3.1
*/
public interface UnivariateDifferentiator {
/** Create an implementation of a differential for a
* {@link UnivariateDifferentiable differentiable function}.
* @param function function to differentiate
* @return differential function
*/
UnivariateDifferential differentiate(UnivariateFunction function);
}

35
src/main/java/org/apache/commons/math3/analysis/differentiation/package-info.java Normal file
View File

@ -0,0 +1,35 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
/**
*
* <p>
* This package holds the main interfaces and basic building block classes
* dealing with differentiation.
* The core class is {@link DerivativeStructure} which holds the value and
* the differentials of a function. This class handles some arbitrary number
* of free parameters and arbitrary derivation order. It is used both as
* the input and the output type for the {@link UnivariateDifferential}
* interface. Any differentiable function should implement this interface.
* </p>
* <p>
* The {@link UnivariateDifferentiator} interface defines a way to differentation
* a simple {@link org.apache.commons.math3.analysis.UnivariateFunction
* univariate function} and get a {@link differential function}.
* </p>
*
*/
package org.apache.commons.math3.analysis.differentiation;

61
src/main/java/org/apache/commons/math3/util/FastMath.java
View File

@ -18,6 +18,9 @@ package org.apache.commons.math3.util;
import java.io.PrintStream;
import org.apache.commons.math3.exception.NotPositiveException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
/**
* Faster, more accurate, portable alternative to {@link Math} and
* {@link StrictMath} for large scale computation.
@ -1145,7 +1148,7 @@ public class FastMath {
/* Normalize the subnormal number. */
bits <<= 1;
while ( (bits & 0x0010000000000000L) == 0) {
--exp;
exp--;
bits <<= 1;
}
}
@ -1165,9 +1168,8 @@ public class FastMath {
xa = aa;
xb = ab;
final double[] lnCoef_last = LN_QUICK_COEF[LN_QUICK_COEF.length - 1];
double ya = lnCoef_last[0];
double yb = lnCoef_last[1];
double ya = LN_QUICK_COEF[LN_QUICK_COEF.length-1][0];
double yb = LN_QUICK_COEF[LN_QUICK_COEF.length-1][1];
for (int i = LN_QUICK_COEF.length - 2; i >= 0; i--) {
/* Multiply a = y * x */
@ -1179,9 +1181,8 @@ public class FastMath {
yb = aa - ya + ab;
/* Add a = y + lnQuickCoef */
final double[] lnCoef_i = LN_QUICK_COEF[i];
aa = ya + lnCoef_i[0];
ab = yb + lnCoef_i[1];
aa = ya + LN_QUICK_COEF[i][0];
ab = yb + LN_QUICK_COEF[i][1];
/* Split y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
@ -1201,7 +1202,7 @@ public class FastMath {
}
// lnm is a log of a number in the range of 1.0 - 2.0, so 0 <= lnm < ln(2)
final double[] lnm = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
double lnm[] = lnMant.LN_MANT[(int)((bits & 0x000ffc0000000000L) >> 42)];
/*
double epsilon = x / Double.longBitsToDouble(bits & 0xfffffc0000000000L);
@ -1212,7 +1213,7 @@ public class FastMath {
// y is the most significant 10 bits of the mantissa
//double y = Double.longBitsToDouble(bits & 0xfffffc0000000000L);
//double epsilon = (x - y) / y;
final double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
double epsilon = (bits & 0x3ffffffffffL) / (TWO_POWER_52 + (bits & 0x000ffc0000000000L));
double lnza = 0.0;
double lnzb = 0.0;
@ -1226,15 +1227,14 @@ public class FastMath {
double xb = ab;
/* Need a more accurate epsilon, so adjust the division. */
final double numer = bits & 0x3ffffffffffL;
final double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
double numer = bits & 0x3ffffffffffL;
double denom = TWO_POWER_52 + (bits & 0x000ffc0000000000L);
aa = numer - xa*denom - xb * denom;
xb += aa / denom;
/* Remez polynomial evaluation */
final double[] lnCoef_last = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1];
double ya = lnCoef_last[0];
double yb = lnCoef_last[1];
double ya = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][0];
double yb = LN_HI_PREC_COEF[LN_HI_PREC_COEF.length-1][1];
for (int i = LN_HI_PREC_COEF.length - 2; i >= 0; i--) {
/* Multiply a = y * x */
@ -1246,9 +1246,8 @@ public class FastMath {
yb = aa - ya + ab;
/* Add a = y + lnHiPrecCoef */
final double[] lnCoef_i = LN_HI_PREC_COEF[i];
aa = ya + lnCoef_i[0];
ab = yb + lnCoef_i[1];
aa = ya + LN_HI_PREC_COEF[i][0];
ab = yb + LN_HI_PREC_COEF[i][1];
/* Split y = a */
tmp = aa * HEX_40000000;
ya = aa + tmp - tmp;
@ -1581,6 +1580,34 @@ public class FastMath {
}
/**
* Raise a double to an int power.
*
* @param d Number to raise.
* @param e Exponent.
* @return d<sup>e</sup>
*/
public static double pow(double d, int e) {
if (e == 0) {
return 1.0;
} else if (e < 0) {
e = -e;
d = 1.0 / d;
}
double result = 1;
double d2p = d;
while (e != 0) {
if ((e & 0x1) != 0) {
result *= d2p;
}
d2p *= d2p;
e = e >> 1;
}
return result;
}
/**
* Computes sin(x) - x, where |x| < 1/16.
* Use a Remez polynomial approximation.

466
src/test/java/org/apache/commons/math3/analysis/differentiation/DSCompilerTest.java Normal file
View File

@ -0,0 +1,466 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import java.lang.reflect.Field;
import java.util.HashMap;
import java.util.Map;
import org.apache.commons.math3.util.ArithmeticUtils;
import org.junit.Assert;
import org.junit.Test;
/**
* Test for class {@link DSCompiler}.
*/
public class DSCompilerTest {
@Test
public void testSize() {
for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) {
long expected = ArithmeticUtils.binomialCoefficient(i + j, i);
Assert.assertEquals(expected, DSCompiler.getCompiler(i, j).getSize());
Assert.assertEquals(expected, DSCompiler.getCompiler(j, i).getSize());
}
}
}
@Test
public void testIndices() {
DSCompiler c = DSCompiler.getCompiler(0, 0);
checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
c = DSCompiler.getCompiler(0, 1);
checkIndices(c.getPartialDerivativeOrders(0), new int[0]);
c = DSCompiler.getCompiler(1, 0);
checkIndices(c.getPartialDerivativeOrders(0), 0);
c = DSCompiler.getCompiler(1, 1);
checkIndices(c.getPartialDerivativeOrders(0), 0);
checkIndices(c.getPartialDerivativeOrders(1), 1);
c = DSCompiler.getCompiler(1, 2);
checkIndices(c.getPartialDerivativeOrders(0), 0);
checkIndices(c.getPartialDerivativeOrders(1), 1);
checkIndices(c.getPartialDerivativeOrders(2), 2);
c = DSCompiler.getCompiler(2, 1);
checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
checkIndices(c.getPartialDerivativeOrders(2), 0, 1);
c = DSCompiler.getCompiler(1, 3);
checkIndices(c.getPartialDerivativeOrders(0), 0);
checkIndices(c.getPartialDerivativeOrders(1), 1);
checkIndices(c.getPartialDerivativeOrders(2), 2);
checkIndices(c.getPartialDerivativeOrders(3), 3);
c = DSCompiler.getCompiler(2, 2);
checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
checkIndices(c.getPartialDerivativeOrders(3), 0, 1);
checkIndices(c.getPartialDerivativeOrders(4), 1, 1);
checkIndices(c.getPartialDerivativeOrders(5), 0, 2);
c = DSCompiler.getCompiler(3, 1);
checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0);
checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1);
c = DSCompiler.getCompiler(1, 4);
checkIndices(c.getPartialDerivativeOrders(0), 0);
checkIndices(c.getPartialDerivativeOrders(1), 1);
checkIndices(c.getPartialDerivativeOrders(2), 2);
checkIndices(c.getPartialDerivativeOrders(3), 3);
checkIndices(c.getPartialDerivativeOrders(4), 4);
c = DSCompiler.getCompiler(2, 3);
checkIndices(c.getPartialDerivativeOrders(0), 0, 0);
checkIndices(c.getPartialDerivativeOrders(1), 1, 0);
checkIndices(c.getPartialDerivativeOrders(2), 2, 0);
checkIndices(c.getPartialDerivativeOrders(3), 3, 0);
checkIndices(c.getPartialDerivativeOrders(4), 0, 1);
checkIndices(c.getPartialDerivativeOrders(5), 1, 1);
checkIndices(c.getPartialDerivativeOrders(6), 2, 1);
checkIndices(c.getPartialDerivativeOrders(7), 0, 2);
checkIndices(c.getPartialDerivativeOrders(8), 1, 2);
checkIndices(c.getPartialDerivativeOrders(9), 0, 3);
c = DSCompiler.getCompiler(3, 2);
checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0);
checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0);
checkIndices(c.getPartialDerivativeOrders(2), 2, 0, 0);
checkIndices(c.getPartialDerivativeOrders(3), 0, 1, 0);
checkIndices(c.getPartialDerivativeOrders(4), 1, 1, 0);
checkIndices(c.getPartialDerivativeOrders(5), 0, 2, 0);
checkIndices(c.getPartialDerivativeOrders(6), 0, 0, 1);
checkIndices(c.getPartialDerivativeOrders(7), 1, 0, 1);
checkIndices(c.getPartialDerivativeOrders(8), 0, 1, 1);
checkIndices(c.getPartialDerivativeOrders(9), 0, 0, 2);
c = DSCompiler.getCompiler(4, 1);
checkIndices(c.getPartialDerivativeOrders(0), 0, 0, 0, 0);
checkIndices(c.getPartialDerivativeOrders(1), 1, 0, 0, 0);
checkIndices(c.getPartialDerivativeOrders(2), 0, 1, 0, 0);
checkIndices(c.getPartialDerivativeOrders(3), 0, 0, 1, 0);
checkIndices(c.getPartialDerivativeOrders(4), 0, 0, 0, 1);
}
@Test
public void testSymmetry() {
for (int i = 0; i < 6; ++i) {
for (int j = 0; j < 6; ++j) {
DSCompiler c = DSCompiler.getCompiler(i, j);
for (int k = 0; k < c.getSize(); ++k) {
Assert.assertEquals(k, c.getPartialDerivativeIndex(c.getPartialDerivativeOrders(k)));
}
}
}
}
@Test public void testMultiplicationRules()
throws SecurityException, NoSuchFieldException, IllegalArgumentException, IllegalAccessException {
Map<String,String> referenceRules = new HashMap<String, String>();
referenceRules.put("(f*g)", "f * g");
referenceRules.put("d(f*g)/dx", "f * dg/dx + df/dx * g");
referenceRules.put("d(f*g)/dy", referenceRules.get("d(f*g)/dx").replaceAll("x", "y"));
referenceRules.put("d(f*g)/dz", referenceRules.get("d(f*g)/dx").replaceAll("x", "z"));
referenceRules.put("d(f*g)/dt", referenceRules.get("d(f*g)/dx").replaceAll("x", "t"));
referenceRules.put("d2(f*g)/dx2", "f * d2g/dx2 + 2 * df/dx * dg/dx + d2f/dx2 * g");
referenceRules.put("d2(f*g)/dy2", referenceRules.get("d2(f*g)/dx2").replaceAll("x", "y"));
referenceRules.put("d2(f*g)/dz2", referenceRules.get("d2(f*g)/dx2").replaceAll("x", "z"));
referenceRules.put("d2(f*g)/dt2", referenceRules.get("d2(f*g)/dx2").replaceAll("x", "t"));
referenceRules.put("d2(f*g)/dxdy", "f * d2g/dxdy + df/dy * dg/dx + df/dx * dg/dy + d2f/dxdy * g");
referenceRules.put("d2(f*g)/dxdz", referenceRules.get("d2(f*g)/dxdy").replaceAll("y", "z"));
referenceRules.put("d2(f*g)/dxdt", referenceRules.get("d2(f*g)/dxdy").replaceAll("y", "t"));
referenceRules.put("d2(f*g)/dydz", referenceRules.get("d2(f*g)/dxdz").replaceAll("x", "y"));
referenceRules.put("d2(f*g)/dydt", referenceRules.get("d2(f*g)/dxdt").replaceAll("x", "y"));
referenceRules.put("d2(f*g)/dzdt", referenceRules.get("d2(f*g)/dxdt").replaceAll("x", "z"));
referenceRules.put("d3(f*g)/dx3", "f * d3g/dx3 +" +
" 3 * df/dx * d2g/dx2 +" +
" 3 * d2f/dx2 * dg/dx +" +
" d3f/dx3 * g");
referenceRules.put("d3(f*g)/dy3", referenceRules.get("d3(f*g)/dx3").replaceAll("x", "y"));
referenceRules.put("d3(f*g)/dz3", referenceRules.get("d3(f*g)/dx3").replaceAll("x", "z"));
referenceRules.put("d3(f*g)/dt3", referenceRules.get("d3(f*g)/dx3").replaceAll("x", "t"));
referenceRules.put("d3(f*g)/dx2dy", "f * d3g/dx2dy +" +
" df/dy * d2g/dx2 +" +
" 2 * df/dx * d2g/dxdy +" +
" 2 * d2f/dxdy * dg/dx +" +
" d2f/dx2 * dg/dy +" +
" d3f/dx2dy * g");
referenceRules.put("d3(f*g)/dxdy2", "f * d3g/dxdy2 +" +
" 2 * df/dy * d2g/dxdy +" +
" d2f/dy2 * dg/dx +" +
" df/dx * d2g/dy2 +" +
" 2 * d2f/dxdy * dg/dy +" +
" d3f/dxdy2 * g");
referenceRules.put("d3(f*g)/dx2dz", referenceRules.get("d3(f*g)/dx2dy").replaceAll("y", "z"));
referenceRules.put("d3(f*g)/dy2dz", referenceRules.get("d3(f*g)/dx2dz").replaceAll("x", "y"));
referenceRules.put("d3(f*g)/dxdz2", referenceRules.get("d3(f*g)/dxdy2").replaceAll("y", "z"));
referenceRules.put("d3(f*g)/dydz2", referenceRules.get("d3(f*g)/dxdz2").replaceAll("x", "y"));
referenceRules.put("d3(f*g)/dx2dt", referenceRules.get("d3(f*g)/dx2dz").replaceAll("z", "t"));
referenceRules.put("d3(f*g)/dy2dt", referenceRules.get("d3(f*g)/dx2dt").replaceAll("x", "y"));
referenceRules.put("d3(f*g)/dz2dt", referenceRules.get("d3(f*g)/dx2dt").replaceAll("x", "z"));
referenceRules.put("d3(f*g)/dxdt2", referenceRules.get("d3(f*g)/dxdy2").replaceAll("y", "t"));
referenceRules.put("d3(f*g)/dydt2", referenceRules.get("d3(f*g)/dxdt2").replaceAll("x", "y"));
referenceRules.put("d3(f*g)/dzdt2", referenceRules.get("d3(f*g)/dxdt2").replaceAll("x", "z"));
referenceRules.put("d3(f*g)/dxdydz", "f * d3g/dxdydz +" +
" df/dz * d2g/dxdy +" +
" df/dy * d2g/dxdz +" +
" d2f/dydz * dg/dx +" +
" df/dx * d2g/dydz +" +
" d2f/dxdz * dg/dy +" +
" d2f/dxdy * dg/dz +" +
" d3f/dxdydz * g");
referenceRules.put("d3(f*g)/dxdydt", referenceRules.get("d3(f*g)/dxdydz").replaceAll("z", "t"));
referenceRules.put("d3(f*g)/dxdzdt", referenceRules.get("d3(f*g)/dxdydt").replaceAll("y", "z"));
referenceRules.put("d3(f*g)/dydzdt", referenceRules.get("d3(f*g)/dxdzdt").replaceAll("x", "y"));
Field multFieldArrayField = DSCompiler.class.getDeclaredField("multIndirection");
multFieldArrayField.setAccessible(true);
for (int i = 0; i < 5; ++i) {
for (int j = 0; j < 4; ++j) {
DSCompiler compiler = DSCompiler.getCompiler(i, j);
int[][][] multIndirection = (int[][][]) multFieldArrayField.get(compiler);
for (int k = 0; k < multIndirection.length; ++k) {
String product = ordersToString(compiler.getPartialDerivativeOrders(k),
"(f*g)", "x", "y", "z", "t");
StringBuilder rule = new StringBuilder();
for (int[] term : multIndirection[k]) {
if (rule.length() > 0) {
rule.append(" + ");
}
if (term[0] > 1) {
rule.append(term[0]).append(" * ");
}
rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[1]),
"f", "x", "y", "z", "t"));
rule.append(" * ");
rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[2]),
"g", "x", "y", "z", "t"));
}
Assert.assertEquals(product, referenceRules.get(product), rule.toString());
}
}
}
}
@Test public void testCompositionRules()
throws SecurityException, NoSuchFieldException, IllegalArgumentException, IllegalAccessException {
// the following reference rules have all been computed independently from the library,
// using only pencil and paper and some search and replace to handle symmetries
Map<String,String> referenceRules = new HashMap<String, String>();
referenceRules.put("(f(g))", "(f(g))");
referenceRules.put("d(f(g))/dx", "d(f(g))/dg * dg/dx");
referenceRules.put("d(f(g))/dy", referenceRules.get("d(f(g))/dx").replaceAll("x", "y"));
referenceRules.put("d(f(g))/dz", referenceRules.get("d(f(g))/dx").replaceAll("x", "z"));
referenceRules.put("d(f(g))/dt", referenceRules.get("d(f(g))/dx").replaceAll("x", "t"));
referenceRules.put("d2(f(g))/dx2", "d2(f(g))/dg2 * dg/dx * dg/dx + d(f(g))/dg * d2g/dx2");
referenceRules.put("d2(f(g))/dy2", referenceRules.get("d2(f(g))/dx2").replaceAll("x", "y"));
referenceRules.put("d2(f(g))/dz2", referenceRules.get("d2(f(g))/dx2").replaceAll("x", "z"));
referenceRules.put("d2(f(g))/dt2", referenceRules.get("d2(f(g))/dx2").replaceAll("x", "t"));
referenceRules.put("d2(f(g))/dxdy", "d2(f(g))/dg2 * dg/dx * dg/dy + d(f(g))/dg * d2g/dxdy");
referenceRules.put("d2(f(g))/dxdz", referenceRules.get("d2(f(g))/dxdy").replaceAll("y", "z"));
referenceRules.put("d2(f(g))/dxdt", referenceRules.get("d2(f(g))/dxdy").replaceAll("y", "t"));
referenceRules.put("d2(f(g))/dydz", referenceRules.get("d2(f(g))/dxdz").replaceAll("x", "y"));
referenceRules.put("d2(f(g))/dydt", referenceRules.get("d2(f(g))/dxdt").replaceAll("x", "y"));
referenceRules.put("d2(f(g))/dzdt", referenceRules.get("d2(f(g))/dxdt").replaceAll("x", "z"));
referenceRules.put("d3(f(g))/dx3", "d3(f(g))/dg3 * dg/dx * dg/dx * dg/dx +" +
" 3 * d2(f(g))/dg2 * dg/dx * d2g/dx2 +" +
" d(f(g))/dg * d3g/dx3");
referenceRules.put("d3(f(g))/dy3", referenceRules.get("d3(f(g))/dx3").replaceAll("x", "y"));
referenceRules.put("d3(f(g))/dz3", referenceRules.get("d3(f(g))/dx3").replaceAll("x", "z"));
referenceRules.put("d3(f(g))/dt3", referenceRules.get("d3(f(g))/dx3").replaceAll("x", "t"));
referenceRules.put("d3(f(g))/dxdy2", "d3(f(g))/dg3 * dg/dx * dg/dy * dg/dy +" +
" 2 * d2(f(g))/dg2 * dg/dy * d2g/dxdy +" +
" d2(f(g))/dg2 * dg/dx * d2g/dy2 +" +
" d(f(g))/dg * d3g/dxdy2");
referenceRules.put("d3(f(g))/dxdz2", referenceRules.get("d3(f(g))/dxdy2").replaceAll("y", "z"));
referenceRules.put("d3(f(g))/dxdt2", referenceRules.get("d3(f(g))/dxdy2").replaceAll("y", "t"));
referenceRules.put("d3(f(g))/dydz2", referenceRules.get("d3(f(g))/dxdz2").replaceAll("x", "y"));
referenceRules.put("d3(f(g))/dydt2", referenceRules.get("d3(f(g))/dxdt2").replaceAll("x", "y"));
referenceRules.put("d3(f(g))/dzdt2", referenceRules.get("d3(f(g))/dxdt2").replaceAll("x", "z"));
referenceRules.put("d3(f(g))/dx2dy", "d3(f(g))/dg3 * dg/dx * dg/dx * dg/dy +" +
" 2 * d2(f(g))/dg2 * dg/dx * d2g/dxdy +" +
" d2(f(g))/dg2 * d2g/dx2 * dg/dy +" +
" d(f(g))/dg * d3g/dx2dy");
referenceRules.put("d3(f(g))/dx2dz", referenceRules.get("d3(f(g))/dx2dy").replaceAll("y", "z"));
referenceRules.put("d3(f(g))/dx2dt", referenceRules.get("d3(f(g))/dx2dy").replaceAll("y", "t"));
referenceRules.put("d3(f(g))/dy2dz", referenceRules.get("d3(f(g))/dx2dz").replaceAll("x", "y"));
referenceRules.put("d3(f(g))/dy2dt", referenceRules.get("d3(f(g))/dx2dt").replaceAll("x", "y"));
referenceRules.put("d3(f(g))/dz2dt", referenceRules.get("d3(f(g))/dx2dt").replaceAll("x", "z"));
referenceRules.put("d3(f(g))/dxdydz", "d3(f(g))/dg3 * dg/dx * dg/dy * dg/dz +" +
" d2(f(g))/dg2 * dg/dy * d2g/dxdz +" +
" d2(f(g))/dg2 * dg/dx * d2g/dydz +" +
" d2(f(g))/dg2 * d2g/dxdy * dg/dz +" +
" d(f(g))/dg * d3g/dxdydz");
referenceRules.put("d3(f(g))/dxdydt", referenceRules.get("d3(f(g))/dxdydz").replaceAll("z", "t"));
referenceRules.put("d3(f(g))/dxdzdt", referenceRules.get("d3(f(g))/dxdydt").replaceAll("y", "z"));
referenceRules.put("d3(f(g))/dydzdt", referenceRules.get("d3(f(g))/dxdzdt").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dx4", "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dx * dg/dx +" +
" 6 * d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dx2 +" +
" 3 * d2(f(g))/dg2 * d2g/dx2 * d2g/dx2 +" +
" 4 * d2(f(g))/dg2 * dg/dx * d3g/dx3 +" +
" d(f(g))/dg * d4g/dx4");
referenceRules.put("d4(f(g))/dy4", referenceRules.get("d4(f(g))/dx4").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dz4", referenceRules.get("d4(f(g))/dx4").replaceAll("x", "z"));
referenceRules.put("d4(f(g))/dt4", referenceRules.get("d4(f(g))/dx4").replaceAll("x", "t"));
referenceRules.put("d4(f(g))/dx3dy", "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dx * dg/dy +" +
" 3 * d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dxdy +" +
" 3 * d3(f(g))/dg3 * dg/dx * d2g/dx2 * dg/dy +" +
" 3 * d2(f(g))/dg2 * d2g/dx2 * d2g/dxdy +" +
" 3 * d2(f(g))/dg2 * dg/dx * d3g/dx2dy +" +
" d2(f(g))/dg2 * d3g/dx3 * dg/dy +" +
" d(f(g))/dg * d4g/dx3dy");
referenceRules.put("d4(f(g))/dx3dz", referenceRules.get("d4(f(g))/dx3dy").replaceAll("y", "z"));
referenceRules.put("d4(f(g))/dx3dt", referenceRules.get("d4(f(g))/dx3dy").replaceAll("y", "t"));
referenceRules.put("d4(f(g))/dxdy3", "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dy * dg/dy +" +
" 3 * d3(f(g))/dg3 * dg/dy * dg/dy * d2g/dxdy +" +
" 3 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dy2 +" +
" 3 * d2(f(g))/dg2 * d2g/dxdy * d2g/dy2 +" +
" 3 * d2(f(g))/dg2 * dg/dy * d3g/dxdy2 +" +
" d2(f(g))/dg2 * dg/dx * d3g/dy3 +" +
" d(f(g))/dg * d4g/dxdy3");
referenceRules.put("d4(f(g))/dxdz3", referenceRules.get("d4(f(g))/dxdy3").replaceAll("y", "z"));
referenceRules.put("d4(f(g))/dxdt3", referenceRules.get("d4(f(g))/dxdy3").replaceAll("y", "t"));
referenceRules.put("d4(f(g))/dy3dz", referenceRules.get("d4(f(g))/dx3dz").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dy3dt", referenceRules.get("d4(f(g))/dx3dt").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dydz3", referenceRules.get("d4(f(g))/dxdz3").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dydt3", referenceRules.get("d4(f(g))/dxdt3").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dz3dt", referenceRules.get("d4(f(g))/dx3dt").replaceAll("x", "z"));
referenceRules.put("d4(f(g))/dzdt3", referenceRules.get("d4(f(g))/dxdt3").replaceAll("x", "z"));
referenceRules.put("d4(f(g))/dx2dy2", "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dy * dg/dy +" +
" 4 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dxdy +" +
" d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dy2 +" +
" 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dxdy +" +
" 2 * d2(f(g))/dg2 * dg/dx * d3g/dxdy2 +" +
" d3(f(g))/dg3 * d2g/dx2 * dg/dy * dg/dy +" +
" 2 * d2(f(g))/dg2 * dg/dy * d3g/dx2dy +" +
" d2(f(g))/dg2 * d2g/dx2 * d2g/dy2 +" +
" d(f(g))/dg * d4g/dx2dy2");
referenceRules.put("d4(f(g))/dx2dz2", referenceRules.get("d4(f(g))/dx2dy2").replaceAll("y", "z"));
referenceRules.put("d4(f(g))/dx2dt2", referenceRules.get("d4(f(g))/dx2dy2").replaceAll("y", "t"));
referenceRules.put("d4(f(g))/dy2dz2", referenceRules.get("d4(f(g))/dx2dz2").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dy2dt2", referenceRules.get("d4(f(g))/dx2dt2").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dz2dt2", referenceRules.get("d4(f(g))/dx2dt2").replaceAll("x", "z"));
referenceRules.put("d4(f(g))/dx2dydz", "d4(f(g))/dg4 * dg/dx * dg/dx * dg/dy * dg/dz +" +
" 2 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dxdz +" +
" d3(f(g))/dg3 * dg/dx * dg/dx * d2g/dydz +" +
" 2 * d3(f(g))/dg3 * dg/dx * d2g/dxdy * dg/dz +" +
" 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dxdz +" +
" 2 * d2(f(g))/dg2 * dg/dx * d3g/dxdydz +" +
" d3(f(g))/dg3 * d2g/dx2 * dg/dy * dg/dz +" +
" d2(f(g))/dg2 * dg/dy * d3g/dx2dz +" +
" d2(f(g))/dg2 * d2g/dx2 * d2g/dydz +" +
" d2(f(g))/dg2 * d3g/dx2dy * dg/dz +" +
" d(f(g))/dg * d4g/dx2dydz");
referenceRules.put("d4(f(g))/dx2dydt", referenceRules.get("d4(f(g))/dx2dydz").replaceAll("z", "t"));
referenceRules.put("d4(f(g))/dx2dzdt", referenceRules.get("d4(f(g))/dx2dydt").replaceAll("y", "z"));
referenceRules.put("d4(f(g))/dxdy2dz", "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dy * dg/dz +" +
" d3(f(g))/dg3 * dg/dy * dg/dy * d2g/dxdz +" +
" 2 * d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dydz +" +
" 2 * d3(f(g))/dg3 * dg/dy * d2g/dxdy * dg/dz +" +
" 2 * d2(f(g))/dg2 * d2g/dxdy * d2g/dydz +" +
" 2 * d2(f(g))/dg2 * dg/dy * d3g/dxdydz +" +
" d3(f(g))/dg3 * dg/dx * d2g/dy2 * dg/dz +" +
" d2(f(g))/dg2 * d2g/dy2 * d2g/dxdz +" +
" d2(f(g))/dg2 * dg/dx * d3g/dy2dz +" +
" d2(f(g))/dg2 * d3g/dxdy2 * dg/dz +" +
" d(f(g))/dg * d4g/dxdy2dz");
referenceRules.put("d4(f(g))/dxdy2dt", referenceRules.get("d4(f(g))/dxdy2dz").replaceAll("z", "t"));
referenceRules.put("d4(f(g))/dy2dzdt", referenceRules.get("d4(f(g))/dx2dzdt").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dxdydz2", "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dz * dg/dz +" +
" 2 * d3(f(g))/dg3 * dg/dy * dg/dz * d2g/dxdz +" +
" 2 * d3(f(g))/dg3 * dg/dx * dg/dz * d2g/dydz +" +
" d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dz2 +" +
" 2 * d2(f(g))/dg2 * d2g/dxdz * d2g/dydz +" +
" d2(f(g))/dg2 * dg/dy * d3g/dxdz2 +" +
" d2(f(g))/dg2 * dg/dx * d3g/dydz2 +" +
" d3(f(g))/dg3 * d2g/dxdy * dg/dz * dg/dz +" +
" 2 * d2(f(g))/dg2 * dg/dz * d3g/dxdydz +" +
" d2(f(g))/dg2 * d2g/dxdy * d2g/dz2 +" +
" d(f(g))/dg * d4g/dxdydz2");
referenceRules.put("d4(f(g))/dxdz2dt", referenceRules.get("d4(f(g))/dxdy2dt").replaceAll("y", "z"));
referenceRules.put("d4(f(g))/dydz2dt", referenceRules.get("d4(f(g))/dxdz2dt").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dxdydt2", referenceRules.get("d4(f(g))/dxdydz2").replaceAll("z", "t"));
referenceRules.put("d4(f(g))/dxdzdt2", referenceRules.get("d4(f(g))/dxdydt2").replaceAll("y", "z"));
referenceRules.put("d4(f(g))/dydzdt2", referenceRules.get("d4(f(g))/dxdzdt2").replaceAll("x", "y"));
referenceRules.put("d4(f(g))/dxdydzdt", "d4(f(g))/dg4 * dg/dx * dg/dy * dg/dz * dg/dt +" +
" d3(f(g))/dg3 * dg/dy * dg/dz * d2g/dxdt +" +
" d3(f(g))/dg3 * dg/dx * dg/dz * d2g/dydt +" +
" d3(f(g))/dg3 * dg/dx * dg/dy * d2g/dzdt +" +
" d3(f(g))/dg3 * dg/dy * d2g/dxdz * dg/dt +" +
" d2(f(g))/dg2 * d2g/dxdz * d2g/dydt +" +
" d2(f(g))/dg2 * dg/dy * d3g/dxdzdt +" +
" d3(f(g))/dg3 * dg/dx * d2g/dydz * dg/dt +" +
" d2(f(g))/dg2 * d2g/dydz * d2g/dxdt +" +
" d2(f(g))/dg2 * dg/dx * d3g/dydzdt +" +
" d3(f(g))/dg3 * d2g/dxdy * dg/dz * dg/dt +" +
" d2(f(g))/dg2 * dg/dz * d3g/dxdydt +" +
" d2(f(g))/dg2 * d2g/dxdy * d2g/dzdt +" +
" d2(f(g))/dg2 * d3g/dxdydz * dg/dt +" +
" d(f(g))/dg * d4g/dxdydzdt");
Field compFieldArrayField = DSCompiler.class.getDeclaredField("compIndirection");
compFieldArrayField.setAccessible(true);
for (int i = 0; i < 5; ++i) {
for (int j = 0; j < 5; ++j) {
DSCompiler compiler = DSCompiler.getCompiler(i, j);
int[][][] compIndirection = (int[][][]) compFieldArrayField.get(compiler);
for (int k = 0; k < compIndirection.length; ++k) {
String product = ordersToString(compiler.getPartialDerivativeOrders(k),
"(f(g))", "x", "y", "z", "t");
StringBuilder rule = new StringBuilder();
for (int[] term : compIndirection[k]) {
if (rule.length() > 0) {
rule.append(" + ");
}
if (term[0] > 1) {
rule.append(term[0]).append(" * ");
}
rule.append(orderToString(term[1], "(f(g))", "g"));
for (int l = 2; l < term.length; ++l) {
rule.append(" * ");
rule.append(ordersToString(compiler.getPartialDerivativeOrders(term[l]),
"g", "x", "y", "z", "t"));
}
}
Assert.assertEquals(product, referenceRules.get(product), rule.toString());
}
}
}
}
private void checkIndices(int[] indices, int ... expected) {
Assert.assertEquals(expected.length, indices.length);
for (int i = 0; i < expected.length; ++i) {
Assert.assertEquals(expected[i], indices[i]);
}
}
private String orderToString(int order, String functionName, String parameterName) {
if (order == 0) {
return functionName;
} else if (order == 1) {
return "d" + functionName + "/d" + parameterName;
} else {
return "d" + order + functionName + "/d" + parameterName + order;
}
}
private String ordersToString(int[] orders, String functionName, String ... parametersNames) {
int sumOrders = 0;
for (int order : orders) {
sumOrders += order;
}
if (sumOrders == 0) {
return functionName;
}
StringBuilder builder = new StringBuilder();
builder.append('d');
if (sumOrders > 1) {
builder.append(sumOrders);
}
builder.append(functionName).append('/');
for (int i = 0; i < orders.length; ++i) {
if (orders[i] > 0) {
builder.append('d').append(parametersNames[i]);
if (orders[i] > 1) {
builder.append(orders[i]);
}
}
}
return builder.toString();
}
}

539
src/test/java/org/apache/commons/math3/analysis/differentiation/DerivativeStructureTest.java Normal file
View File

@ -0,0 +1,539 @@
/*
* Licensed to the Apache Software Foundation (ASF) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* The ASF licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.apache.commons.math3.analysis.differentiation;
import java.util.Arrays;
import java.util.List;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.util.ArithmeticUtils;
import org.apache.commons.math3.util.FastMath;
import org.junit.Assert;
import org.junit.Test;
/**
* Test for class {@link DerivativeStructure}.
*/
public class DerivativeStructureTest {
@Test(expected=NumberIsTooLargeException.class)
public void testWrongVariableIndex() {
new DerivativeStructure(3, 1, 3, 1.0);
}
@Test(expected=DimensionMismatchException.class)
public void testMissingOrders() {
new DerivativeStructure(3, 1, 0, 1.0).getPartialDerivative(0, 1);
}
@Test(expected=NumberIsTooLargeException.class)
public void testTooLargeOrder() {
new DerivativeStructure(3, 1, 0, 1.0).getPartialDerivative(1, 1, 2);
}
@Test
public void testVariableWithoutDerivative0() {
DerivativeStructure v = new DerivativeStructure(1, 0, 0, 1.0);
Assert.assertEquals(1.0, v.getValue(), 1.0e-15);
}
@Test(expected=NumberIsTooLargeException.class)
public void testVariableWithoutDerivative1() {
DerivativeStructure v = new DerivativeStructure(1, 0, 0, 1.0);
Assert.assertEquals(1.0, v.getPartialDerivative(1), 1.0e-15);
}
@Test
public void testVariable() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0),
1.0, 1.0, 0.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0),
2.0, 0.0, 1.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0),
3.0, 0.0, 0.0, 1.0);
}
}
@Test
public void testConstant() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
checkF0F1(new DerivativeStructure(3, maxOrder, FastMath.PI),
FastMath.PI, 0.0, 0.0, 0.0);
}
}
@Test
public void testPrimitiveAdd() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).add(5), 6.0, 1.0, 0.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).add(5), 7.0, 0.0, 1.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).add(5), 8.0, 0.0, 0.0, 1.0);
}
}
@Test
public void testAdd() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
DerivativeStructure xyz = x.add(y.add(z));
checkF0F1(xyz, x.getValue() + y.getValue() + z.getValue(), 1.0, 1.0, 1.0);
}
}
@Test
public void testPrimitiveSubtract() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).subtract(5), -4.0, 1.0, 0.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).subtract(5), -3.0, 0.0, 1.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).subtract(5), -2.0, 0.0, 0.0, 1.0);
}
}
@Test
public void testSubtract() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
DerivativeStructure xyz = x.subtract(y.subtract(z));
checkF0F1(xyz, x.getValue() - (y.getValue() - z.getValue()), 1.0, -1.0, 1.0);
}
}
@Test
public void testPrimitiveMultiply() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).multiply(5), 5.0, 5.0, 0.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).multiply(5), 10.0, 0.0, 5.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).multiply(5), 15.0, 0.0, 0.0, 5.0);
}
}
@Test
public void testMultiply() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
DerivativeStructure xyz = x.multiply(y.multiply(z));
for (int i = 0; i <= maxOrder; ++i) {
for (int j = 0; j <= maxOrder; ++j) {
for (int k = 0; k <= maxOrder; ++k) {
if (i + j + k <= maxOrder) {
Assert.assertEquals((i == 0 ? x.getValue() : (i == 1 ? 1.0 : 0.0)) *
(j == 0 ? y.getValue() : (j == 1 ? 1.0 : 0.0)) *
(k == 0 ? z.getValue() : (k == 1 ? 1.0 : 0.0)),
xyz.getPartialDerivative(i, j, k),
1.0e-15);
}
}
}
}
}
}
@Test
public void testNegate() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
checkF0F1(new DerivativeStructure(3, maxOrder, 0, 1.0).negate(), -1.0, -1.0, 0.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 1, 2.0).negate(), -2.0, 0.0, -1.0, 0.0);
checkF0F1(new DerivativeStructure(3, maxOrder, 2, 3.0).negate(), -3.0, 0.0, 0.0, -1.0);
}
}
@Test
public void testReciprocal() {
for (double x = 0.1; x < 1.2; x += 0.1) {
DerivativeStructure r = new DerivativeStructure(1, 6, 0, x).reciprocal();
Assert.assertEquals(1 / x, r.getValue(), 1.0e-15);
for (int i = 1; i < r.getOrder(); ++i) {
double expected = ArithmeticUtils.pow(-1, i) * ArithmeticUtils.factorial(i) /
FastMath.pow(x, i + 1);
Assert.assertEquals(expected, r.getPartialDerivative(i), 1.0e-15 * FastMath.abs(expected));
}
}
}
@Test
public void testPower() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
for (int n = 0; n < 10; ++n) {
DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
DerivativeStructure y = new DerivativeStructure(3, maxOrder, 1, 2.0);
DerivativeStructure z = new DerivativeStructure(3, maxOrder, 2, 3.0);
List<DerivativeStructure> list = Arrays.asList(x, y, z,
x.add(y).add(z),
x.multiply(y).multiply(z));
if (n == 0) {
for (DerivativeStructure ds : list) {
checkEquals(ds.getField().getOne(), ds.pow(n), 1.0e-15);
}
} else if (n == 1) {
for (DerivativeStructure ds : list) {
checkEquals(ds, ds.pow(n), 1.0e-15);
}
} else {
for (DerivativeStructure ds : list) {
DerivativeStructure p = ds.getField().getOne();
for (int i = 0; i < n; ++i) {
p = p.multiply(ds);
}
checkEquals(p, ds.pow(n), 1.0e-15);
}
}
}
}
}
@Test
public void testExpression() {
double epsilon = 2.5e-13;
for (double x = 0; x < 2; x += 0.2) {
DerivativeStructure dsX = new DerivativeStructure(3, 5, 0, x);
for (double y = 0; y < 2; y += 0.2) {
DerivativeStructure dsY = new DerivativeStructure(3, 5, 1, y);
for (double z = 0; z >- 2; z -= 0.2) {
DerivativeStructure dsZ = new DerivativeStructure(3, 5, 2, z);
// f(x, y, z) = x + 5 x y - 2 z + (8 z x - y)^3
DerivativeStructure ds =
new DerivativeStructure(1, dsX,
5, dsX.multiply(dsY),
-2, dsZ,
1, new DerivativeStructure(8, dsZ.multiply(dsX),
-1, dsY).pow(3));
double f = x + 5 * x * y - 2 * z + FastMath.pow(8 * z * x - y, 3);
Assert.assertEquals(f, ds.getValue(),
FastMath.abs(epsilon * f));
// df/dx = 1 + 5 y + 24 (8 z x - y)^2 z
double dfdx = 1 + 5 * y + 24 * z * FastMath.pow(8 * z * x - y, 2);
Assert.assertEquals(dfdx, ds.getPartialDerivative(1, 0, 0),
FastMath.abs(epsilon * dfdx));
// df/dxdy = 5 + 48 z*(y - 8 z x)
double dfdxdy = 5 + 48 * z * (y - 8 * z * x);
Assert.assertEquals(dfdxdy, ds.getPartialDerivative(1, 1, 0),
FastMath.abs(epsilon * dfdxdy));
// df/dxdydz = 48 (y - 16 z x)
double dfdxdydz = 48 * (y - 16 * z * x);
Assert.assertEquals(dfdxdydz, ds.getPartialDerivative(1, 1, 1),
FastMath.abs(epsilon * dfdxdydz));
}
}
}
}
@Test
public void testCompositionOneVariableX() {
double epsilon = 1.0e-13;
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.1) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
for (double y = 0.1; y < 1.2; y += 0.1) {
DerivativeStructure dsY = new DerivativeStructure(1, maxOrder, y);
DerivativeStructure f = dsX.divide(dsY).sqrt();
double f0 = FastMath.sqrt(x / y);
Assert.assertEquals(f0, f.getValue(), FastMath.abs(epsilon * f0));
if (f.getOrder() > 0) {
double f1 = 1 / (2 * FastMath.sqrt(x * y));
Assert.assertEquals(f1, f.getPartialDerivative(1), FastMath.abs(epsilon * f1));
if (f.getOrder() > 1) {
double f2 = -f1 / (2 * x);
Assert.assertEquals(f2, f.getPartialDerivative(2), FastMath.abs(epsilon * f2));
if (f.getOrder() > 2) {
double f3 = (f0 + x / (2 * y * f0)) / (4 * x * x * x);
Assert.assertEquals(f3, f.getPartialDerivative(3), FastMath.abs(epsilon * f3));
}
}
}
}
}
}
}
@Test
public void testTrigo() {
double epsilon = 2.0e-12;
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.1) {
DerivativeStructure dsX = new DerivativeStructure(3, maxOrder, 0, x);
for (double y = 0.1; y < 1.2; y += 0.1) {
DerivativeStructure dsY = new DerivativeStructure(3, maxOrder, 1, y);
for (double z = 0.1; z < 1.2; z += 0.1) {
DerivativeStructure dsZ = new DerivativeStructure(3, maxOrder, 2, z);
DerivativeStructure f = dsX.divide(dsY.cos().add(dsZ.tan())).sin();
double a = FastMath.cos(y) + FastMath.tan(z);
double f0 = FastMath.sin(x / a);
Assert.assertEquals(f0, f.getValue(), FastMath.abs(epsilon * f0));
if (f.getOrder() > 0) {
double dfdx = FastMath.cos(x / a) / a;
Assert.assertEquals(dfdx, f.getPartialDerivative(1, 0, 0), FastMath.abs(epsilon * dfdx));
double dfdy = x * FastMath.sin(y) * dfdx / a;
Assert.assertEquals(dfdy, f.getPartialDerivative(0, 1, 0), FastMath.abs(epsilon * dfdy));
double cz = FastMath.cos(z);
double cz2 = cz * cz;
double dfdz = -x * dfdx / (a * cz2);
Assert.assertEquals(dfdz, f.getPartialDerivative(0, 0, 1), FastMath.abs(epsilon * dfdz));
if (f.getOrder() > 1) {
double df2dx2 = -(f0 / (a * a));
Assert.assertEquals(df2dx2, f.getPartialDerivative(2, 0, 0), FastMath.abs(epsilon * df2dx2));
double df2dy2 = x * FastMath.cos(y) * dfdx / a -
x * x * FastMath.sin(y) * FastMath.sin(y) * f0 / (a * a * a * a) +
2 * FastMath.sin(y) * dfdy / a;
Assert.assertEquals(df2dy2, f.getPartialDerivative(0, 2, 0), FastMath.abs(epsilon * df2dy2));
double c4 = cz2 * cz2;
double df2dz2 = x * (2 * a * (1 - a * cz * FastMath.sin(z)) * dfdx - x * f0 / a ) / (a * a * a * c4);
Assert.assertEquals(df2dz2, f.getPartialDerivative(0, 0, 2), FastMath.abs(epsilon * df2dz2));
double df2dxdy = dfdy / x - x * FastMath.sin(y) * f0 / (a * a * a);
Assert.assertEquals(df2dxdy, f.getPartialDerivative(1, 1, 0), FastMath.abs(epsilon * df2dxdy));
}
}
}
}
}
}
}
@Test
public void testSqrtDefinition() {
double[] epsilon = new double[] { 5.0e-16, 5.0e-16, 2.0e-15, 5.0e-14, 2.0e-12 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure sqrt1 = dsX.pow(0.5);
DerivativeStructure sqrt2 = dsX.sqrt();
DerivativeStructure zero = sqrt1.subtract(sqrt2);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testSqrtPow2() {
double[] epsilon = new double[] { 1.0e-16, 3.0e-16, 2.0e-15, 6.0e-14, 6.0e-12 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure rebuiltX = dsX.multiply(dsX).sqrt();
DerivativeStructure zero = rebuiltX.subtract(dsX);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testCbrtDefinition() {
double[] epsilon = new double[] { 4.0e-16, 9.0e-16, 6.0e-15, 2.0e-13, 4.0e-12 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure cbrt1 = dsX.pow(1.0 / 3.0);
DerivativeStructure cbrt2 = dsX.cbrt();
DerivativeStructure zero = cbrt1.subtract(cbrt2);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testCbrtPow3() {
double[] epsilon = new double[] { 1.0e-16, 5.0e-16, 8.0e-15, 3.0e-13, 4.0e-11 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure rebuiltX = dsX.multiply(dsX.multiply(dsX)).cbrt();
DerivativeStructure zero = rebuiltX.subtract(dsX);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testExp() {
double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16, 1.0e-16 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
double refExp = FastMath.exp(x);
DerivativeStructure exp = new DerivativeStructure(1, maxOrder, 0, x).exp();
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(refExp, exp.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testLog() {
double[] epsilon = new double[] { 1.0e-16, 1.0e-16, 3.0e-14, 7.0e-13, 3.0e-11 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure log = new DerivativeStructure(1, maxOrder, 0, x).log();
Assert.assertEquals(FastMath.log(x), log.getValue(), epsilon[0]);
for (int n = 1; n <= maxOrder; ++n) {
double refDer = -ArithmeticUtils.factorial(n - 1) / FastMath.pow(-x, n);
Assert.assertEquals(refDer, log.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testLogExp() {
double[] epsilon = new double[] { 2.0e-16, 2.0e-16, 3.0e-16, 2.0e-15, 6.0e-15 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure rebuiltX = dsX.exp().log();
DerivativeStructure zero = rebuiltX.subtract(dsX);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0.0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testTangentDefinition() {
double[] epsilon = new double[] { 5.0e-16, 2.0e-15, 3.0e-14, 5.0e-13, 2.0e-11 };
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.001) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, 0, x);
DerivativeStructure tan1 = dsX.sin().divide(dsX.cos());
DerivativeStructure tan2 = dsX.tan();
DerivativeStructure zero = tan1.subtract(tan2);
for (int n = 0; n <= maxOrder; ++n) {
Assert.assertEquals(0, zero.getPartialDerivative(n), epsilon[n]);
}
}
}
}
@Test
public void testCompositionOneVariableY() {
double epsilon = 1.0e-13;
for (int maxOrder = 0; maxOrder < 5; ++maxOrder) {
for (double x = 0.1; x < 1.2; x += 0.1) {
DerivativeStructure dsX = new DerivativeStructure(1, maxOrder, x);
for (double y = 0.1; y < 1.2; y += 0.1) {
DerivativeStructure dsY = new DerivativeStructure(1, maxOrder, 0, y);
DerivativeStructure f = dsX.divide(dsY).sqrt();
double f0 = FastMath.sqrt(x / y);
Assert.assertEquals(f0, f.getValue(), FastMath.abs(epsilon * f0));
if (f.getOrder() > 0) {
double f1 = -x / (2 * y * y * f0);
Assert.assertEquals(f1, f.getPartialDerivative(1), FastMath.abs(epsilon * f1));
if (f.getOrder() > 1) {
double f2 = (f0 - x / (4 * y * f0)) / (y * y);
Assert.assertEquals(f2, f.getPartialDerivative(2), FastMath.abs(epsilon * f2));
if (f.getOrder() > 2) {
double f3 = (x / (8 * y * f0) - 2 * f0) / (y * y * y);
Assert.assertEquals(f3, f.getPartialDerivative(3), FastMath.abs(epsilon * f3));
}
}
}
}
}
}
}
@Test
public void testField() {
for (int maxOrder = 1; maxOrder < 5; ++maxOrder) {
DerivativeStructure x = new DerivativeStructure(3, maxOrder, 0, 1.0);
checkF0F1(x.getField().getZero(), 0.0, 0.0, 0.0, 0.0);
checkF0F1(x.getField().getOne(), 1.0, 0.0, 0.0, 0.0);
Assert.assertEquals(maxOrder, x.getField().getZero().getOrder());
Assert.assertEquals(3, x.getField().getZero().getFreeParameters());
Assert.assertEquals(DerivativeStructure.class, x.getField().getRuntimeClass());
}
}
private void checkF0F1(DerivativeStructure ds, double value, double...derivatives) {
// check dimension
Assert.assertEquals(derivatives.length, ds.getFreeParameters());
// check value, directly and also as 0th order derivative
Assert.assertEquals(value, ds.getValue(), 1.0e-15);
Assert.assertEquals(value, ds.getPartialDerivative(new int[ds.getFreeParameters()]), 1.0e-15);
// check first order derivatives
for (int i = 0; i < derivatives.length; ++i) {
int[] orders = new int[derivatives.length];
orders[i] = 1;
Assert.assertEquals(derivatives[i], ds.getPartialDerivative(orders), 1.0e-15);
}
}
private void checkEquals(DerivativeStructure ds1, DerivativeStructure ds2, double epsilon) {
// check dimension
Assert.assertEquals(ds1.getFreeParameters(), ds2.getFreeParameters());
Assert.assertEquals(ds1.getOrder(), ds2.getOrder());
int[] derivatives = new int[ds1.getFreeParameters()];
int sum = 0;
while (true) {
if (sum <= ds1.getOrder()) {
Assert.assertEquals(ds1.getPartialDerivative(derivatives),
ds2.getPartialDerivative(derivatives),
epsilon);
}
boolean increment = true;
sum = 0;
for (int i = derivatives.length - 1; i >= 0; --i) {
if (increment) {
if (derivatives[i] == ds1.getOrder()) {
derivatives[i] = 0;
} else {
derivatives[i]++;
increment = false;
}
}
sum += derivatives[i];
}
if (increment) {
return;
}
}
}
}