/*
* 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.analysis.UnivariateFunction;
import org.apache.commons.math3.analysis.UnivariateMatrixFunction;
import org.apache.commons.math3.analysis.UnivariateVectorFunction;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.NotPositiveException;
import org.apache.commons.math3.exception.NumberIsTooLargeException;
import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.util.FastMath;
/** Univariate functions differentiator using finite differences.
* <p>
* This class creates some wrapper objects around regular
* {@link UnivariateFunction univariate functions} (or {@link
* UnivariateVectorFunction univariate vector functions} or {@link
* UnivariateMatrixFunction univariate matrix functions}). These
* wrapper objects compute derivatives in addition to function
* values.
* </p>
* <p>
* The wrapper objects work by calling the underlying function on
* a sampling grid around the current point and performing polynomial
* interpolation. A finite differences scheme with n points is
* theoretically able to compute derivatives up to order n-1, but
* it is generally better to have a slight margin. The step size must
* also be small enough in order for the polynomial approximation to
* be good in the current point neighborhood, but it should not be too
* small because numerical instability appears quickly (there are several
* differences of close points). Choosing the number of points and
* the step size is highly problem dependent.
* </p>
* <p>
* As an example of good and bad settings, lets consider the quintic
* polynomial function {@code f(x) = (x-1)*(x-0.5)*x*(x+0.5)*(x+1)}.
* Since it is a polynomial, finite differences with at least 6 points
* should theoretically recover the exact same polynomial and hence
* compute accurate derivatives for any order. However, due to numerical
* errors, we get the following results for a 7 points finite differences
* for abscissae in the [-10, 10] range:
* <ul>
* <li>step size = 0.25, second order derivative error about 9.97e-10</li>
* <li>step size = 0.25, fourth order derivative error about 5.43e-8</li>
* <li>step size = 1.0e-6, second order derivative error about 148</li>
* <li>step size = 1.0e-6, fourth order derivative error about 6.35e+14</li>
* </ul>
* <p>
* This example shows that the small step size is really bad, even simply
* for second order derivative!</p>
*
* @since 3.1
*/
public class FiniteDifferencesDifferentiator
implements UnivariateFunctionDifferentiator, UnivariateVectorFunctionDifferentiator,
UnivariateMatrixFunctionDifferentiator, Serializable {
/** Serializable UID. */
private static final long serialVersionUID = 20120917L;
/** Number of points to use. */
private final int nbPoints;
/** Step size. */
private final double stepSize;
/** Half sample span. */
private final double halfSampleSpan;
/** Lower bound for independent variable. */
private final double tMin;
/** Upper bound for independent variable. */
private final double tMax;
/**
* Build a differentiator with number of points and step size when independent variable is unbounded.
* <p>
* Beware that wrong settings for the finite differences differentiator
* can lead to highly unstable and inaccurate results, especially for
* high derivation orders. Using very small step sizes is often a
* <em>bad</em> idea.
* </p>
* @param nbPoints number of points to use
* @param stepSize step size (gap between each point)
* @exception NotPositiveException if {@code stepsize <= 0} (note that
* {@link NotPositiveException} extends {@link NumberIsTooSmallException})
* @exception NumberIsTooSmallException {@code nbPoint <= 1}
*/
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize)
throws NotPositiveException, NumberIsTooSmallException {
this(nbPoints, stepSize, Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY);
}
/**
* Build a differentiator with number of points and step size when independent variable is bounded.
* <p>
* When the independent variable is bounded (tLower < t < tUpper), the sampling
* points used for differentiation will be adapted to ensure the constraint holds
* even near the boundaries. This means the sample will not be centered anymore in
* these cases. At an extreme case, computing derivatives exactly at the lower bound
* will lead the sample to be entirely on the right side of the derivation point.
* </p>
* <p>
* Note that the boundaries are considered to be excluded for function evaluation.
* </p>
* <p>
* Beware that wrong settings for the finite differences differentiator
* can lead to highly unstable and inaccurate results, especially for
* high derivation orders. Using very small step sizes is often a
* <em>bad</em> idea.
* </p>
* @param nbPoints number of points to use
* @param stepSize step size (gap between each point)
* @param tLower lower bound for independent variable (may be {@code Double.NEGATIVE_INFINITY}
* if there are no lower bounds)
* @param tUpper upper bound for independent variable (may be {@code Double.POSITIVE_INFINITY}
* if there are no upper bounds)
* @exception NotPositiveException if {@code stepsize <= 0} (note that
* {@link NotPositiveException} extends {@link NumberIsTooSmallException})
* @exception NumberIsTooSmallException {@code nbPoint <= 1}
* @exception NumberIsTooLargeException {@code stepSize * (nbPoints - 1) >= tUpper - tLower}
*/
public FiniteDifferencesDifferentiator(final int nbPoints, final double stepSize,
final double tLower, final double tUpper)
throws NotPositiveException, NumberIsTooSmallException, NumberIsTooLargeException {
if (nbPoints <= 1) {
throw new NumberIsTooSmallException(stepSize, 1, false);
}
this.nbPoints = nbPoints;
if (stepSize <= 0) {
throw new NotPositiveException(stepSize);
}
this.stepSize = stepSize;
halfSampleSpan = 0.5 * stepSize * (nbPoints - 1);
if (2 * halfSampleSpan >= tUpper - tLower) {
throw new NumberIsTooLargeException(2 * halfSampleSpan, tUpper - tLower, false);
}
final double safety = FastMath.ulp(halfSampleSpan);
this.tMin = tLower + halfSampleSpan + safety;
this.tMax = tUpper - halfSampleSpan - safety;
}
/**
* Get the number of points to use.
* @return number of points to use
*/
public int getNbPoints() {
return nbPoints;
}
/**
* Get the step size.
* @return step size
*/
public double getStepSize() {
return stepSize;
}
/**
* Evaluate derivatives from a sample.
* <p>
* Evaluation is done using divided differences.
* </p>
* @param t evaluation abscissa value and derivatives
* @param t0 first sample point abscissa
* @param y function values sample {@code y[i] = f(t[i]) = f(t0 + i * stepSize)}
* @return value and derivatives at {@code t}
* @exception NumberIsTooLargeException if the requested derivation order
* is larger or equal to the number of points
*/
private DerivativeStructure evaluate(final DerivativeStructure t, final double t0,
final double[] y)
throws NumberIsTooLargeException {
// create divided differences diagonal arrays
final double[] top = new double[nbPoints];
final double[] bottom = new double[nbPoints];
for (int i = 0; i < nbPoints; ++i) {
// update the bottom diagonal of the divided differences array
bottom[i] = y[i];
for (int j = 1; j <= i; ++j) {
bottom[i - j] = (bottom[i - j + 1] - bottom[i - j]) / (j * stepSize);
}
// update the top diagonal of the divided differences array
top[i] = bottom[0];
}
// evaluate interpolation polynomial (represented by top diagonal) at t
final int order = t.getOrder();
final int parameters = t.getFreeParameters();
final double[] derivatives = t.getAllDerivatives();
final double dt0 = t.getValue() - t0;
DerivativeStructure interpolation = new DerivativeStructure(parameters, order, 0.0);
DerivativeStructure monomial = null;
for (int i = 0; i < nbPoints; ++i) {
if (i == 0) {
// start with monomial(t) = 1
monomial = new DerivativeStructure(parameters, order, 1.0);
} else {
// monomial(t) = (t - t0) * (t - t1) * ... * (t - t(i-1))
derivatives[0] = dt0 - (i - 1) * stepSize;
final DerivativeStructure deltaX = new DerivativeStructure(parameters, order, derivatives);
monomial = monomial.multiply(deltaX);
}
interpolation = interpolation.add(monomial.multiply(top[i]));
}
return interpolation;
}
/** {@inheritDoc}
* <p>The returned object cannot compute derivatives to arbitrary orders. The
* value function will throw a {@link NumberIsTooLargeException} if the requested
* derivation order is larger or equal to the number of points.
* </p>
*/
public UnivariateDifferentiableFunction differentiate(final UnivariateFunction function) {
return new UnivariateDifferentiableFunction() {
/** {@inheritDoc} */
public double value(final double x) throws MathIllegalArgumentException {
return function.value(x);
}
/** {@inheritDoc} */
public DerivativeStructure value(final DerivativeStructure t)
throws MathIllegalArgumentException {
// check we can achieve the requested derivation order with the sample
if (t.getOrder() >= nbPoints) {
throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
}
// compute sample position, trying to be centered if possible
final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
// compute sample points
final double[] y = new double[nbPoints];
for (int i = 0; i < nbPoints; ++i) {
y[i] = function.value(t0 + i * stepSize);
}
// evaluate derivatives
return evaluate(t, t0, y);
}
};
}
/** {@inheritDoc}
* <p>The returned object cannot compute derivatives to arbitrary orders. The
* value function will throw a {@link NumberIsTooLargeException} if the requested
* derivation order is larger or equal to the number of points.
* </p>
*/
public UnivariateDifferentiableVectorFunction differentiate(final UnivariateVectorFunction function) {
return new UnivariateDifferentiableVectorFunction() {
/** {@inheritDoc} */
public double[]value(final double x) throws MathIllegalArgumentException {
return function.value(x);
}
/** {@inheritDoc} */
public DerivativeStructure[] value(final DerivativeStructure t)
throws MathIllegalArgumentException {
// check we can achieve the requested derivation order with the sample
if (t.getOrder() >= nbPoints) {
throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
}
// compute sample position, trying to be centered if possible
final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
// compute sample points
double[][] y = null;
for (int i = 0; i < nbPoints; ++i) {
final double[] v = function.value(t0 + i * stepSize);
if (i == 0) {
y = new double[v.length][nbPoints];
}
for (int j = 0; j < v.length; ++j) {
y[j][i] = v[j];
}
}
// evaluate derivatives
final DerivativeStructure[] value = new DerivativeStructure[y.length];
for (int j = 0; j < value.length; ++j) {
value[j] = evaluate(t, t0, y[j]);
}
return value;
}
};
}
/** {@inheritDoc}
* <p>The returned object cannot compute derivatives to arbitrary orders. The
* value function will throw a {@link NumberIsTooLargeException} if the requested
* derivation order is larger or equal to the number of points.
* </p>
*/
public UnivariateDifferentiableMatrixFunction differentiate(final UnivariateMatrixFunction function) {
return new UnivariateDifferentiableMatrixFunction() {
/** {@inheritDoc} */
public double[][] value(final double x) throws MathIllegalArgumentException {
return function.value(x);
}
/** {@inheritDoc} */
public DerivativeStructure[][] value(final DerivativeStructure t)
throws MathIllegalArgumentException {
// check we can achieve the requested derivation order with the sample
if (t.getOrder() >= nbPoints) {
throw new NumberIsTooLargeException(t.getOrder(), nbPoints, false);
}
// compute sample position, trying to be centered if possible
final double t0 = FastMath.max(FastMath.min(t.getValue(), tMax), tMin) - halfSampleSpan;
// compute sample points
double[][][] y = null;
for (int i = 0; i < nbPoints; ++i) {
final double[][] v = function.value(t0 + i * stepSize);
if (i == 0) {
y = new double[v.length][v[0].length][nbPoints];
}
for (int j = 0; j < v.length; ++j) {
for (int k = 0; k < v[j].length; ++k) {
y[j][k][i] = v[j][k];
}
}
}
// evaluate derivatives
final DerivativeStructure[][] value = new DerivativeStructure[y.length][y[0].length];
for (int j = 0; j < value.length; ++j) {
for (int k = 0; k < y[j].length; ++k) {
value[j][k] = evaluate(t, t0, y[j][k]);
}
}
return value;
}
};
}
}