com.opengamma.strata.math.impl.interpolation.CubicSplineNakSolver Maven / Gradle / Ivy

/*
 * Copyright (C) 2013 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.strata.math.impl.interpolation;

import java.util.Arrays;

import com.opengamma.strata.collect.array.DoubleArray;
import com.opengamma.strata.collect.array.DoubleMatrix;

/**
 * Solves cubic spline problem with Not-A-Knot endpoint conditions, where the third derivative
 * at the endpoints is the same as that of their adjacent points.
 */
public class CubicSplineNakSolver extends CubicSplineSolver {

  @Override
  public DoubleMatrix solve(final double[] xValues, final double[] yValues) {

    final double[] intervals = getDiffs(xValues);

    return getSplineCoeffs(xValues, yValues, intervals, matrixEqnSolver(getMatrix(intervals), getVector(yValues, intervals)));
  }

  @Override
  public DoubleMatrix[] solveWithSensitivity(final double[] xValues, final double[] yValues) {
    final double[] intervals = getDiffs(xValues);
    final double[][] toBeInv = getMatrix(intervals);
    final double[] vector = getVector(yValues, intervals);
    final double[][] vecSensitivity = getVectorSensitivity(intervals);

    return getSplineCoeffsWithSensitivity(xValues, yValues, intervals, toBeInv, vector, vecSensitivity);
  }

  @Override
  public DoubleMatrix[] solveMultiDim(final double[] xValues, final DoubleMatrix yValuesMatrix) {
    final int dim = yValuesMatrix.rowCount();
    final DoubleMatrix[] coefMatrix = new DoubleMatrix[dim];

    for (int i = 0; i < dim; ++i) {
      coefMatrix[i] = solve(xValues, yValuesMatrix.row(i).toArray());
    }

    return coefMatrix;
  }

  @Override
  public DoubleArray getKnotsMat1D(final double[] xValues) {
    final int nData = xValues.length;
    if (nData == 2) {
      return DoubleArray.of(xValues[0], xValues[nData - 1]);
    }
    if (nData == 3) {
      return DoubleArray.of(xValues[0], xValues[nData - 1]);
    }
    return DoubleArray.copyOf(xValues);
  }

  /**
   * @param xValues X values of Data
   * @param yValues Y values of Data
   * @param intervals {xValues[1]-xValues[0], xValues[2]-xValues[1],...}
   * @param solnVector Values of second derivative at knots
   * @return Coefficient matrix whose i-th row vector is (a_0,a_1,...) for i-th intervals, where a_0,a_1,... are coefficients of f(x) = a_0 + a_1 x^1 + ....
   */
  private DoubleMatrix getSplineCoeffs(final double[] xValues, final double[] yValues, final double[] intervals, final double[] solnVector) {
    final int nData = xValues.length;

    if (nData == 2) {
      final double[][] res = new double[][] {{
        yValues[1] / intervals[0] - yValues[0] / intervals[0] - intervals[0] * solnVector[0] / 2. - intervals[0] * solnVector[1] / 6. + intervals[0] * solnVector[0] / 6., yValues[0] } };
      return DoubleMatrix.copyOf(res);
    }
    if (nData == 3) {
      final double[][] res = new double[][] {{solnVector[0] / 2., yValues[1] / intervals[0] - yValues[0] / intervals[0] - intervals[0] * solnVector[0] / 2., yValues[0] } };
      return DoubleMatrix.copyOf(res);
    }
    return getCommonSplineCoeffs(xValues, yValues, intervals, solnVector);
  }

