* This implementation is based on the Akima implementation in the CubicSpline * class in the Math.NET Numerics library. The method referenced is * CubicSpline.InterpolateAkimaSorted *

*

* The {@link #interpolate(double[], double[]) interpolate} method returns a * {@link PolynomialSplineFunction} consisting of n cubic polynomials, defined * over the subintervals determined by the x values, {@code x[0] < x[i] ... < x[n]}. * The Akima algorithm requires that {@code n >= 5}. *

*/ export class AkimaSplineInterpolator implements UnivariateInterpolator { /** The minimum number of points that are needed to compute the function. */ private static readonly MINIMUM_NUMBER_POINTS = 5; /** * Computes an interpolating function for the data set. * * @param xvals the arguments for the interpolation points * @param yvals the values for the interpolation points * @return a function which interpolates the data set */ public interpolate(xvals: number[], yvals: number[]) : PolynomialSplineFunction { if (!xvals || !yvals) { throw new Error("Null argument."); } if (xvals.length != yvals.length) { throw new Error("Dimension mismatch for xvals and yvals."); } if (xvals.length < AkimaSplineInterpolator.MINIMUM_NUMBER_POINTS) { throw new Error("Number of points is too small."); } MathArrays.checkOrder(xvals); let numberOfDiffAndWeightElements = xvals.length - 1; let differences: number[] = Array(numberOfDiffAndWeightElements); let weights: number[] = Array(numberOfDiffAndWeightElements); for (let i = 0; i < differences.length; i++) { differences[i] = (yvals[i + 1] - yvals[i]) / (xvals[i + 1] - xvals[i]); } for (let i = 1; i < weights.length; i++) { weights[i] = Math.abs(differences[i] - differences[i - 1]); } // Prepare Hermite interpolation scheme. let firstDerivatives: number[] = Array(xvals.length); for (let i = 2; i < firstDerivatives.length - 2; i++) { let wP = weights[i + 1]; let wM = weights[i - 1]; if (Math.abs(wP) < EPSILON && Math.abs(wM) < EPSILON) { let xv = xvals[i]; let xvP = xvals[i + 1]; let xvM = xvals[i - 1]; firstDerivatives[i] = (((xvP - xv) * differences[i - 1]) + ((xv - xvM) * differences[i])) / (xvP - xvM); } else { firstDerivatives[i] = ((wP * differences[i - 1]) + (wM * differences[i])) / (wP + wM); } } firstDerivatives[0] = this.differentiateThreePoint(xvals, yvals, 0, 0, 1, 2); firstDerivatives[1] = this.differentiateThreePoint(xvals, yvals, 1, 0, 1, 2); firstDerivatives[xvals.length - 2] = this.differentiateThreePoint(xvals, yvals, xvals.length - 2, xvals.length - 3, xvals.length - 2, xvals.length - 1); firstDerivatives[xvals.length - 1] = this.differentiateThreePoint(xvals, yvals, xvals.length - 1, xvals.length - 3, xvals.length - 2, xvals.length - 1); return this.interpolateHermiteSorted(xvals, yvals, firstDerivatives); } /** * Three point differentiation helper, modeled off of the same method in the * Math.NET CubicSpline class. This is used by both the Apache Math and the * Math.NET Akima Cubic Spline algorithms * * @param xvals x values to calculate the numerical derivative with * @param yvals y values to calculate the numerical derivative with * @param indexOfDifferentiation index of the elemnt we are calculating the derivative around * @param indexOfFirstSample index of the first element to sample for the three point method * @param indexOfSecondsample index of the second element to sample for the three point method * @param indexOfThirdSample index of the third element to sample for the three point method * @return the derivative */ private differentiateThreePoint(xvals: number[], yvals: number[], indexOfDifferentiation: number, indexOfFirstSample: number, indexOfSecondsample: number, indexOfThirdSample: number) : number { let x0 = yvals[indexOfFirstSample]; let x1 = yvals[indexOfSecondsample]; let x2 = yvals[indexOfThirdSample]; let t = xvals[indexOfDifferentiation] - xvals[indexOfFirstSample]; let t1 = xvals[indexOfSecondsample] - xvals[indexOfFirstSample]; let t2 = xvals[indexOfThirdSample] - xvals[indexOfFirstSample]; let a = (x2 - x0 - (t2 / t1 * (x1 - x0))) / (t2 * t2 - t1 * t2); let b = (x1 - x0 - a * t1 * t1) / t1; return (2 * a * t) + b; } /** * Creates a Hermite cubic spline interpolation from the set of (x,y) value * pairs and their derivatives. This is modeled off of the * InterpolateHermiteSorted method in the Math.NET CubicSpline class. * * @param xvals x values for interpolation * @param yvals y values for interpolation * @param firstDerivatives first derivative values of the function * @return polynomial that fits the function */ private interpolateHermiteSorted(xvals: number[], yvals: number[], firstDerivatives: number[]) : PolynomialSplineFunction { if (xvals.length != yvals.length) { throw new Error("Dimension mismatch"); } if (xvals.length != firstDerivatives.length) { throw new Error("Dimension mismatch"); } let minimumLength = 2; if (xvals.length < minimumLength) { throw new Error("Not enough points."); } let size = xvals.length - 1; let polynomials: PolynomialFunction[] = Array(size); let coefficients: number[] = Array(4); for (let i = 0; i < polynomials.length; i++) { let w = xvals[i + 1] - xvals[i]; let w2 = w * w; let yv = yvals[i]; let yvP = yvals[i + 1]; let fd = firstDerivatives[i]; let fdP = firstDerivatives[i + 1]; coefficients[0] = yv; coefficients[1] = firstDerivatives[i]; coefficients[2] = (3 * (yvP - yv) / w - 2 * fd - fdP) / w; coefficients[3] = (2 * (yv - yvP) / w + fd + fdP) / w2; polynomials[i] = new PolynomialFunction(coefficients); } return new PolynomialSplineFunction(xvals, polynomials); } } // end class AkimaSplineInterpolator //------------------------------------------------------------------------------ /** * Computes a natural (also known as "free", "unclamped") cubic spline interpolation for the data set. *

* The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction} * consisting of n cubic polynomials, defined over the subintervals determined by the x values, * {@code x[0] < x[i] ... < x[n].} The x values are referred to as "knot points."

*

* The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest * knot point and strictly less than the largest knot point is computed by finding the subinterval to which * x belongs and computing the value of the corresponding polynomial at `x - x[i] ` where * `i` is the index of the subinterval. See {@link PolynomialSplineFunction} for more details. *

*

* The interpolating polynomials satisfy:

*
1. The value of the PolynomialSplineFunction at each of the input x values equals the * corresponding y value.
2. *
3. Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials * "match up" at the knot points, as do their first and second derivatives).
4. *

*

* The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires, * Numerical Analysis, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131. *

* */ export class SplineInterpolator implements UnivariateInterpolator { /** * Computes an interpolating function for the data set. * @param x the arguments for the interpolation points * @param y the values for the interpolation points * @return a function which interpolates the data set */ public interpolate(x: number[], y: number[]) : PolynomialSplineFunction { if (x.length != y.length) { throw new Error("Dimension msmatch."); } if (x.length < 3) { throw new Error("Number of points is too small."); } // Number of intervals. The number of data points is n + 1. let n = x.length - 1; MathArrays.checkOrder(x); // Differences between knot points let h: number[] = Array(n); for (let i = 0; i < n; i++) { h[i] = x[i + 1] - x[i]; } let mu: number[] = Array(n); let z: number[] = Array(n + 1); mu[0] = 0; z[0] = 0; let g = 0; for (let i = 1; i < n; i++) { g = 2 * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1]; mu[i] = h[i] / g; z[i] = (3 * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) / (h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g; } // cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants) let b: number[] = Array(n); let c: number[] = Array(n + 1); let d: number[] = Array(n); z[n] = 0; c[n] = 0; for (let j = n -1; j >=0; j--) { c[j] = z[j] - mu[j] * c[j + 1]; b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2 * c[j]) / 3; d[j] = (c[j + 1] - c[j]) / (3 * h[j]); } let polynomials: PolynomialFunction[] = Array(n); let coefficients: number[] = Array(4); for (let i = 0; i < n; i++) { coefficients[0] = y[i]; coefficients[1] = b[i]; coefficients[2] = c[i]; coefficients[3] = d[i]; polynomials[i] = new PolynomialFunction(coefficients); } return new PolynomialSplineFunction(x, polynomials); } } // end class SplineInterpolator //------------------------------------------------------------------------------ /** * Implements a linear function for interpolation of real univariate functions. */ export class LinearInterpolator implements UnivariateInterpolator { /** * Computes a linear interpolating function for the data set. * * @param x the arguments for the interpolation points * @param y the values for the interpolation points * @return a function which interpolates the data set */ public interpolate(x: number[], y: number[]) : PolynomialSplineFunction { if (x.length != y.length) { throw new Error("Dimension msmatch."); } if (x.length < 2) { throw new Error("Number of points is too small."); } // Number of intervals. The number of data points is n + 1. let n = x.length - 1; MathArrays.checkOrder(x); // Slope of the lines between the datapoints. let m: number[] = Array(n); for (let i = 0; i < n; i++) { m[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]); } let polynomials: PolynomialFunction[] = Array(n); let coefficients: number[] = Array(2); for (let i = 0; i < n; i++) { coefficients[0] = y[i]; coefficients[1] = m[i]; polynomials[i] = new PolynomialFunction(coefficients); } return new PolynomialSplineFunction(x, polynomials); } } // end class LinearInterpolator //------------------------------------------------------------------------------ /** * Represents a polynomial spline function. *

* A polynomial spline function consists of a set of * interpolating polynomials and an ascending array of domain * knot points, determining the intervals over which the spline function * is defined by the constituent polynomials. The polynomials are assumed to * have been computed to match the values of another function at the knot * points. The value consistency constraints are not currently enforced by * `PolynomialSplineFunction` itself, but are assumed to hold among * the polynomials and knot points passed to the constructor.

*

* N.B.: The polynomials in the `polynomials` property must be * centered on the knot points to compute the spline function values. * See below.

*

* The domain of the polynomial spline function is * `[smallest knot, largest knot]`. Attempts to evaluate the * function at values outside of this range generate IllegalArgumentExceptions. *

*

* The value of the polynomial spline function for an argument `x` * is computed as follows: *

*
1. The knot array is searched to find the segment to which `x` * belongs. If `x` is less than the smallest knot point or greater * than the largest one, an `IllegalArgumentException` * is thrown.
2. *
3. Let `j` be the index of the largest knot point that is less * than or equal to `x`. The value returned is * {@code polynomials[j](x - knot[j])}
*/ class PolynomialSplineFunction { /** * Spline segment interval delimiters (knots). * Size is n + 1 for n segments. */ private knots: number[]; /** * The polynomial functions that make up the spline. The first element * determines the value of the spline over the first subinterval, the * second over the second, etc. Spline function values are determined by * evaluating these functions at {@code (x - knot[i])} where i is the * knot segment to which x belongs. */ private polynomials: PolynomialFunction[]; /** * Number of spline segments. It is equal to the number of polynomials and * to the number of partition points - 1. */ private n: number; /** * Construct a polynomial spline function with the given segment delimiters * and interpolating polynomials. * The constructor copies both arrays and assigns the copies to the knots * and polynomials properties, respectively. * * @param knots Spline segment interval delimiters. * @param polynomials Polynomial functions that make up the spline. */ public constructor(knots: number[], polynomials: PolynomialFunction[]) { if (!knots || !polynomials) { throw new Error("Null argument."); } if (knots.length < 2) { throw new Error("Illegal argument knots."); } if (knots.length - 1 != polynomials.length) { throw new Error("Dimension mismatch."); } MathArrays.checkOrder(knots); this.n = knots.length -1; this.knots = knots.slice(); this.polynomials = polynomials.slice(); } /** * Compute the value for the function. * See {@link PolynomialSplineFunction} for details on the algorithm for * computing the value of the function. * * @param v Point for which the function value should be computed. * @return the value. */ public value(v: number) : number { // Patched to allow an extended function domain. // if (v < this.knots[0] || v > this.knots[this.n]) { // throw new Error("Argument value out of range."); // } let i = Arrays.binarySearch(this.knots, v); if (i < 0) { i = -i - 2; } // This will handle the case where v is the last knot value // There are only n-1 polynomials, so if v is the last knot // then we will use the last polynomial to calculate the value. // Patched to allow an extended function domain. // if ( i >= this.polynomials.length ) { // i--; // } i = Math.max(0, Math.min(i, this.polynomials.length - 1)); return this.polynomials[i].value(v - this.knots[i]); } } // end class PolynomialSplineFunction //------------------------------------------------------------------------------ /** * Immutable representation of a real polynomial function with real coefficients. *

* Horner's Method * is used to evaluate the function.

*/ class PolynomialFunction { /** * The coefficients of the polynomial, ordered by degree -- i.e., * coefficients[0] is the constant term and coefficients[n] is the * coefficient of x^n where n is the degree of the polynomial. */ private coefficients: number[]; /** * Construct a polynomial with the given coefficients. The first element * of the coefficients array is the constant term. Higher degree * coefficients follow in sequence. The degree of the resulting polynomial * is the index of the last non-null element of the array, or 0 if all elements * are null. *

* The constructor makes a copy of the input array and assigns the copy to * the coefficients property.

* * @param c Polynomial coefficients. */ constructor(c: number[]) { let n = c.length; if (n == 0) { throw new Error("Empty polynomials coefficients array"); } while ((n > 1) && (c[n - 1] == 0)) { --n; } this.coefficients = c.slice(); } /** * Compute the value of the function for the given argument. *

* The value returned is

* {@code coefficients[n] * x^n + ... + coefficients[1] * x + coefficients[0]} *

* * @param x Argument for which the function value should be computed. * @return the value of the polynomial at the given point. * * @see org.hipparchus.analysis.UnivariateFunction#value(double) */ public value(x: number) : number { return PolynomialFunction.evaluate(this.coefficients, x); } /** * Uses Horner's Method to evaluate the polynomial with the given coefficients at * the argument. * * @param coefficients Coefficients of the polynomial to evaluate. * @param argument Input value. * @return the value of the polynomial. */ protected static evaluate(coefficients: number[], argument: number) : number { let n = coefficients.length; if (n == 0) { throw new Error("Empty polynomials coefficients array."); } let result = coefficients[n - 1]; for (let j = n - 2; j >= 0; j--) { result = argument * result + coefficients[j]; } return result; } } // end class PolynomialFunction //------------------------------------------------------------------------------ /** * Arrays utilities. */ class MathArrays { /** * Check that the given array is sorted in strictly increasing order. * * @param val Values. */ public static checkOrder(val: number[]) { let previous = val[0]; let max = val.length; for (let index = 1; index < max; index++) { if (val[index] <= previous) { throw new Error("Non-monotonic sequence exception."); } previous = val[index]; } } } // end class MathArrays //------------------------------------------------------------------------------ // Corresponds to java.util.Arrays. class Arrays { public static binarySearch(a: number[], key: number) : number { let low = 0; let high = a.length - 1; while (low <= high) { let mid = (low + high) >>> 1; let midVal = a[mid]; if (midVal < key) low = mid + 1; else if (midVal > key) high = mid - 1; else if (midVal == key) return mid; else { // values might be NaN throw new Error("Invalid number encountered in binary search."); } } return -(low + 1); // key not found. } } // end class Arrays } // end namespace ApacheCommonsMath