com.opengamma.strata.math.impl.interpolation.BicubicSplineInterpolator Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of strata-math Show documentation
Mathematic support for Strata
There is a newer version: 2.12.46
/*
 * Copyright (C) 2013 - present by OpenGamma Inc. and the OpenGamma group of companies
 *
 * Please see distribution for license.
 */
package com.opengamma.strata.math.impl.interpolation;

import static com.opengamma.strata.math.impl.matrix.MatrixAlgebraFactory.OG_ALGEBRA;

import com.opengamma.strata.collect.ArgChecker;
import com.opengamma.strata.collect.array.DoubleArray;
import com.opengamma.strata.collect.array.DoubleMatrix;
import com.opengamma.strata.math.impl.function.PiecewisePolynomialFunction1D;

/**
 *  Given a set of data (x0Values_i, x1Values_j, yValues_{ij}), derive the piecewise bicubic function,
 *  f(x0,x1) = sum_{i=0}^{3} sum_{j=0}^{3} coefMat_{ij} (x0-x0Values_i)^{3-i} (x1-x1Values_j)^{3-j},
 *  for the region x0Values_i < x0 < x0Values_{i+1}, x1Values_j < x1 < x1Values_{j+1}  such that
 *  f(x0Values_a, x1Values_b) = yValues_{ab} where a={i,i+1}, b={j,j+1}. 
 *  1D piecewise polynomial interpolation methods are called to determine first derivatives and cross derivative at data points
 *  Note that the value of the cross derivative at {ij} is not "accurate" if yValues_{ij} = 0.
 */
public class BicubicSplineInterpolator extends PiecewisePolynomialInterpolator2D {

  private static final double ERROR = 1.e-13;

  private PiecewisePolynomialInterpolator[] _method;
  private static DoubleMatrix INV_MAT;

