class BicubicSplineInterpolator extends PiecewisePolynomialInterpolator2D {

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

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

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

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

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

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

    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];
        }
      }
    }

    DoubleMatrix2D[][] coefMat = new DoubleMatrix2D[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 DoubleMatrix1D diffs = new DoubleMatrix1D(diffsVec);
        final double[] ansVec = ((DoubleMatrix1D) OG_ALGEBRA.multiply(new DoubleMatrix2D(s_invMat), diffs)).getData();

        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[l + m * (order)] / Math.pow((x0Values[i + 1] - x0Values[i]), l) / Math.pow((x1Values[j + 1] - x1Values[j]), m);
            ArgumentChecker.isFalse(Double.isNaN(coefMatTmp[order - l - 1][order - m - 1]), "Too large/small input");
            ArgumentChecker.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);
        ArgumentChecker.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] = new DoubleMatrix2D(coefMatTmp);
      }
    }

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

  public  BicubicSplineInterpolator(final PiecewisePolynomialInterpolator[] method);
  public  BicubicSplineInterpolator(final PiecewisePolynomialInterpolator method);

}