/**
* 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 org.apache.commons.lang.NotImplementedException;
import org.apache.commons.lang.Validate;
import com.opengamma.analytics.math.function.PiecewisePolynomialFunction1D;
import com.opengamma.analytics.math.function.PiecewisePolynomialWithSensitivityFunction1D;
import com.opengamma.analytics.math.interpolation.data.ArrayInterpolator1DDataBundle;
import com.opengamma.analytics.math.interpolation.data.Interpolator1DDataBundle;
import com.opengamma.analytics.math.interpolation.data.Interpolator1DPiecewisePoynomialDataBundle;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.util.ArgumentChecker;
import com.opengamma.util.ParallelArrayBinarySort;
/**
* Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) for use in yield curves. The yield curve r(t) is such that P(0,t) = exp(-t*r(t)) is the
* discount factor for time t. Here we actually interpolate on the quantity f(t) = t*r(t) (which must be monotonically increasing) rather than r(t) itself.
* However the inputs are still the set {t_i,r_i}, and the interpolate method returns r(t) rather than f(t). If t_0 != 0 an extra data point at zero is inserted
* such that t_0 = 0 (the value of r_0 is irrelevant)
*
*/
public class PCHIPYieldCurveInterpolator1D extends Interpolator1D {
/** Serialization version */
private static final long serialVersionUID = 1L;
// TODO have options on method
private static final PiecewisePolynomialFunction1D FUNC = new PiecewisePolynomialWithSensitivityFunction1D();
@Override
public Double interpolate(final Interpolator1DDataBundle data, final Double value) {
final double eps = 1e-10;
Validate.notNull(value, "value");
Validate.notNull(data, "data bundle");
Validate.isTrue(data instanceof Interpolator1DPiecewisePoynomialDataBundle);
final Interpolator1DPiecewisePoynomialDataBundle polyData = (Interpolator1DPiecewisePoynomialDataBundle) data;
ArgumentChecker.isFalse(value < 0, "value must be zero or positive");
if (value == 0) {
return 0.0;
}
// TODO direct evaluation using coefficients
final double t = Math.max(eps, value); // Avoid divide by zero
final DoubleMatrix1D res = FUNC.evaluate(polyData.getPiecewisePolynomialResultsWithSensitivity(), t);
return res.getEntry(0) / t;
}
@Override
public double firstDerivative(final Interpolator1DDataBundle data, final Double value) {
final double eps = 1e-10;
Validate.notNull(value, "value");
Validate.notNull(data, "data bundle");
Validate.isTrue(data instanceof Interpolator1DPiecewisePoynomialDataBundle);
final Interpolator1DPiecewisePoynomialDataBundle polyData = (Interpolator1DPiecewisePoynomialDataBundle) data;
ArgumentChecker.isFalse(value < 0, "value must be zero or positive");
if (value == 0) {
return 0.0;
}
final double t = Math.max(eps, value);
final DoubleMatrix1D resValue = FUNC.evaluate(polyData.getPiecewisePolynomialResultsWithSensitivity(), t);
final DoubleMatrix1D resDerivative = FUNC.differentiate(polyData.getPiecewisePolynomialResultsWithSensitivity(), t);
return (resDerivative.getEntry(0) - resValue.getEntry(0) / t) / t;
}
@Override
public double[] getNodeSensitivitiesForValue(final Interpolator1DDataBundle data, final Double value) {
throw new NotImplementedException();
}
@Override
public Interpolator1DDataBundle getDataBundle(final double[] x, final double[] y) {
final int n = x.length;
ArgumentChecker.isTrue(n == y.length, "x and y different lengths");
final double[] xSrt = new double[n];
final double[] ySrt = new double[n];
System.arraycopy(x, 0, xSrt, 0, n);
System.arraycopy(y, 0, ySrt, 0, n);
ParallelArrayBinarySort.parallelBinarySort(xSrt, ySrt);
return getDataBundleFromSortedArrays(xSrt, ySrt);
}
@Override
public Interpolator1DDataBundle getDataBundleFromSortedArrays(final double[] x, final double[] y) {
final int n = x.length;
ArgumentChecker.isTrue(n == y.length, "x and y different lengths");
ArgumentChecker.isTrue(n == y.length, "x and y different lengths");
ArgumentChecker.isTrue(x[0] >= 0, "first x-values cannot be negative");
double[] xx = x;
double[] xy;
if (xx[0] > 0) {
final double[] temp = new double[n + 1];
System.arraycopy(xx, 0, temp, 1, n);
xx = temp;
xy = new double[n + 1];
for (int i = 1; i <= n; i++) {
xy[i] = xx[i] * y[i - 1];
}
} else {
xy = new double[n];
for (int i = 1; i < n; i++) {
xy[i] = x[i] * y[i];
}
}
// final PiecewisePolynomialResult poly = BASE.interpolate(xx, xy);
return new Interpolator1DPiecewisePoynomialDataBundle(new ArrayInterpolator1DDataBundle(xx, xy, true), new PiecewiseCubicHermiteSplineInterpolatorWithSensitivity());
}
}