/**
* Copyright (C) 2012 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.analytics.math.differentiation;
import org.apache.commons.lang.Validate;
import com.opengamma.analytics.math.MathException;
import com.opengamma.analytics.math.function.Function1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.analytics.math.matrix.MatrixAlgebra;
import com.opengamma.analytics.math.matrix.OGMatrixAlgebra;
import com.opengamma.util.ArgumentChecker;
/**
*
*/
public class MatrixFieldFirstOrderDifferentiator implements Differentiator<DoubleMatrix1D, DoubleMatrix2D, DoubleMatrix2D[]> {
private static final MatrixAlgebra MA = new OGMatrixAlgebra();
private static final double DEFAULT_EPS = 1e-5;
private final double _eps;
private final double _twoEps;
private final double _oneOverTwpEps;
public MatrixFieldFirstOrderDifferentiator() {
_eps = DEFAULT_EPS;
_twoEps = 2 * DEFAULT_EPS;
_oneOverTwpEps = 1.0 / _twoEps;
}
public MatrixFieldFirstOrderDifferentiator(final double eps) {
ArgumentChecker.isTrue(eps > 1e-15, "eps of {} is below machine tolerance of 1e-15. Please choose a higher value or use default", eps);
_eps = eps;
_twoEps = 2 * eps;
_oneOverTwpEps = 1.0 / _twoEps;
}
@Override
public Function1D<DoubleMatrix1D, DoubleMatrix2D[]> differentiate(final Function1D<DoubleMatrix1D, DoubleMatrix2D> function) {
Validate.notNull(function);
return new Function1D<DoubleMatrix1D, DoubleMatrix2D[]>() {
@SuppressWarnings("synthetic-access")
@Override
public DoubleMatrix2D[] evaluate(final DoubleMatrix1D x) {
Validate.notNull(x, "x");
final int n = x.getNumberOfElements();
final DoubleMatrix2D[] res = new DoubleMatrix2D[n];
final double[] xData = x.getData();
for (int i = 0; i < n; i++) {
final double oldValue = xData[i];
xData[i] += _eps;
final DoubleMatrix2D up = function.evaluate(x);
xData[i] -= _twoEps;
final DoubleMatrix2D down = function.evaluate(x);
res[i] = (DoubleMatrix2D) MA.scale(MA.subtract(up, down), _oneOverTwpEps); //TODO have this in one operation
xData[i] = oldValue;
}
return res;
}
};
}
@Override
public Function1D<DoubleMatrix1D, DoubleMatrix2D[]> differentiate(final Function1D<DoubleMatrix1D, DoubleMatrix2D> function, final Function1D<DoubleMatrix1D, Boolean> domain) {
Validate.notNull(function);
Validate.notNull(domain);
final double[] wFwd = new double[] {-3. / _twoEps, 4. / _twoEps, -1. / _twoEps };
final double[] wCent = new double[] {-1. / _twoEps, 0., 1. / _twoEps };
final double[] wBack = new double[] {1. / _twoEps, -4. / _twoEps, 3. / _twoEps };
return new Function1D<DoubleMatrix1D, DoubleMatrix2D[]>() {
@SuppressWarnings("synthetic-access")
@Override
public DoubleMatrix2D[] evaluate(final DoubleMatrix1D x) {
Validate.notNull(x, "x");
ArgumentChecker.isTrue(domain.evaluate(x), "point {} is not in the function domain", x.toString());
final int n = x.getNumberOfElements();
final double[] xData = x.getData();
double oldValue;
final DoubleMatrix2D[] y = new DoubleMatrix2D[3];
final DoubleMatrix2D[] res = new DoubleMatrix2D[n];
double[] w;
for (int i = 0; i < n; i++) {
oldValue = xData[i];
xData[i] += _eps;
if (!domain.evaluate(x)) {
xData[i] = oldValue - _twoEps;
if (!domain.evaluate(x)) {
throw new MathException("cannot get derivative at point " + x.toString() + " in direction " + i);
}
y[0] = function.evaluate(x);
xData[i] = oldValue;
y[2] = function.evaluate(x);
xData[i] = oldValue - _eps;
y[1] = function.evaluate(x);
w = wBack;
} else {
final DoubleMatrix2D temp = function.evaluate(x);
xData[i] = oldValue - _eps;
if (!domain.evaluate(x)) {
y[1] = temp;
xData[i] = oldValue;
y[0] = function.evaluate(x);
xData[i] = oldValue + _twoEps;
y[2] = function.evaluate(x);
w = wFwd;
} else {
y[2] = temp;
xData[i] = oldValue - _eps;
y[0] = function.evaluate(x);
w = wCent;
}
}
res[i] = (DoubleMatrix2D) MA.add(MA.scale(y[0], w[0]), MA.scale(y[2], w[2]));
if (w[1] != 0) {
res[i] = (DoubleMatrix2D) MA.add(res[i], MA.scale(y[1], w[1]));
}
xData[i] = oldValue;
}
return res;
}
};
}
}