/**
* 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.opengamma.analytics.math.interpolation.data.Interpolator1DPiecewisePoynomialWithExtraKnotsDataBundle;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.util.ArgumentChecker;
import com.opengamma.util.ParallelArrayBinarySort;
/**
* Shape-preserving C2 cubic spline interpolation based on
* S. Pruess "Shape preserving C2 cubic spline interpolation"
* IMA Journal of Numerical Analysis (1993) 13 (4): 493-507.
* where two extra knots are introduced between adjacent data points
* As the position of the new knots are data dependent, the matrix form of yValues producing multi-splines is not relevant
*/
public class ShapePreservingCubicSplineInterpolator extends PiecewisePolynomialInterpolator {
private static final double INF = 1. / 0.;
private static final double ERROR = 1.e-12;
@Override
public PiecewisePolynomialResult interpolate(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 > 2, "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");
}
for (int i = 0; i < nDataPts; ++i) {
for (int j = i + 1; j < nDataPts; ++j) {
ArgumentChecker.isFalse(xValues[i] == xValues[j], "xValues should be distinct");
}
}
final double[] xValuesSrt = Arrays.copyOf(xValues, nDataPts);
final double[] yValuesSrt = Arrays.copyOf(yValues, nDataPts);
ParallelArrayBinarySort.parallelBinarySort(xValuesSrt, yValuesSrt);
final double[] intervals = intervalsCalculator(xValuesSrt);
final double[] slopes = slopesCalculator(yValuesSrt, intervals);
final double[] beta = betaCalculator(slopes);
double[] first = firstDiffFinder(intervals, slopes);
double[] rValues = rValuesCalculator(slopes, first);
boolean correctSign = false;
int it = 0;
while (correctSign == false) {
correctSign = signChecker(beta, rValues);
if (correctSign == false) {
first = firstDiffSweep(intervals, slopes, beta, first);
rValues = rValuesCalculator(slopes, first);
}
++it;
if (it > 10) {
throw new IllegalArgumentException("Spline is not found!");
}
}
final double[] second = secondDiffFinder(intervals, beta, rValues);
final double[] tau = tauFinder(intervals, slopes, beta, first, second);
final double[] knots = knotsProvider(xValuesSrt, intervals, tau);
final double[][] coefMatrix = solve(yValuesSrt, intervals, slopes, first, second, tau);
for (int i = 0; i < coefMatrix.length; ++i) {
double ref = 0.;
final double interval = knots[i + 1] - knots[i];
for (int j = 0; j < 4; ++j) {
ref += coefMatrix[i][j] * Math.pow(interval, 3 - j);
ArgumentChecker.isFalse(Double.isNaN(coefMatrix[i][j]), "Too large input");
ArgumentChecker.isFalse(Double.isInfinite(coefMatrix[i][j]), "Too large input");
}
final double yVal = i == coefMatrix.length - 1 ? yValues[nDataPts - 1] : coefMatrix[i + 1][3];
final double bound = Math.max(Math.abs(ref) + Math.abs(yVal), 1.e-1);
ArgumentChecker.isTrue(Math.abs(ref - yVal) < ERROR * bound, "Input is too large/small or data points are too close");
}
return new PiecewisePolynomialResult(new DoubleMatrix1D(knots), new DoubleMatrix2D(coefMatrix), 4, 1);
}
@Override
public PiecewisePolynomialResult interpolate(final double[] xValues, final double[][] yValuesMatrix) {
throw new IllegalArgumentException("Method with multidimensional yValues is not supported");
}
/**
* Since this interpolation method introduces new breakpoints in certain cases, {@link PiecewisePolynomialResultsWithSensitivity} is not well-defined
* Instead the node sensitivity is computed in {@link MonotoneConvexSplineInterpolator1D} via {@link Interpolator1DPiecewisePoynomialWithExtraKnotsDataBundle}
* @param xValues The xValues
* @param yValues The yValues
* @return NotImplementedException
*/
@Override
public PiecewisePolynomialResultsWithSensitivity interpolateWithSensitivity(final double[] xValues, final double[] yValues) {
throw new NotImplementedException();
}
/**
* @param xValues
* @return Intervals of xValues, ( xValues_i - xValues_{i-1} )
*/
private double[] intervalsCalculator(final double[] xValues) {
final int nDataPts = xValues.length;
final double[] intervals = new double[nDataPts - 1];
for (int i = 0; i < nDataPts - 1; ++i) {
intervals[i] = xValues[i + 1] - xValues[i];
}
return intervals;
}
/**
* @param yValues Y values of data
* @param intervals Intervals of x data
* @return Slopes
*/
private double[] slopesCalculator(final double[] yValues, final double[] intervals) {
final int nDataPts = yValues.length;
final double[] slopes = new double[nDataPts - 1];
for (int i = 0; i < nDataPts - 1; ++i) {
slopes[i] = (yValues[i + 1] - yValues[i]) / intervals[i];
}
return slopes;
}
/**
* @param intervals
* @param slopes
* @return First derivative at knots
*/
private double[] firstDiffFinder(final double[] intervals, final double[] slopes) {
final int nInts = intervals.length;
final double[] res = new double[nInts + 1];
res[0] = endpointFirst(intervals[0], intervals[1], slopes[0], slopes[1]);
res[nInts] = endpointFirst(intervals[nInts - 1], intervals[nInts - 2], slopes[nInts - 1], slopes[nInts - 2]);
for (int i = 1; i < nInts; ++i) {
if (Math.signum(slopes[i]) != Math.signum(slopes[i - 1]) | (slopes[i] == 0 | slopes[i - 1] == 0)) {
res[i] = 0.;
} else {
final double den1 = 2. * intervals[i] + intervals[i - 1];
final double den2 = intervals[i] + 2. * intervals[i - 1];
res[i] = 3. * (intervals[i] + intervals[i - 1]) / (den1 / slopes[i - 1] + den2 / slopes[i]);
}
}
return res;
}
/**
* Estimate first derivatives at endpoints
* @param ints1 First (last) interval
* @param ints2 Second (second last) Interval
* @param grads1 First (last) slope
* @param grads2 Second (second last) slope
* @return Slope at the first (last) data point
*/
private double endpointFirst(final double ints1, final double ints2, final double grads1, final double grads2) {
final double val = (2. * ints1 + ints2) * grads1 / (ints1 + ints2) - ints1 * grads2 / (ints1 + ints2);
if (Math.signum(val) != Math.signum(grads1)) {
return 0.;
}
if (Math.signum(grads1) != Math.signum(grads2) && Math.abs(val) > 3. * Math.abs(grads1)) {
return 3. * grads1;
}
return val;
}
/**
* @param slopes
* @return sign(slopes_{i + 1} - slopes_i)
*/
private double[] betaCalculator(final double[] slopes) {
final int nSlopes = slopes.length;
final double[] res = new double[nSlopes + 1];
for (int i = 0; i < nSlopes - 1; ++i) {
res[i + 1] = Math.signum(slopes[i + 1] - slopes[i]);
}
res[0] = res[1];
res[nSlopes] = res[nSlopes - 1];
return res;
}
/**
* In the notation i =1,2,...,N-1,
* R_{2*i-1} = 6*slopes_i - 4*first_i - 2*first_{i+1}
* R_{2*i} = - 6*slopes_i + 2*first_i + 4*first_{i+1}
* @param slopes
* @param first First derivatives
* @return R functions
*/
private double[] rValuesCalculator(final double[] slopes, final double[] first) {
final int nData = first.length;
final double[] res = new double[2 * nData];
for (int i = 1; i < nData - 1; ++i) {
res[2 * i] = 2. * first[i - 1] + 4. * first[i] - 6. * slopes[i - 1];
res[2 * i - 1] = -4. * first[i - 1] - 2. * first[i] + 6. * slopes[i - 1];
if (Math.abs(res[2 * i]) <= 0.1 * ERROR) {
res[2 * i] = 0.;
}
if (Math.abs(res[2 * i - 1]) <= 0.1 * ERROR) {
res[2 * i - 1] = 0.;
}
}
res[0] = INF;
res[2 * nData - 1] = INF;
return res;
}
/**
* Check beta_i * R_{2*i-1} \geq 0 and beta_{i+1} *R_{2*i} \geq 0
* @param beta
* @param rValues
* @return True if the two inequalities are satisfied
*/
private boolean signChecker(final double[] beta, final double[] rValues) {
final int nData = beta.length;
for (int i = 1; i < nData - 2; ++i) {
if (beta[i] * rValues[2 * i + 1] < 0 | beta[i + 1] * rValues[2 * i + 2] < 0) {
return false;
}
}
return true;
}
/**
* As the double sweep algorithm determines only intervals [minTmp, maxTmp] which should contain first derivative, choice of the values of first derivatives is not unique.
* Infeasibility is returned in some cases, which is resolved by another choice of first derivatives within the allowed intervals, [minTmp, maxTmp].
* @param intervals
* @param slopes
* @param beta
* @param first
* @return First derivatives satisfying convexity conditions
*/
private double[] firstDiffSweep(final double[] intervals, final double[] slopes, final double[] beta, final double[] first) {
final int nData = intervals.length + 1;
final double[] res = new double[nData];
double minVal = 0.;
double maxVal = 0.;
if (beta[0] > 0) {
minVal = 3. * slopes[0] - 2. * slopes[1];
maxVal = slopes[0];
} else {
minVal = slopes[0];
maxVal = 3. * slopes[0] - 2. * slopes[1];
}
for (int i = 1; i < nData; ++i) {
double minTmp = 0.;
double maxTmp = 0.;
if (beta[i - 1] == 0.) {
if (beta[i] == 0.) {
minTmp = -INF;
maxTmp = INF;
} else {
if (beta[i] > 0.) {
minTmp = 0.5 * (3. * slopes[i - 1] - maxVal);
maxTmp = INF;
} else {
minTmp = -INF;
maxTmp = Math.min(slopes[i - 1], 0.5 * (3. * slopes[i - 1] - minVal));
}
}
} else {
if (beta[i - 1] > 0.) {
if (beta[i] == 0.) {
minTmp = slopes[i - 1];
maxTmp = 3. * slopes[i - 1] - 2. * minVal;
} else {
if (beta[i] > 0.) {
minTmp = Math.max(slopes[i - 1], 0.5 * (3. * slopes[i - 1] - maxVal));
maxTmp = 3. * slopes[i - 1] - 2. * minVal;
} else {
minTmp = -INF;
maxTmp = Math.min(3. * slopes[i - 1] - 2. * minVal, 0.5 * (3. * slopes[i - 1] - minVal));
}
}
} else {
if (beta[i] == 0) {
minTmp = 3. * slopes[i - 1] - 2. * minVal;
maxTmp = INF;
} else {
if (beta[i] > 0.) {
minTmp = Math.max(3. * slopes[i - 1] - 2. * maxVal, 0.5 * (3. * slopes[i - 1] - maxVal));
maxTmp = INF;
} else {
minTmp = 3. * slopes[i - 1] - 2. * maxVal;
maxTmp = Math.min(slopes[i - 1], 0.5 * (3. * slopes[i - 1] - minVal));
}
}
}
}
minVal = minTmp;
maxVal = maxTmp;
if (minTmp != -INF && maxTmp != INF) {
res[i] = 0.5 * (minTmp + maxTmp);
} else {
if (first[i] < minTmp) {
res[i] = minTmp != INF ? minTmp : first[i];
} else {
if (first[i] > maxTmp) {
res[i] = maxTmp != -INF ? maxTmp : first[i];
} else {
res[i] = first[i];
}
}
}
}
if (minVal > maxVal) {
res[nData - 1] = first[nData - 1];
} else {
if (first[nData - 1] < minVal) {
res[nData - 1] = minVal != INF ? minVal : first[nData - 1];
} else {
if (first[nData - 1] > maxVal) {
res[nData - 1] = maxVal != -INF ? maxVal : first[nData - 1];
} else {
res[nData - 1] = first[nData - 1];
}
}
}
double minTmp = 0.;
double maxTmp = 0.;
for (int i = nData - 2; i > -1; --i) {
minTmp = 0.;
maxTmp = 0.;
if (beta[i] == 0.) {
if (beta[i + 1] == 0.) {
minTmp = -INF;
maxTmp = INF;
} else {
if (beta[i + 1] > 0.) {
minTmp = 3. * slopes[i] - 2. * res[i + 1];
maxTmp = INF;
} else {
minTmp = -INF;
maxTmp = 3. * slopes[i] - 2. * res[i + 1];
}
}
} else {
if (beta[i] > 0.) {
if (beta[i + 1] == 0.) {
minTmp = -INF;
maxTmp = 0.5 * (3. * slopes[i] - res[i + 1]);
} else {
if (beta[i + 1] > 0.) {
minTmp = 3. * slopes[i] - 2. * res[i + 1];
maxTmp = 0.5 * (3. * slopes[i] - res[i + 1]);
} else {
minTmp = -INF;
maxTmp = Math.min(3. * slopes[i] - 2. * res[i + 1], 0.5 * (3. * slopes[i] - res[i + 1]));
}
}
} else {
if (beta[i + 1] == 0.) {
minTmp = 0.5 * (3. * slopes[i] - res[i + 1]);
maxTmp = INF;
} else {
if (beta[i + 1] > 0.) {
minTmp = Math.max(3. * slopes[i] - 2. * res[i + 1], 0.5 * (3. * slopes[i] - res[i + 1]));
maxTmp = INF;
} else {
minTmp = 0.5 * (3. * slopes[i] - res[i + 1]);
maxTmp = 3. * slopes[i] - 2. * res[i + 1];
}
}
}
}
if (minTmp > maxTmp) {
throw new IllegalArgumentException("Local monotonicity can not be preserved");
}
if (res[i] < minTmp) {
res[i] = minTmp != INF ? minTmp : res[i];
} else {
if (res[i] > maxTmp) {
res[i] = maxTmp != -INF ? maxTmp : res[i];
}
}
}
return res;
}
/**
* @param intervals
* @param beta
* @param rValues
* @return Second derivatives satisfying local monotonicity
*/
private double[] secondDiffFinder(final double[] intervals, final double[] beta, final double[] rValues) {
final int nData = intervals.length + 1;
final double[] res = new double[nData];
for (int i = 1; i < nData - 1; ++i) {
res[i] = (rValues[2 * i + 1] > 0 && rValues[2 * i] > 0) ? beta[i] * Math.min(beta[i] * rValues[2 * i + 1] / intervals[i], beta[i] * rValues[2 * i] / intervals[i - 1]) : 0.;
}
res[0] = rValues[1] > 0 ? beta[0] * rValues[1] / intervals[0] : 0.;
res[nData - 1] = rValues[2 * (nData - 1)] > 0 ? beta[nData - 1] * rValues[2 * (nData - 1)] / intervals[nData - 2] : 0.;
return res;
}
/**
* Extra knots are introduced at xValues[i] + tau[i] * intervals[i] and xValues[i + 1] - tau[i] * intervals[i]
* @param intervals
* @param slopes
* @param beta
* @param first
* @param second
* @return tau
*/
private double[] tauFinder(final double[] intervals, final double[] slopes, final double[] beta, final double[] first, final double[] second) {
final int nData = intervals.length + 1;
final double[] res = new double[nData - 1];
Arrays.fill(res, 1. / 3.);
for (int i = 1; i < nData - 2; ++i) {
boolean ineq1 = false;
boolean ineq2 = false;
final double bound1 = 6. * slopes[i] * beta[i];
final double bound2 = 6. * slopes[i] * beta[i + 1];
double ref1 = (4. * first[i] + 2. * first[i + 1] + intervals[i] * second[i] * res[i] * (2. - res[i]) - intervals[i] * second[i + 1] * res[i] * (1. - res[i])) * beta[i];
double ref2 = (2. * first[i] + 4. * first[i + 1] + intervals[i] * second[i] * res[i] * (1. - res[i]) - intervals[i] * second[i + 1] * res[i] * (2. - res[i])) * beta[i + 1];
while (ineq1 == false) {
if (ref1 - ERROR * Math.abs(ref1) <= bound1 + ERROR * Math.abs(bound1)) {
ineq1 = true;
} else {
res[i] *= 0.8;
}
if (res[i] < ERROR / 100.) {
throw new IllegalArgumentException("Spline is not found");
}
ref1 = (4. * first[i] + 2. * first[i + 1] + intervals[i] * second[i] * res[i] * (2. - res[i]) - intervals[i] * second[i + 1] * res[i] * (1. - res[i])) * beta[i];
}
while (ineq2 == false) {
if (ref2 + ERROR * Math.abs(ref2) >= bound2 - ERROR * Math.abs(bound2)) {
ineq2 = true;
} else {
res[i] *= 0.8;
}
if (res[i] < ERROR / 100.) {
throw new IllegalArgumentException("Spline is not found");
}
ref2 = (2. * first[i] + 4. * first[i + 1] + intervals[i] * second[i] * res[i] * (1. - res[i]) - intervals[i] * second[i + 1] * res[i] * (2. - res[i])) * beta[i + 1];
}
}
return res;
}
/**
* @param xValues
* @param intervals
* @param tau
* @return {... , xValues[i], xValues[i] + tau[i] * intervals[i], xValues[i + 1] - tau[i] * intervals[i], xValues[i + 1], ...}
*/
private double[] knotsProvider(final double[] xValues, final double[] intervals, final double[] tau) {
final int nData = xValues.length;
final double[] res = new double[3 * nData - 2];
for (int i = 0; i < nData - 1; ++i) {
res[3 * i] = xValues[i];
res[3 * i + 1] = xValues[i] + tau[i] * intervals[i];
res[3 * i + 2] = xValues[i + 1] - tau[i] * intervals[i];
}
res[3 * (nData - 1)] = xValues[nData - 1];
return res;
}
/**
* Determine value, first derivative, second derivative and third derivative of interpolant at extra knots
* @param yValues
* @param intervals
* @param slopes
* @param first
* @param second
* @param tau
* @return Coefficient matrix whose i-th row vector is {a_n, a_{n-1}, ... } of f(x) = a_n * (x-x_i)^n + a_{n-1} * (x-x_i)^{n-1} +... for the i-th interval
*/
private double[][] solve(final double[] yValues, final double[] intervals, final double[] slopes, final double[] first, final double[] second, final double[] tau) {
final int nData = yValues.length;
final double[][] res = new double[3 * (nData - 1)][4];
final double[] secNewKnots1 = new double[nData - 1];
final double[] secNewKnots2 = new double[nData - 1];
for (int i = 0; i < nData - 1; ++i) {
secNewKnots1[i] = 6. * slopes[i] / intervals[i] / (1. - tau[i]) - 4. * first[i] / intervals[i] / (1. - tau[i]) - 2. * first[i + 1] / intervals[i] / (1. - tau[i]) - tau[i] * (2. - tau[i]) *
second[i] / (1. - tau[i]) + tau[i] * second[i + 1];
secNewKnots2[i] = -6. * slopes[i] / intervals[i] / (1. - tau[i]) + 4. * first[i + 1] / intervals[i] / (1. - tau[i]) + 2. * first[i] / intervals[i] / (1. - tau[i]) - tau[i] * (2. - tau[i]) *
second[i + 1] / (1. - tau[i]) + tau[i] * second[i];
}
for (int i = 0; i < nData - 1; ++i) {
res[3 * i][0] = (secNewKnots1[i] - second[i]) / 6. / tau[i] / intervals[i];
res[3 * i][1] = 0.5 * second[i];
res[3 * i][2] = first[i];
res[3 * i][3] = yValues[i];
res[3 * i + 1][0] = (secNewKnots2[i] - secNewKnots1[i]) / 6. / (1. - 2. * tau[i]) / intervals[i];
res[3 * i + 1][1] = 0.5 * secNewKnots1[i];
res[3 * i + 1][2] = first[i] + (second[i] + secNewKnots1[i]) * tau[i] * intervals[i] * 0.5;
res[3 * i + 1][3] = yValues[i] + tau[i] * intervals[i] * first[i] + (2. * second[i] + secNewKnots1[i]) * tau[i] * intervals[i] * tau[i] * intervals[i] / 6.;
res[3 * i + 2][0] = (second[i + 1] - secNewKnots2[i]) / 6. / tau[i] / intervals[i];
res[3 * i + 2][1] = 0.5 * secNewKnots2[i];
res[3 * i + 2][2] = first[i + 1] - (second[i + 1] + secNewKnots2[i]) * tau[i] * intervals[i] * 0.5;
res[3 * i + 2][3] = yValues[i + 1] - tau[i] * intervals[i] * first[i + 1] + (2. * second[i + 1] + secNewKnots2[i]) * tau[i] * intervals[i] * tau[i] * intervals[i] / 6.;
}
return res;
}
}