/**
* Copyright (C) 2009 - present by OpenGamma Inc. and the OpenGamma group of companies
*
* Please see distribution for license.
*/
package com.opengamma.analytics.math.interpolation.data;
import static com.opengamma.analytics.math.matrix.MatrixAlgebraFactory.OG_ALGEBRA;
import org.apache.commons.lang.ObjectUtils;
import com.opengamma.analytics.math.linearalgebra.InverseTridiagonalMatrixCalculator;
import com.opengamma.analytics.math.linearalgebra.TridiagonalMatrix;
import com.opengamma.analytics.math.matrix.DoubleMatrix1D;
import com.opengamma.analytics.math.matrix.DoubleMatrix2D;
import com.opengamma.util.ArgumentChecker;
/**
*
*/
public class Interpolator1DCubicSplineDataBundle implements Interpolator1DDataBundle {
private final Interpolator1DDataBundle _underlyingData;
private double[] _secondDerivatives;
private double[][] _secondDerivativesSensitivities;
private final double _leftFirstDev;
private final double _rightFirstDev;
private final boolean _leftNatural;
private final boolean _rightNatural;
public Interpolator1DCubicSplineDataBundle(final Interpolator1DDataBundle underlyingData) {
ArgumentChecker.notNull(underlyingData, "underlying data");
_underlyingData = underlyingData;
_leftFirstDev = 0;
_rightFirstDev = 0;
_leftNatural = true;
_rightNatural = true;
}
/**
* Data bundle for a cubic spline
* @param underlyingData the data
* @param leftGrad The gradient of the function at the left most knot. <b>Note: </b>to leave this unspecified (i.e. natural with zero second derivative),
* set the value to Double.POSITIVE_INFINITY
* @param rightGrad The gradient of the function at the right most knot. <b>Note: </b>to leave this unspecified (i.e. natural with zero second derivative),
* set the value to Double.POSITIVE_INFINITY
*/
public Interpolator1DCubicSplineDataBundle(final Interpolator1DDataBundle underlyingData, final double leftGrad, final double rightGrad) {
ArgumentChecker.notNull(underlyingData, "underlying data");
_underlyingData = underlyingData;
if (Double.isInfinite(leftGrad)) {
_leftFirstDev = 0;
_leftNatural = true;
} else {
_leftFirstDev = leftGrad;
_leftNatural = false;
}
if (Double.isInfinite(rightGrad)) {
_rightFirstDev = 0;
_rightNatural = true;
} else {
_rightFirstDev = leftGrad;
_rightNatural = false;
}
}
private double[] calculateSecondDerivative() {
final double[] x = getKeys();
final double[] y = getValues();
final int n = x.length;
final double[] deltaX = new double[n - 1];
final double[] deltaYOverDeltaX = new double[n - 1];
final double[] oneOverDeltaX = new double[n - 1];
for (int i = 0; i < n - 1; i++) {
deltaX[i] = x[i + 1] - x[i];
oneOverDeltaX[i] = 1.0 / deltaX[i];
deltaYOverDeltaX[i] = (y[i + 1] - y[i]) * oneOverDeltaX[i];
}
final DoubleMatrix2D inverseTriDiag = getInverseTridiagonalMatrix(deltaX);
final DoubleMatrix1D rhsVector = getRHSVector(deltaYOverDeltaX);
return ((DoubleMatrix1D) OG_ALGEBRA.multiply(inverseTriDiag, rhsVector)).getData();
}
@Override
public boolean containsKey(final Double key) {
return _underlyingData.containsKey(key);
}
@Override
public Double firstKey() {
return _underlyingData.firstKey();
}
@Override
public Double firstValue() {
return _underlyingData.firstValue();
}
@Override
public Double get(final Double key) {
return _underlyingData.get(key);
}
@Override
public InterpolationBoundedValues getBoundedValues(final Double key) {
return _underlyingData.getBoundedValues(key);
}
@Override
public double[] getKeys() {
return _underlyingData.getKeys();
}
@Override
public int getLowerBoundIndex(final Double value) {
return _underlyingData.getLowerBoundIndex(value);
}
@Override
public Double getLowerBoundKey(final Double value) {
return _underlyingData.getLowerBoundKey(value);
}
@Override
public double[] getValues() {
return _underlyingData.getValues();
}
@Override
public Double higherKey(final Double key) {
return _underlyingData.higherKey(key);
}
@Override
public Double higherValue(final Double key) {
return _underlyingData.higherValue(key);
}
@Override
public Double lastKey() {
return _underlyingData.lastKey();
}
@Override
public Double lastValue() {
return _underlyingData.lastValue();
}
@Override
public int size() {
return _underlyingData.size();
}
public double[] getSecondDerivatives() {
if (_secondDerivatives == null) {
_secondDerivatives = calculateSecondDerivative();
}
return _secondDerivatives;
}
//TODO not ideal that it recomputes the inverse matrix
public double[][] getSecondDerivativesSensitivities() {
if (_secondDerivativesSensitivities == null) {
final double[] x = getKeys();
final double[] y = getValues();
final int n = x.length;
final double[] deltaX = new double[n - 1];
final double[] deltaYOverDeltaX = new double[n - 1];
final double[] oneOverDeltaX = new double[n - 1];
for (int i = 0; i < n - 1; i++) {
deltaX[i] = x[i + 1] - x[i];
oneOverDeltaX[i] = 1.0 / deltaX[i];
deltaYOverDeltaX[i] = (y[i + 1] - y[i]) * oneOverDeltaX[i];
}
final DoubleMatrix2D inverseTriDiag = getInverseTridiagonalMatrix(deltaX);
final DoubleMatrix2D rhsMatrix = getRHSMatrix(oneOverDeltaX);
_secondDerivativesSensitivities = ((DoubleMatrix2D) OG_ALGEBRA.multiply(inverseTriDiag, rhsMatrix)).getData();
}
return _secondDerivativesSensitivities;
}
private DoubleMatrix2D getRHSMatrix(final double[] oneOverDeltaX) {
final int n = oneOverDeltaX.length + 1;
final double[][] res = new double[n][n];
for (int i = 1; i < n - 1; i++) {
res[i][i - 1] = oneOverDeltaX[i - 1];
res[i][i] = -oneOverDeltaX[i] - oneOverDeltaX[i - 1];
res[i][i + 1] = oneOverDeltaX[i];
}
if (!_leftNatural) {
res[0][0] = oneOverDeltaX[0];
res[0][1] = -oneOverDeltaX[0];
}
if (!_rightNatural) {
res[n - 1][n - 1] = -oneOverDeltaX[n - 2];
res[n - 2][n - 2] = oneOverDeltaX[n - 2];
}
return new DoubleMatrix2D(res);
}
private DoubleMatrix1D getRHSVector(final double[] deltaYOverDeltaX) {
final int n = deltaYOverDeltaX.length + 1;
final double[] res = new double[n];
for (int i = 1; i < n - 1; i++) {
res[i] = deltaYOverDeltaX[i] - deltaYOverDeltaX[i - 1];
}
if (!_leftNatural) {
res[0] = _leftFirstDev - deltaYOverDeltaX[0];
}
if (!_rightNatural) {
res[n - 1] = _rightFirstDev - deltaYOverDeltaX[n - 2];
}
return new DoubleMatrix1D(res);
}
private DoubleMatrix2D getInverseTridiagonalMatrix(final double[] deltaX) {
final InverseTridiagonalMatrixCalculator invertor = new InverseTridiagonalMatrixCalculator();
final int n = deltaX.length + 1;
final double[] a = new double[n];
final double[] b = new double[n - 1];
final double[] c = new double[n - 1];
for (int i = 1; i < n - 1; i++) {
a[i] = (deltaX[i - 1] + deltaX[i]) / 3.0;
b[i] = deltaX[i] / 6.0;
c[i - 1] = deltaX[i - 1] / 6.0;
}
// Boundary condition
if (_leftNatural) {
a[0] = 1.0;
b[0] = 0.0;
} else {
a[0] = -deltaX[0] / 3.0;
b[0] = deltaX[0] / 6.0;
}
if (_rightNatural) {
a[n - 1] = 1.0;
c[n - 2] = 0.0;
} else {
a[n - 1] = deltaX[n - 2] / 3.0;
c[n - 2] = deltaX[n - 2] / 6.0;
}
final TridiagonalMatrix tridiagonal = new TridiagonalMatrix(a, b, c);
return invertor.evaluate(tridiagonal);
}
@Override
public void setYValueAtIndex(final int index, final double y) {
ArgumentChecker.notNegative(index, "index");
if (index >= size()) {
throw new IllegalArgumentException("Index was greater than number of data points");
}
_underlyingData.setYValueAtIndex(index, y);
_secondDerivatives = null;
_secondDerivativesSensitivities = null;
}
@Override
public int hashCode() {
final int prime = 31;
int result = 1;
long temp;
temp = Double.doubleToLongBits(_leftFirstDev);
result = prime * result + (int) (temp ^ (temp >>> 32));
result = prime * result + 1237;
temp = Double.doubleToLongBits(_rightFirstDev);
result = prime * result + (int) (temp ^ (temp >>> 32));
result = prime * result + 1237;
result = prime * result + ((_underlyingData == null) ? 0 : _underlyingData.hashCode());
return result;
}
@Override
public boolean equals(final Object obj) {
if (this == obj) {
return true;
}
if (obj == null) {
return false;
}
if (getClass() != obj.getClass()) {
return false;
}
final Interpolator1DCubicSplineDataBundle other = (Interpolator1DCubicSplineDataBundle) obj;
if (!ObjectUtils.equals(_underlyingData, other._underlyingData)) {
return false;
}
if (Double.doubleToLongBits(_leftFirstDev) != Double.doubleToLongBits(other._leftFirstDev)) {
return false;
}
if (_leftNatural != other._leftNatural) {
return false;
}
if (Double.doubleToLongBits(_rightFirstDev) != Double.doubleToLongBits(other._rightFirstDev)) {
return false;
}
if (_rightNatural != other._rightNatural) {
return false;
}
return true;
}
}