bValues = bValuesCalculator(slopes, first);
intervalsA = getIntervalsA(intervals, slopes, first, bValues);
intervalsB = getIntervalsB(intervals, slopes, first, aValues);
}
final double[] second = secondDerivativeCalculator(initialSecond, intervalsA, intervalsB);
firstWithSensitivity[0] = new DoubleMatrix1D(first);
secondWithSensitivity[0] = new DoubleMatrix1D(second);
/*
* Centered finite difference method is used for computing node sensitivity
*/
int nExtra = (nDataPts == yValuesLen) ? 0 : 1;
final double[] yValuesUp = Arrays.copyOf(yValues, nDataPts + 2 * nExtra);
final double[] yValuesDw = Arrays.copyOf(yValues, nDataPts + 2 * nExtra);
final double[][] tmpFirst = new double[nDataPts][nDataPts];
final double[][] tmpSecond = new double[nDataPts][nDataPts];
for (int l = nExtra; l < nDataPts + nExtra; ++l) {
final double den = Math.abs(yValues[l]) < SMALL ? EPS : yValues[l] * EPS;
yValuesUp[l] = Math.abs(yValues[l]) < SMALL ? EPS : yValues[l] * (1. + EPS);
yValuesDw[l] = Math.abs(yValues[l]) < SMALL ? -EPS : yValues[l] * (1. - EPS);
final double[] yValuesSrtUp = Arrays.copyOfRange(yValuesUp, nExtra, nDataPts + nExtra);
final double[] yValuesSrtDw = Arrays.copyOfRange(yValuesDw, nExtra, nDataPts + nExtra);
final DoubleMatrix1D[] yValuesUpDw = new DoubleMatrix1D[] {new DoubleMatrix1D(yValuesUp), new DoubleMatrix1D(yValuesDw) };
final DoubleMatrix1D[] yValuesSrtUpDw = new DoubleMatrix1D[] {new DoubleMatrix1D(yValuesSrtUp), new DoubleMatrix1D(yValuesSrtDw) };
final DoubleMatrix1D[] firstSecondUpDw = new DoubleMatrix1D[4];
for (int ii = 0; ii < 2; ++ii) {
final double[] slopesUpDw = _solver.slopesCalculator(yValuesSrtUpDw[ii].getData(), intervals);
final PiecewisePolynomialResult resultUpDw = _method.interpolate(xValues, yValuesUpDw[ii].getData());
final double[] initialFirstUpDw = _function.differentiate(resultUpDw, xValues).getData()[0];
final double[] initialSecondUpDw = _function.differentiateTwice(resultUpDw, xValues).getData()[0];
double[] firstUpDw = firstDerivativeCalculator(yValuesSrtUpDw[ii].getData(), intervals, slopesUpDw, initialFirstUpDw);
boolean modFirstUpDw = false;
double[] aValuesUpDw = aValuesCalculator(slopesUpDw, firstUpDw);
double[] bValuesUpDw = bValuesCalculator(slopesUpDw, firstUpDw);
double[][] intervalsAUpDw = getIntervalsA(intervals, slopesUpDw, firstUpDw, bValuesUpDw);
double[][] intervalsBUpDw = getIntervalsB(intervals, slopesUpDw, firstUpDw, aValuesUpDw);
while (modFirstUpDw == false) {
k = 0;
for (int i = 0; i < nDataPts - 2; ++i) {
if (firstUpDw[i + 1] > 0.) {
if (intervalsAUpDw[i + 1][1] + Math.abs(intervalsAUpDw[i + 1][1]) * ERROR < intervalsBUpDw[i][0] - Math.abs(intervalsBUpDw[i][0]) * ERROR |
intervalsAUpDw[i + 1][0] - Math.abs(intervalsAUpDw[i + 1][0]) * ERROR > intervalsBUpDw[i][1] + Math.abs(intervalsBUpDw[i][1]) * ERROR) {
++k;
firstUpDw[i + 1] = firstDerivativesRecalculator(intervals, slopesUpDw, aValuesUpDw, bValuesUpDw, i + 1);
}
}
}
if (k == 0) {
modFirstUpDw = true;
}
aValuesUpDw = aValuesCalculator(slopesUpDw, firstUpDw);
bValuesUpDw = bValuesCalculator(slopesUpDw, firstUpDw);
intervalsAUpDw = getIntervalsA(intervals, slopesUpDw, firstUpDw, bValuesUpDw);
intervalsBUpDw = getIntervalsB(intervals, slopesUpDw, firstUpDw, aValuesUpDw);
}
final double[] secondUpDw = secondDerivativeCalculator(initialSecondUpDw, intervalsAUpDw, intervalsBUpDw);
firstSecondUpDw[ii] = new DoubleMatrix1D(firstUpDw);
firstSecondUpDw[2 + ii] = new DoubleMatrix1D(secondUpDw);
}
for (int j = 0; j < nDataPts; ++j) {
tmpFirst[j][l - nExtra] = 0.5 * (firstSecondUpDw[0].getData()[j] - firstSecondUpDw[1].getData()[j]) / den;
tmpSecond[j][l - nExtra] = 0.5 * (firstSecondUpDw[2].getData()[j] - firstSecondUpDw[3].getData()[j]) / den;
}
yValuesUp[l] = yValues[l];
yValuesDw[l] = yValues[l];
}
for (int i = 0; i < nDataPts; ++i) {
firstWithSensitivity[i + 1] = new DoubleMatrix1D(tmpFirst[i]);
secondWithSensitivity[i + 1] = new DoubleMatrix1D(tmpSecond[i]);
}
final DoubleMatrix2D[] resMatrix = _solver.solveWithSensitivity(yValuesSrt, intervals, slopes, slopesSensitivity, firstWithSensitivity, secondWithSensitivity);
for (int l = 0; l < nDataPts; ++l) {
DoubleMatrix2D m = resMatrix[l];
final int rows = m.getNumberOfRows();
final int cols = m.getNumberOfColumns();
for (int i = 0; i < rows; ++i) {
for (int j = 0; j < cols; ++j) {
ArgumentChecker.isTrue(Doubles.isFinite(m.getEntry(i, j)), "Matrix contains a NaN or infinite");
}
}
}
final DoubleMatrix2D coefMatrix = resMatrix[0];
final DoubleMatrix2D[] coefSenseMatrix = new DoubleMatrix2D[nDataPts - 1];
System.arraycopy(resMatrix, 1, coefSenseMatrix, 0, nDataPts - 1);
final int nCoefs = coefMatrix.getNumberOfColumns();
return new PiecewisePolynomialResultsWithSensitivity(new DoubleMatrix1D(xValues), coefMatrix, nCoefs, 1, coefSenseMatrix);
}