  static {
    double[][] invMat = new double[16][16];
    invMat[0] = new double[] {1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
    invMat[1] = new double[] {0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
    invMat[2] = new double[] {-3., 3., 0., 0., -2., -1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
    invMat[3] = new double[] {2., -2., 0., 0., 1., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
    invMat[4] = new double[] {0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.};
    invMat[5] = new double[] {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.};
    invMat[6] = new double[] {0., 0., 0., 0., 0., 0., 0., 0., -3., 3., 0., 0., -2., -1., 0., 0.};
    invMat[7] = new double[] {0., 0., 0., 0., 0., 0., 0., 0., 2., -2., 0., 0., 1., 1., 0., 0.};
    invMat[8] = new double[] {-3., 0., 3., 0., 0., 0., 0., 0., -2., 0., -1., 0., 0., 0., 0., 0.};
    invMat[9] = new double[] {0., 0., 0., 0., -3., 0., 3., 0., 0., 0., 0., 0., -2., 0., -1., 0.};
    invMat[10] = new double[] {9., -9., -9., 9., 6., 3., -6., -3., 6., -6., 3., -3., 4., 2., 2., 1.};
    invMat[11] = new double[] {-6., 6., 6., -6., -3., -3., 3., 3., -4., 4., -2., 2., -2., -2., -1., -1.};
    invMat[12] = new double[] {2., 0., -2., 0., 0., 0., 0., 0., 1., 0., 1., 0., 0., 0., 0., 0.};
    invMat[13] = new double[] {0., 0., 0., 0., 2., 0., -2., 0., 0., 0., 0., 0., 1., 0., 1., 0.};
    invMat[14] = new double[] {-6., 6., 6., -6., -4., -2., 4., 2., -3., 3., -3., 3., -2., -1., -2., -1.};
    invMat[15] = new double[] {4., -4., -4., 4., 2., 2., -2., -2., 2., -2., 2., -2., 1., 1., 1., 1.};
    INV_MAT = DoubleMatrix.ofUnsafe(invMat);
  }

  /**
   * Constructor which can take different methods for x0 and x1.
   * @param method Choose 2 of {@link PiecewisePolynomialInterpolator}
   */
  public BicubicSplineInterpolator(final PiecewisePolynomialInterpolator[] method) {
    ArgChecker.notNull(method, "method");
    ArgChecker.isTrue(method.length == 2, "two methods should be chosen");

    _method = new PiecewisePolynomialInterpolator[2];
    for (int i = 0; i < 2; ++i) {
      _method[i] = method[i];
    }
  }

  /**
   * Constructor using the same interpolation method for x0 and x1.
   * @param method {@link PiecewisePolynomialInterpolator}
   */
  public BicubicSplineInterpolator(final PiecewisePolynomialInterpolator method) {
    _method = new PiecewisePolynomialInterpolator[] {method, method };
  }

  @Override
  public PiecewisePolynomialResult2D interpolate(final double[] x0Values, final double[] x1Values, final double[][] yValues) {

    ArgChecker.notNull(x0Values, "x0Values");
    ArgChecker.notNull(x1Values, "x1Values");
    ArgChecker.notNull(yValues, "yValues");

    final int nData0 = x0Values.length;
    final int nData1 = x1Values.length;

    DoubleMatrix yValuesMatrix = DoubleMatrix.copyOf(yValues);
    final PiecewisePolynomialFunction1D func = new PiecewisePolynomialFunction1D();
    double[][] diff0 = new double[nData1][nData0];
    double[][] diff1 = new double[nData0][nData1];
    double[][] cross = new double[nData0][nData1];

    PiecewisePolynomialResult result0 = _method[0].interpolate(x0Values, OG_ALGEBRA.getTranspose(yValuesMatrix).toArray());
    diff0 = func.differentiate(result0, x0Values).toArray();

    PiecewisePolynomialResult result1 = _method[1].interpolate(x1Values, yValuesMatrix.toArray());
    diff1 = func.differentiate(result1, x1Values).toArray();

    final int order = 4;

    for (int i = 0; i < nData0; ++i) {
      for (int j = 0; j < nData1; ++j) {
        if (yValues[i][j] == 0.) {
          if (diff0[j][i] == 0.) {
            cross[i][j] = diff1[i][j];
          } else {
            if (diff1[i][j] == 0.) {
              cross[i][j] = diff0[j][i];
            } else {
              cross[i][j] = Math.signum(diff0[j][i] * diff1[i][j]) * Math.sqrt(Math.abs(diff0[j][i] * diff1[i][j]));
            }
          }
        } else {
          cross[i][j] = diff0[j][i] * diff1[i][j] / yValues[i][j];
        }
      }
    }

    DoubleMatrix[][] coefMat = new DoubleMatrix[nData0 - 1][nData1 - 1];
    for (int i = 0; i < nData0 - 1; ++i) {
      for (int j = 0; j < nData1 - 1; ++j) {
        double[] diffsVec = new double[16];
        for (int l = 0; l < 2; ++l) {
          for (int m = 0; m < 2; ++m) {
            diffsVec[l + 2 * m] = yValues[i + l][j + m];
          }
        }
        for (int l = 0; l < 2; ++l) {
          for (int m = 0; m < 2; ++m) {
            diffsVec[4 + l + 2 * m] = diff0[j + m][i + l];
          }
        }
        for (int l = 0; l < 2; ++l) {
          for (int m = 0; m < 2; ++m) {
            diffsVec[8 + l + 2 * m] = diff1[i + l][j + m];
          }
        }
        for (int l = 0; l < 2; ++l) {
          for (int m = 0; m < 2; ++m) {
            diffsVec[12 + l + 2 * m] = cross[i + l][j + m];
          }
        }
        final DoubleArray diffs = DoubleArray.copyOf(diffsVec);
        final DoubleArray ansVec = ((DoubleArray) OG_ALGEBRA.multiply(INV_MAT, diffs));

        double ref = 0.;
        double[][] coefMatTmp = new double[order][order];
        for (int l = 0; l < order; ++l) {
          for (int m = 0; m < order; ++m) {
            coefMatTmp[order - l - 1][order - m - 1] =
                ansVec.get(l + m * (order)) /
                    Math.pow((x0Values[i + 1] - x0Values[i]), l) /
                    Math.pow((x1Values[j + 1] - x1Values[j]), m);
            ArgChecker.isFalse(Double.isNaN(coefMatTmp[order - l - 1][order - m - 1]), "Too large/small input");
            ArgChecker.isFalse(Double.isInfinite(coefMatTmp[order - l - 1][order - m - 1]), "Too large/small input");
            ref += coefMatTmp[order - l - 1][order - m - 1] * Math.pow((x0Values[i + 1] - x0Values[i]), l) *
                Math.pow((x1Values[j + 1] - x1Values[j]), m);
          }
        }
        final double bound = Math.max(Math.abs(ref) + Math.abs(yValues[i + 1][j + 1]), 0.1);
        ArgChecker.isTrue(Math.abs(ref - yValues[i + 1][j + 1]) < ERROR * bound, "Input is too large/small or data points are too close");
        coefMat[i][j] = DoubleMatrix.copyOf(coefMatTmp);
      }
    }

    return new PiecewisePolynomialResult2D(
        DoubleArray.copyOf(x0Values),
        DoubleArray.copyOf(x1Values),
        coefMat,
        new int[] {order, order});
  }

}