/**
* Copyright (C) 2013 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.analytics.math.function;
import org.apache.commons.lang.NotImplementedException;
import com.opengamma.analytics.math.FunctionUtils;
import com.opengamma.analytics.math.interpolation.PiecewisePolynomialResultsWithSensitivity;
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;
/**
* Give a class {@link PiecewisePolynomialResultsWithSensitivity}, Compute node sensitivity of function value, first derivative value and second derivative value
*/
public class PiecewisePolynomialWithSensitivityFunction1D extends PiecewisePolynomialFunction1D {
private static final MatrixAlgebra MA = new OGMatrixAlgebra();
/**
* @param pp {@link PiecewisePolynomialResultsWithSensitivity}
* @param xKey
* @return Node sensitivity value at x=xKey
*/
public DoubleMatrix1D nodeSensitivity(final PiecewisePolynomialResultsWithSensitivity pp, final double xKey) {
ArgumentChecker.notNull(pp, "null pp");
ArgumentChecker.isFalse(Double.isNaN(xKey), "xKey containing NaN");
ArgumentChecker.isFalse(Double.isInfinite(xKey), "xKey containing Infinity");
if (pp.getDimensions() > 1) {
throw new NotImplementedException();
}
final double[] knots = pp.getKnots().getData();
final int nKnots = knots.length;
int interval = FunctionUtils.getLowerBoundIndex(knots, xKey);
if (interval == nKnots - 1) {
interval--; // there is 1 less interval that knots
}
final double s = xKey - knots[interval];
final DoubleMatrix2D a = pp.getCoefficientSensitivity(interval);
final int nCoefs = a.getNumberOfRows();
DoubleMatrix1D res = a.getRowVector(0);
for (int i = 1; i < nCoefs; i++) {
res = (DoubleMatrix1D) MA.scale(res, s);
res = (DoubleMatrix1D) MA.add(res, a.getRowVector(i));
}
return res;
}
/**
* @param pp {@link PiecewisePolynomialResultsWithSensitivity}
* @param xKeys
* @return Node sensitivity value at x=xKeys
*/
public DoubleMatrix1D[] nodeSensitivity(final PiecewisePolynomialResultsWithSensitivity pp, final double[] xKeys) {
ArgumentChecker.notNull(pp, "null pp");
ArgumentChecker.notNull(xKeys, "null xKeys");
final int nKeys = xKeys.length;
final DoubleMatrix1D[] res = new DoubleMatrix1D[nKeys];
for (int i = 0; i < nKeys; ++i) {
ArgumentChecker.isFalse(Double.isNaN(xKeys[i]), "xKey containing NaN");
ArgumentChecker.isFalse(Double.isInfinite(xKeys[i]), "xKey containing Infinity");
}
if (pp.getDimensions() > 1) {
throw new NotImplementedException();
}
final double[] knots = pp.getKnots().getData();
final int nKnots = knots.length;
for (int j = 0; j < nKeys; ++j) {
final double xKey = xKeys[j];
int interval = FunctionUtils.getLowerBoundIndex(knots, xKey);
if (interval == nKnots - 1) {
interval--; // there is 1 less interval that knots
}
final double s = xKey - knots[interval];
final DoubleMatrix2D a = pp.getCoefficientSensitivity(interval);
final int nCoefs = a.getNumberOfRows();
res[j] = a.getRowVector(0);
for (int i = 1; i < nCoefs; i++) {
res[j] = (DoubleMatrix1D) MA.scale(res[j], s);
res[j] = (DoubleMatrix1D) MA.add(res[j], a.getRowVector(i));
}
}
return res;
}
/**
* @param pp {@link PiecewisePolynomialResultsWithSensitivity}
* @param xKey
* @return Node sensitivity of derivative value at x=xKey
*/
public DoubleMatrix1D differentiateNodeSensitivity(final PiecewisePolynomialResultsWithSensitivity pp, final double xKey) {
ArgumentChecker.notNull(pp, "null pp");
ArgumentChecker.isFalse(Double.isNaN(xKey), "xKey containing NaN");
ArgumentChecker.isFalse(Double.isInfinite(xKey), "xKey containing Infinity");
if (pp.getDimensions() > 1) {
throw new NotImplementedException();
}
final int nCoefs = pp.getOrder();
ArgumentChecker.isFalse(nCoefs < 2, "Polynomial degree is too low");
final double[] knots = pp.getKnots().getData();
final int nKnots = knots.length;
int interval = FunctionUtils.getLowerBoundIndex(knots, xKey);
if (interval == nKnots - 1) {
interval--; // there is 1 less interval that knots
}
final double s = xKey - knots[interval];
final DoubleMatrix2D a = pp.getCoefficientSensitivity(interval);
DoubleMatrix1D res = (DoubleMatrix1D) MA.scale(a.getRowVector(0), nCoefs - 1);
for (int i = 1; i < nCoefs - 1; i++) {
res = (DoubleMatrix1D) MA.scale(res, s);
res = (DoubleMatrix1D) MA.add(res, MA.scale(a.getRowVector(i), nCoefs - 1 - i));
}
return res;
}
/**
* @param pp {@link PiecewisePolynomialResultsWithSensitivity}
* @param xKeys
* @return Node sensitivity of derivative value at x=xKeys
*/
public DoubleMatrix1D[] differentiateNodeSensitivity(final PiecewisePolynomialResultsWithSensitivity pp, final double[] xKeys) {
ArgumentChecker.notNull(pp, "null pp");
if (pp.getDimensions() > 1) {
throw new NotImplementedException();
}
final int nCoefs = pp.getOrder();
ArgumentChecker.isFalse(nCoefs < 2, "Polynomial degree is too low");
final int nIntervals = pp.getNumberOfIntervals();
final DoubleMatrix2D[] diffSense = new DoubleMatrix2D[nIntervals];
final DoubleMatrix2D[] senseMat = pp.getCoefficientSensitivityAll();
final int nData = senseMat[0].getNumberOfColumns();
for (int i = 0; i < nIntervals; ++i) {
final double[][] tmp = new double[nCoefs - 1][nData];
for (int j = 0; j < nCoefs - 1; ++j) {
for (int k = 0; k < nData; ++k) {
tmp[j][k] = (nCoefs - 1 - j) * senseMat[i].getData()[j][k];
}
}
diffSense[i] = new DoubleMatrix2D(tmp);
}
PiecewisePolynomialResultsWithSensitivity ppDiff = new PiecewisePolynomialResultsWithSensitivity(pp.getKnots(), pp.getCoefMatrix(), nCoefs - 1, pp.getDimensions(), diffSense);
return nodeSensitivity(ppDiff, xKeys);
}
/**
* @param pp {@link PiecewisePolynomialResultsWithSensitivity}
* @param xKey
* @return Node sensitivity of second derivative value at x=xKey
*/
public DoubleMatrix1D differentiateTwiceNodeSensitivity(final PiecewisePolynomialResultsWithSensitivity pp, final double xKey) {
ArgumentChecker.notNull(pp, "null pp");
ArgumentChecker.isFalse(Double.isNaN(xKey), "xKey containing NaN");
ArgumentChecker.isFalse(Double.isInfinite(xKey), "xKey containing Infinity");
if (pp.getDimensions() > 1) {
throw new NotImplementedException();
}
final int nCoefs = pp.getOrder();
ArgumentChecker.isFalse(nCoefs < 3, "Polynomial degree is too low");
final double[] knots = pp.getKnots().getData();
final int nKnots = knots.length;
int interval = FunctionUtils.getLowerBoundIndex(knots, xKey);
if (interval == nKnots - 1) {
interval--; // there is 1 less interval that knots
}
final double s = xKey - knots[interval];
final DoubleMatrix2D a = pp.getCoefficientSensitivity(interval);
DoubleMatrix1D res = (DoubleMatrix1D) MA.scale(a.getRowVector(0), (nCoefs - 1) * (nCoefs - 2));
for (int i = 1; i < nCoefs - 2; i++) {
res = (DoubleMatrix1D) MA.scale(res, s);
res = (DoubleMatrix1D) MA.add(res, MA.scale(a.getRowVector(i), (nCoefs - 1 - i) * (nCoefs - 2 - i)));
}
return res;
}
/**
* @param pp {@link PiecewisePolynomialResultsWithSensitivity}
* @param xKeys
* @return Node sensitivity of second derivative value at x=xKeys
*/
public DoubleMatrix1D[] differentiateTwiceNodeSensitivity(final PiecewisePolynomialResultsWithSensitivity pp, final double[] xKeys) {
ArgumentChecker.notNull(pp, "null pp");
if (pp.getDimensions() > 1) {
throw new NotImplementedException();
}
final int nCoefs = pp.getOrder();
ArgumentChecker.isFalse(nCoefs < 3, "Polynomial degree is too low");
final int nIntervals = pp.getNumberOfIntervals();
final DoubleMatrix2D[] diffSense = new DoubleMatrix2D[nIntervals];
final DoubleMatrix2D[] senseMat = pp.getCoefficientSensitivityAll();
final int nData = senseMat[0].getNumberOfColumns();
for (int i = 0; i < nIntervals; ++i) {
final double[][] tmp = new double[nCoefs - 2][nData];
for (int j = 0; j < nCoefs - 2; ++j) {
for (int k = 0; k < nData; ++k) {
tmp[j][k] = (nCoefs - 1 - j) * (nCoefs - 2 - j) * senseMat[i].getData()[j][k];
}
}
diffSense[i] = new DoubleMatrix2D(tmp);
}
PiecewisePolynomialResultsWithSensitivity ppDiff = new PiecewisePolynomialResultsWithSensitivity(pp.getKnots(), pp.getCoefMatrix(), nCoefs - 2, pp.getDimensions(), diffSense);
return nodeSensitivity(ppDiff, xKeys);
}
}