@Override protected double doFirstDerivative(double xValue) { int lowerIndex = lowerBoundIndex(xValue, xValues); int higherIndex = lowerIndex + 1; RealPolynomialFunction1D[] quadFirstDerivative = quadraticsFirstDerivative.get(); // at start of curve, or only one interval if (lowerIndex == 0 || intervalCount == 1) { RealPolynomialFunction1D quadraticFirstDerivative = quadFirstDerivative[0]; double x = xValue - xValues[1]; return quadraticFirstDerivative.applyAsDouble(x); } // at end of curve if (higherIndex >= intervalCount) { RealPolynomialFunction1D quadraticFirstDerivative = quadFirstDerivative[intervalCount - 2]; double x = xValue - xValues[intervalCount - 1]; return quadraticFirstDerivative.applyAsDouble(x); } RealPolynomialFunction1D quadratic1 = quadratics[lowerIndex - 1]; RealPolynomialFunction1D quadratic2 = quadratics[higherIndex - 1]; RealPolynomialFunction1D quadratic1FirstDerivative = quadFirstDerivative[lowerIndex - 1]; RealPolynomialFunction1D quadratic2FirstDerivative = quadFirstDerivative[higherIndex - 1]; double w = WEIGHT_FUNCTION.getWeight((xValues[higherIndex] - xValue) / (xValues[higherIndex] - xValues[lowerIndex])); return w * quadratic1FirstDerivative.applyAsDouble(xValue - xValues[lowerIndex]) + (1 - w) * quadratic2FirstDerivative.applyAsDouble(xValue - xValues[higherIndex]) + (quadratic2.applyAsDouble(xValue - xValues[higherIndex]) - quadratic1.applyAsDouble(xValue - xValues[lowerIndex])) / (xValues[higherIndex] - xValues[lowerIndex]); }
Bound(DoubleArray xValues, DoubleArray yValues) { super(xValues, yValues); this.xValues = xValues.toArrayUnsafe(); this.yValues = yValues.toArrayUnsafe(); this.intervalCount = xValues.size() - 1; this.quadratics = quadratics(this.xValues, this.yValues, this.intervalCount); this.quadraticsFirstDerivative = Suppliers.memoize(() -> quadraticsFirstDerivative(this.xValues, this.yValues, this.intervalCount)); }
@Override protected DoubleArray doParameterSensitivity(double xValue) { int lowerIndex = lowerBoundIndex(xValue, xValues); int higherIndex = lowerIndex + 1; int n = xValues.length; double[] temp = quadraticSensitivities(xValues, xValue, 1); result[0] = temp[0]; result[1] = temp[1]; double[] temp = quadraticSensitivities(xValues, xValue, n - 2); result[n - 3] = temp[0]; result[n - 2] = temp[1]; return DoubleArray.ofUnsafe(result); double[] temp1 = quadraticSensitivities(xValues, xValue, lowerIndex); double[] temp2 = quadraticSensitivities(xValues, xValue, higherIndex); double w = WEIGHT_FUNCTION.getWeight((xValues[higherIndex] - xValue) / (xValues[higherIndex] - xValues[lowerIndex])); result[lowerIndex - 1] = w * temp1[0];
@Override protected double doInterpolate(double xValue) { // x-value is less than the x-value of the last node (lowerIndex < intervalCount) int lowerIndex = lowerBoundIndex(xValue, xValues); int higherIndex = lowerIndex + 1; // at start of curve if (lowerIndex == 0) { RealPolynomialFunction1D quadratic = quadratics[0]; double x = xValue - xValues[1]; return quadratic.applyAsDouble(x); } // at end of curve if (higherIndex == intervalCount) { RealPolynomialFunction1D quadratic = quadratics[intervalCount - 2]; double x = xValue - xValues[intervalCount - 1]; return quadratic.applyAsDouble(x); } // normal case RealPolynomialFunction1D quadratic1 = quadratics[lowerIndex - 1]; RealPolynomialFunction1D quadratic2 = quadratics[higherIndex - 1]; double w = WEIGHT_FUNCTION.getWeight((xValues[higherIndex] - xValue) / (xValues[higherIndex] - xValues[lowerIndex])); return w * quadratic1.applyAsDouble(xValue - xValues[lowerIndex]) + (1 - w) * quadratic2.applyAsDouble(xValue - xValues[higherIndex]); }
@Override public BoundCurveInterpolator bind(DoubleArray xValues, DoubleArray yValues) { return new Bound(xValues, yValues); }
private static RealPolynomialFunction1D[] quadratics(double[] x, double[] y, int intervalCount) { if (intervalCount == 1) { double a = y[1]; double b = (y[1] - y[0]) / (x[1] - x[0]); return new RealPolynomialFunction1D[] {new RealPolynomialFunction1D(a, b)}; } RealPolynomialFunction1D[] quadratic = new RealPolynomialFunction1D[intervalCount - 1]; for (int i = 1; i < intervalCount; i++) { quadratic[i - 1] = quadratic(x, y, i); } return quadratic; }
private static RealPolynomialFunction1D[] quadraticsFirstDerivative(double[] x, double[] y, int intervalCount) { if (intervalCount == 1) { double b = (y[1] - y[0]) / (x[1] - x[0]); return new RealPolynomialFunction1D[] {new RealPolynomialFunction1D(b)}; } else { RealPolynomialFunction1D[] quadraticFirstDerivative = new RealPolynomialFunction1D[intervalCount - 1]; for (int i = 1; i < intervalCount; i++) { quadraticFirstDerivative[i - 1] = quadraticFirstDerivative(x, y, i); } return quadraticFirstDerivative; } }
@Override public BoundCurveInterpolator bind( BoundCurveExtrapolator extrapolatorLeft, BoundCurveExtrapolator extrapolatorRight) { return new Bound(this, extrapolatorLeft, extrapolatorRight); } }