@Override protected double doFirstDerivative(double xValue) { // x-value is less than the x-value of the last node (lowerIndex < intervalCount) int low = lowerBoundIndex(xValue, xValues); int high = low + 1; int n = dataSize - 1; if (low == n) { low = n - 1; high = n; } double delta = xValues[high] - xValues[low]; if (Math.abs(delta) < EPS) { throw new MathException("x data points were not distinct"); } double a = (xValues[high] - xValue) / delta; double b = (xValue - xValues[low]) / delta; double[] y2 = calculateSecondDerivative(xValues, yValues, dataSize, leftFirstDev, rightFirstDev, leftNatural, rightNatural); return (yValues[high] - yValues[low]) / delta + ((-3. * a * a + 1.) * y2[low] + (3. * b * b - 1.) * y2[high]) * delta / 6.; }
@Override protected DoubleArray doParameterSensitivity(double xValue) { // x-value is less than the x-value of the last node (lowerIndex < intervalCount) int low = lowerBoundIndex(xValue, xValues); double[] result = new double[dataSize]; if (low == dataSize - 1) { result[dataSize - 1] = 1.0; return DoubleArray.ofUnsafe(result); } int high = low + 1; double delta = xValues[high] - xValues[low]; double a = (xValues[high] - xValue) / delta; double b = (xValue - xValues[low]) / delta; double c = a * (a * a - 1) * delta * delta / 6.; double d = b * (b * b - 1) * delta * delta / 6.; double[][] y2Sensitivities = getSecondDerivativesSensitivities(xValues, yValues, dataSize, leftNatural, rightNatural); for (int i = 0; i < dataSize; i++) { result[i] = c * y2Sensitivities[low][i] + d * y2Sensitivities[high][i]; } result[low] += a; result[high] += b; return DoubleArray.ofUnsafe(result); }
private static double[] calculateSecondDerivative( double[] xValues, double[] yValues, int dataSize, double leftFirstDev, double rightFirstDev, boolean leftNatural, boolean rightNatural) { double[] deltaX = new double[dataSize - 1]; double[] deltaYOverDeltaX = new double[dataSize - 1]; double[] oneOverDeltaX = new double[dataSize - 1]; for (int i = 0; i < dataSize - 1; i++) { deltaX[i] = xValues[i + 1] - xValues[i]; oneOverDeltaX[i] = 1.0 / deltaX[i]; deltaYOverDeltaX[i] = (yValues[i + 1] - yValues[i]) * oneOverDeltaX[i]; } DoubleMatrix inverseTriDiag = getInverseTridiagonalMatrix(deltaX, leftNatural, rightNatural); DoubleArray rhsVector = getRightVector(deltaYOverDeltaX, leftFirstDev, rightFirstDev, leftNatural, rightNatural); return ((DoubleArray) OG_ALGEBRA.multiply(inverseTriDiag, rhsVector)).toArray(); }
private static double[][] getSecondDerivativesSensitivities( double[] xValues, double[] yValues, int dataSize, boolean leftNatural, boolean rightNatural) { double[] deltaX = new double[dataSize - 1]; double[] deltaYOverDeltaX = new double[dataSize - 1]; double[] oneOverDeltaX = new double[dataSize - 1]; for (int i = 0; i < dataSize - 1; i++) { deltaX[i] = xValues[i + 1] - xValues[i]; oneOverDeltaX[i] = 1.0 / deltaX[i]; deltaYOverDeltaX[i] = (yValues[i + 1] - yValues[i]) * oneOverDeltaX[i]; } DoubleMatrix inverseTriDiag = getInverseTridiagonalMatrix(deltaX, leftNatural, rightNatural); DoubleMatrix rhsMatrix = getRightMatrix(oneOverDeltaX, leftNatural, rightNatural); return ((DoubleMatrix) OG_ALGEBRA.multiply(inverseTriDiag, rhsMatrix)).toArray(); }
@Override protected double doInterpolate(double xValue) { // x-value is less than the x-value of the last node (lowerIndex < intervalCount) int low = lowerBoundIndex(xValue, xValues); int high = low + 1; int n = dataSize - 1; if (low == n) { return yValues[n]; } double delta = xValues[high] - xValues[low]; if (Math.abs(delta) < EPS) { throw new MathException("x data points were not distinct"); } double a = (xValues[high] - xValue) / delta; double b = (xValue - xValues[low]) / delta; double[] y2 = calculateSecondDerivative(xValues, yValues, dataSize, leftFirstDev, rightFirstDev, leftNatural, rightNatural); return a * yValues[low] + b * yValues[high] + (a * (a * a - 1) * y2[low] + b * (b * b - 1) * y2[high]) * delta * delta / 6.; }
@Override public BoundCurveInterpolator bind(DoubleArray xValues, DoubleArray yValues) { return new Bound(xValues, yValues); }
@Override public BoundCurveInterpolator bind( BoundCurveExtrapolator extrapolatorLeft, BoundCurveExtrapolator extrapolatorRight) { return new Bound(this, extrapolatorLeft, extrapolatorRight); } }