PiecewisePolynomialWithSensitivityFunction1D.java example

Explorer
OG-Platform-master
- projects
- sesame
/**
 * 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);
  }
}