}
}
PiecewisePolynomialInterpolator method = new NaturalSplineInterpolator();
MeshgridInterpolator2D interp = new MeshgridInterpolator2D(method);
PiecewisePolynomialFunction1D func = new PiecewisePolynomialFunction1D();
final int n0Keys = 10 * n0Data + 1;
final int n1Keys = 10 * n1Data + 1;
final double eps = 1.e-6;
double[] x0Keys = new double[n0Keys];
double[] x1Keys = new double[n1Keys];
double[] x0KeysUp = new double[n0Keys];
double[] x1KeysUp = new double[n1Keys];
double[] x0KeysDown = new double[n0Keys];
double[] x1KeysDown = new double[n1Keys];
for (int i = 0; i < n0Keys; ++i) {
x0Keys[i] = x0Values[0] + (x0Values[n0Data - 1] - x0Values[0]) / (n0Keys - 1) * i;
x0KeysUp[i] = x0Keys[i] == 0. ? eps : x0Keys[i] * (1. + eps);
x0KeysDown[i] = x0Keys[i] == 0. ? -eps : x0Keys[i] * (1. - eps);
}
for (int i = 0; i < n1Keys; ++i) {
x1Keys[i] = x1Values[0] + (x1Values[n1Data - 1] - x1Values[0]) / (n1Keys - 1) * i;
x1KeysUp[i] = x1Keys[i] == 0. ? eps : x1Keys[i] * (1. + eps);
x1KeysDown[i] = x1Keys[i] == 0. ? -eps : x1Keys[i] * (1. - eps);
}
final double[][] knotsRes = interp.interpolate(x0Values, x1Values, yValues, x0Values, x1Values).getData();
final double[][] res = interp.interpolate(x0Values, x1Values, yValues, x0Keys, x1Keys).getData();
final double[][] resX0 = interp.differentiateX0(x0Values, x1Values, yValues, x0Keys, x1Keys).getData();
final double[][] resX1 = interp.differentiateX1(x0Values, x1Values, yValues, x0Keys, x1Keys).getData();
final double[][] resX0Twice = interp.differentiateX0Twice(x0Values, x1Values, yValues, x0Keys, x1Keys).getData();
final double[][] resX1Twice = interp.differentiateX1Twice(x0Values, x1Values, yValues, x0Keys, x1Keys).getData();
final double[][] resTran = OG_ALGEBRA.getTranspose(interp.interpolate(x1Values, x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData(), x1Keys, x0Keys)).getData();
final double[][] resX0Tran = OG_ALGEBRA.getTranspose(interp.differentiateX1(x1Values, x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData(), x1Keys, x0Keys)).getData();
final double[][] resX1Tran = OG_ALGEBRA.getTranspose(interp.differentiateX0(x1Values, x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData(), x1Keys, x0Keys)).getData();
final double[][] resX0TwiceTran = OG_ALGEBRA.getTranspose(interp.differentiateX1Twice(x1Values, x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData(), x1Keys, x0Keys))
.getData();
final double[][] resX1TwiceTran = OG_ALGEBRA.getTranspose(interp.differentiateX0Twice(x1Values, x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData(), x1Keys, x0Keys))
.getData();
for (int i = 0; i < n0Data; ++i) {
for (int j = 0; j < n1Data; ++j) {
double ref = Math.abs(yValues[i][j]) < 0.1 * EPS ? 0.1 : Math.abs(yValues[i][j]);
assertEquals(knotsRes[i][j], yValues[i][j], EPS * ref);
}
}
{
double ref = Math.abs(yValues[1][2]) < 0.1 * EPS ? 0.1 : Math.abs(yValues[1][2]);
assertEquals(interp.interpolate(x0Values, x1Values, yValues, x0Values[1], x1Values[2]), yValues[1][2], EPS * ref);
}
for (int i = 0; i < n0Keys; ++i) {
for (int j = 0; j < n1Keys; ++j) {
double ref = Math.abs(res[i][j]) < 0.1 * EPS ? 0.1 : Math.abs(res[i][j]);
assertEquals(res[i][j], resTran[i][j], EPS * ref);
ref = Math.abs(resX0[i][j]) < 0.1 * EPS ? 0.1 : Math.abs(resX0[i][j]);
assertEquals(resX0[i][j], resX0Tran[i][j], EPS * ref);
ref = Math.abs(resX1[i][j]) < 0.1 * EPS ? 0.1 : Math.abs(resX1[i][j]);
assertEquals(resX1[i][j], resX1Tran[i][j], EPS * ref);
ref = Math.abs(resX0Twice[i][j]) < 0.1 * EPS ? 0.1 : Math.abs(resX0Twice[i][j]);
assertEquals(resX0Twice[i][j], resX0TwiceTran[i][j], EPS * ref);
ref = Math.abs(resX1Twice[i][j]) < 0.1 * EPS ? 0.1 : Math.abs(resX1Twice[i][j]);
assertEquals(resX1Twice[i][j], resX1TwiceTran[i][j], EPS * ref);
}
}
{
final double val = func.differentiate(method.interpolate(x1Values, yValues[1]), x1Keys[2]).getData()[0];
final double resVal = interp.differentiateX1(x0Values, x1Values, yValues, x0Values[1], x1Keys[2]);
final double ref = Math.abs(val) < 0.1 * EPS ? 0.1 : Math.abs(val);
assertEquals(resVal, val, EPS * ref);
}
{
final double val = func.differentiate(method.interpolate(x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData()[1]), x0Keys[2]).getData()[0];
final double resVal = interp.differentiateX0(x0Values, x1Values, yValues, x0Keys[2], x1Values[1]);
final double ref = Math.abs(val) < 0.1 * EPS ? 0.1 : Math.abs(val);
assertEquals(resVal, val, EPS * ref);
}
{
final double val = func.differentiateTwice(method.interpolate(x1Values, yValues[1]), x1Keys[2]).getData()[0];
final double resVal = interp.differentiateX1Twice(x0Values, x1Values, yValues, x0Values[1], x1Keys[2]);
final double ref = Math.abs(val) < 0.1 * EPS ? 0.1 : Math.abs(val);
assertEquals(resVal, val, EPS * ref);
}
{
final double val = func.differentiateTwice(method.interpolate(x0Values, OG_ALGEBRA.getTranspose(new DoubleMatrix2D(yValues)).getData()[1]), x0Keys[2]).getData()[0];
final double resVal = interp.differentiateX0Twice(x0Values, x1Values, yValues, x0Keys[2], x1Values[1]);
final double ref = Math.abs(val) < 0.1 * EPS ? 0.1 : Math.abs(val);
assertEquals(resVal, val, EPS * ref);
}