Java tutorial
/** * Copyright (C) 2013 - present by OpenGamma Inc. and the OpenGamma group of companies * * Please see distribution for license. */ package com.opengamma.analytics.math.interpolation; import java.util.Arrays; import org.apache.commons.lang.NotImplementedException; import com.google.common.primitives.Doubles; import com.opengamma.analytics.math.FunctionUtils; import com.opengamma.analytics.math.matrix.DoubleMatrix1D; import com.opengamma.analytics.math.matrix.DoubleMatrix2D; import com.opengamma.util.ArgumentChecker; import com.opengamma.util.ParallelArrayBinarySort; /** * C1 cubic interpolation preserving monotonicity based on * Fritsch, F. N.; Carlson, R. E. (1980) * "Monotone Piecewise Cubic Interpolation", SIAM Journal on Numerical Analysis 17 (2): 238246. * Fritsch, F. N. and Butland, J. (1984) * "A method for constructing local monotone piecewise cubic interpolants", SIAM Journal on Scientific and Statistical Computing 5 (2): 300-304. * * For interpolation without node sensitivity, use {@link PiecewiseCubicHermiteSplineInterpolator} */ public class PiecewiseCubicHermiteSplineInterpolatorWithSensitivity extends PiecewisePolynomialInterpolator { @Override public PiecewisePolynomialResultsWithSensitivity interpolateWithSensitivity(final double[] xValues, final double[] yValues) { ArgumentChecker.notNull(xValues, "xValues"); ArgumentChecker.notNull(yValues, "yValues"); ArgumentChecker.isTrue(xValues.length == yValues.length, "xValues length = yValues length"); ArgumentChecker.isTrue(xValues.length > 1, "Data points should be more than 1"); final int nDataPts = xValues.length; for (int i = 0; i < nDataPts; ++i) { ArgumentChecker.isFalse(Double.isNaN(xValues[i]), "xData containing NaN"); ArgumentChecker.isFalse(Double.isInfinite(xValues[i]), "xData containing Infinity"); ArgumentChecker.isFalse(Double.isNaN(yValues[i]), "yData containing NaN"); ArgumentChecker.isFalse(Double.isInfinite(yValues[i]), "yData containing Infinity"); } double[] xValuesSrt = Arrays.copyOf(xValues, nDataPts); double[] yValuesSrt = Arrays.copyOf(yValues, nDataPts); ParallelArrayBinarySort.parallelBinarySort(xValuesSrt, yValuesSrt); for (int i = 1; i < nDataPts; ++i) { ArgumentChecker.isFalse(xValuesSrt[i - 1] == xValuesSrt[i], "xValues should be distinct"); } final DoubleMatrix2D[] temp = solve(xValuesSrt, yValuesSrt); // check the matrices // TODO remove some of these tests ArgumentChecker.noNulls(temp, "error in solve - some matrices are null"); int n = temp.length; ArgumentChecker.isTrue(n == nDataPts, "wrong number of matricies"); for (int k = 0; k < n; k++) { DoubleMatrix2D m = temp[k]; final int rows = m.getNumberOfRows(); final int cols = m.getNumberOfColumns(); for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols; ++j) { ArgumentChecker.isTrue(Doubles.isFinite(m.getEntry(i, j)), "Matrix contains a NaN or infinite"); } } } DoubleMatrix2D coefMatrix = temp[0]; DoubleMatrix2D[] coefMatrixSense = new DoubleMatrix2D[n - 1]; System.arraycopy(temp, 1, coefMatrixSense, 0, n - 1); return new PiecewisePolynomialResultsWithSensitivity(new DoubleMatrix1D(xValuesSrt), coefMatrix, 4, 1, coefMatrixSense); } /** * @param xValues X values of data * @param yValues Y values of data * @return Coefficient matrix whose i-th row vector is {a3, a2, a1, a0} of f(x) = a3 * (x-x_i)^3 + a2 * (x-x_i)^2 +... for the i-th interval */ private DoubleMatrix2D[] solve(final double[] xValues, final double[] yValues) { final int n = xValues.length; double[][] coeff = new double[n - 1][4]; double[] h = new double[n - 1]; double[] delta = new double[n - 1]; DoubleMatrix2D[] res = new DoubleMatrix2D[n]; for (int i = 0; i < n - 1; ++i) { h[i] = xValues[i + 1] - xValues[i]; delta[i] = (yValues[i + 1] - yValues[i]) / h[i]; } if (n == 2) { // TODO check this - should be yValues coeff[0][2] = delta[0]; coeff[0][3] = xValues[0]; } else { SlopeFinderResults temp = slopeFinder(h, delta, yValues); final double[] d = temp.getSlopes().getData(); final double[][] dDy = temp.getSlopeJacobian().getData(); // form up the coefficient matrix for (int i = 0; i < n - 1; ++i) { coeff[i][0] = (d[i] - 2 * delta[i] + d[i + 1]) / h[i] / h[i]; // b coeff[i][1] = (3 * delta[i] - 2. * d[i] - d[i + 1]) / h[i]; // c coeff[i][2] = d[i]; coeff[i][3] = yValues[i]; } // // TODO this would all be a lot nicer if we had multiplication of sparse matrices // double[][] dDy = new double[n][]; // for (int i = 0; i < n; i++) { // final double[] vec = new double[n]; // vec[0] = -dDDelta[0][i] / h[0]; // vec[n - 1] = dDDelta[n - 2][i] / h[n - 2]; // for (int j = 1; j < n - 1; j++) { // vec[j] = -dDDelta[j][i] / h[j] + dDDelta[j - 1][i] / h[j - 1]; // } // dDy[i] = vec; // } double[][] bDy = new double[n - 1][n]; double[][] cDy = new double[n - 1][n]; for (int i = 0; i < n - 1; i++) { final double invH = 1 / h[i]; final double invH2 = invH * invH; final double invH3 = invH * invH2; cDy[i][i] = -3 * invH2; cDy[i][i + 1] = 3 * invH2; bDy[i][i] = 2 * invH3; bDy[i][i + 1] = -2 * invH3; for (int j = 0; j < n; j++) { cDy[i][j] -= (2 * dDy[i][j] + dDy[i + 1][j]) * invH; bDy[i][j] += (dDy[i][j] + dDy[i + 1][j]) * invH2; } } // Now we have to pack this into an array of DoubleMatrix2D - my kingdom for a tensor class res[0] = new DoubleMatrix2D(coeff); for (int k = 0; k < n - 1; k++) { double[][] coeffSense = new double[4][]; coeffSense[0] = bDy[k]; coeffSense[1] = cDy[k]; coeffSense[2] = dDy[k]; coeffSense[3] = new double[n]; coeffSense[3][k] = 1.0; res[k + 1] = new DoubleMatrix2D(coeffSense); } } return res; } private class SlopeFinderResults { private final DoubleMatrix1D _d; private final DoubleMatrix2D _dDy; public SlopeFinderResults(final DoubleMatrix1D d, final DoubleMatrix2D dDy) { // this is a private class - don't do the normal checks on inputs _d = d; _dDy = dDy; } public DoubleMatrix1D getSlopes() { return _d; } public DoubleMatrix2D getSlopeJacobian() { return _dDy; } } /** * Finds the the first derivatives at knots and their sensitivity to delta * @param h * @param delta * @return slope finder results */ private SlopeFinderResults slopeFinder(final double[] h, final double[] delta, final double[] y) { final int n = y.length; final double[] invDelta = new double[n - 1]; final double[] invDelta2 = new double[n - 1]; final double[] invH = new double[n - 1]; for (int i = 0; i < (n - 1); i++) { invDelta[i] = 1 / delta[i]; invDelta2[i] = invDelta[i] * invDelta[i]; invH[i] = 1 / h[i]; } final double[] d = new double[n]; // TODO it would be better if this were a sparse matrix final double[][] jac = new double[n][n]; // internal points for (int i = 1; i < n - 1; ++i) { if (delta[i] * delta[i - 1] > 0.) { final double w1 = 2. * h[i] + h[i - 1]; final double w2 = h[i] + 2. * h[i - 1]; final double w12 = w1 + w2; d[i] = w12 / (w1 * invDelta[i - 1] + w2 * invDelta[i]); final double z1 = d[i] * d[i] / w12; jac[i][i - 1] = -w1 * invH[i - 1] * invDelta2[i - 1] * z1; jac[i][i] = (w1 * invH[i - 1] * invDelta2[i - 1] - w2 * invH[i] * invDelta2[i]) * z1; jac[i][i + 1] = w2 * invH[i] * invDelta2[i] * z1; } else if (delta[i] == 0 ^ delta[i - 1] == 0) { // d is zero, so we don't explicitly set it final double w1 = 2. * h[i] + h[i - 1]; final double w2 = h[i] + 2. * h[i - 1]; final double w12 = w1 + w2; final double z2 = 0.5 * w12 / FunctionUtils.square(w1 * delta[i] + w2 * delta[i - 1]); jac[i][i - 1] = -w1 * invH[i - 1] * delta[i] * delta[i] * z2; jac[i][i] = (w1 * invH[i - 1] * delta[i] * delta[i] - w2 * invH[i] * delta[i - 1] * delta[i - 1]) * z2; jac[i][i + 1] = w2 * invH[i] * delta[i - 1] * delta[i - 1] * z2; } } // fill in end points double[] temp = endpointSlope(h[0], h[1], delta[0], delta[1], false); d[0] = temp[0]; for (int i = 0; i < 3; i++) { jac[0][i] = temp[i + 1]; } temp = endpointSlope(h[n - 2], h[n - 3], delta[n - 2], delta[n - 3], true); d[n - 1] = temp[0]; for (int i = 1; i < 4; i++) { jac[n - 1][n - i] = temp[i]; } return new SlopeFinderResults(new DoubleMatrix1D(d), new DoubleMatrix2D(jac)); } /** * First derivative at end point and its sensitivity to delta * @param h1 * @param h2 * @param y1 * @param y2 * @param y3 * @return array of length 4 - the first element contains d, while the other three are sensitivities to the ys */ private double[] endpointSlope(final double h1, final double h2, final double del1, final double del2, final boolean rightSide) { final double[] res = new double[4]; if (del1 == 0.0) { // quick exist for particular edge case // d and dDy3 are both zero - no need to explicitly set if (del2 == 0) { res[1] = -(2 * h1 + h2) / h1 / (h1 + h2); if (h1 > 2 * h2) { res[2] = 3 / h1; } else { res[2] = (h1 + h2) / h1 / h2; } } else { res[1] = -1.5 / h1; res[2] = -res[1]; } if (rightSide) { res[1] = -res[1]; res[2] = -res[2]; } return res; } // This value is used in the clauses - may not be the returned value final double d = ((2. * h1 + h2) * del1 - h1 * del2) / (h1 + h2); if (Math.signum(d) != Math.signum(del1)) { // again d is now set to zero if (Math.abs(d) < 1e-15) { res[1] = -(2 * h1 + h2) / h1 / (h1 + h2); res[2] = (h1 + h2) / h1 / h2; res[3] = -h1 / h2 / (h1 + h2); } } else if (Math.signum(del1) != Math.signum(del2) && Math.abs(d) > 3. * Math.abs(del1)) { res[0] = 3 * del1; res[1] = -3 / h1; res[2] = -res[1]; } else { res[0] = d; res[1] = -(2 * h1 + h2) / h1 / (h1 + h2); res[2] = (h1 + h2) / h1 / h2; res[3] = -h1 / h2 / (h1 + h2); } if (rightSide) { for (int i = 1; i < 4; i++) { res[i] = -res[i]; } } return res; } @Override public PiecewisePolynomialResult interpolate(final double[] xValues, final double[] yValues) { throw new NotImplementedException(); } @Override public PiecewisePolynomialResult interpolate(double[] xValues, double[][] yValuesMatrix) { throw new NotImplementedException(); } }