  private DoubleMatrix[] getSplineCoeffsWithSensitivity(final double[] xValues, final double[] yValues, final double[] intervals, final double[][] toBeInv, final double[] vector,
      final double[][] vecSensitivity) {
    final int nData = xValues.length;

    if (nData == 2) {
      final DoubleMatrix[] res = new DoubleMatrix[nData];
      res[0] = DoubleMatrix.of(1, 1, yValues[1] / intervals[0] - yValues[0] / intervals[0], yValues[0]);
      res[1] = DoubleMatrix.of(2, 2, -1d / intervals[0], 1d / intervals[0], 1d, 0d);
      return res;
    }
    if (nData == 3) {
      final DoubleMatrix[] res = new DoubleMatrix[2];
      final DoubleArray[] soln = combinedMatrixEqnSolver(toBeInv, vector, vecSensitivity);
      final double[][] coef = new double[][] {{soln[0].get(0) / 2.,
          yValues[1] / intervals[0] - yValues[0] / intervals[0] - intervals[0] * soln[0].get(0) / 2., yValues[0]}};
      res[0] = DoubleMatrix.copyOf(coef);
      final double[][] coefSense = new double[3][0];
      coefSense[0] = new double[] {soln[1].get(0) / 2., soln[2].get(0) / 2., soln[3].get(0) / 2.};
      coefSense[1] = new double[]
          {-1. / intervals[0] - intervals[0] * soln[1].get(0) / 2., 1. / intervals[0] - intervals[0] * soln[2].get(0) / 2.,
              -intervals[0] * soln[3].get(0) / 2.};
      coefSense[2] = new double[] {1., 0., 0. };
      res[1] = DoubleMatrix.copyOf(coefSense);
      return res;
    }
    final DoubleMatrix[] res = new DoubleMatrix[nData];
    final DoubleArray[] soln = combinedMatrixEqnSolver(toBeInv, vector, vecSensitivity);
    res[0] = getCommonSplineCoeffs(xValues, yValues, intervals, soln[0].toArray());
    final double[][] solnMatrix = new double[nData][nData];
    for (int i = 0; i < nData; ++i) {
      for (int j = 0; j < nData; ++j) {
        solnMatrix[i][j] = soln[j + 1].get(i);
      }
    }
    final DoubleMatrix[] tmp = getCommonSensitivityCoeffs(intervals, solnMatrix);
    System.arraycopy(tmp, 0, res, 1, nData - 1);
    return res;
  }

  /**
   * Cubic spline is obtained by solving a linear problem Ax=b where A is a square matrix and x,b are vector
   * @param yValues Y Values of data
   * @param intervals {xValues[1]-xValues[0], xValues[2]-xValues[1],...}
   * @return Vector b
   */
  private double[] getVector(final double[] yValues, final double[] intervals) {

    final int nData = yValues.length;
    double[] res = new double[nData];

    if (nData == 3) {
      for (int i = 0; i < nData; ++i) {
        res[i] = 2. * yValues[2] / (intervals[0] + intervals[1]) - 2. * yValues[0] / (intervals[0] + intervals[1]) - 2. * yValues[1] / (intervals[0]) + 2. * yValues[0] / (intervals[0]);
      }
    } else {
      res = getCommonVectorElements(yValues, intervals);
    }
    return res;
  }

  private double[][] getVectorSensitivity(final double[] intervals) {

    final int nData = intervals.length + 1;
    double[][] res = new double[nData][nData];

    if (nData == 3) {
      for (int i = 0; i < nData; ++i) {
        res[i][0] = -2. / (intervals[0] + intervals[1]) + 2. / (intervals[0]);
        res[i][1] = -2. / (intervals[0]);
        res[i][2] = 2. / (intervals[0] + intervals[1]);
      }
    } else {
      res = getCommonVectorSensitivity(intervals);
    }
    return res;
  }

  /**
   * Cubic spline is obtained by solving a linear problem Ax=b where A is a square matrix and x,b are vector
   * @param intervals {xValues[1]-xValues[0], xValues[2]-xValues[1],...}
   * @return Matrix A
   */
  private double[][] getMatrix(final double[] intervals) {

    final int nData = intervals.length + 1;
    double[][] res = new double[nData][nData];

    for (int i = 0; i < nData; ++i) {
      Arrays.fill(res[i], 0.);
    }

    if (nData == 2) {
      res[0][1] = intervals[0];
      res[1][0] = intervals[0];
      return res;
    }
    if (nData == 3) {
      res[0][0] = intervals[1];
      res[1][1] = intervals[1];
      res[2][2] = intervals[1];
      return res;
    }
    res = getCommonMatrixElements(intervals);
    res[0][0] = -intervals[1];
    res[0][1] = intervals[0] + intervals[1];
    res[0][2] = -intervals[0];
    res[nData - 1][nData - 3] = -intervals[nData - 2];
    res[nData - 1][nData - 2] = intervals[nData - 3] + intervals[nData - 2];
    res[nData - 1][nData - 1] = -intervals[nData - 3];
    return res;
  }

}