HermitePolynomialApproximation.js

/*global define*/
define([
        './defaultValue',
        './defined',
        './DeveloperError',
        './Math'
    ], function(
        defaultValue,
        defined,
        DeveloperError,
        CesiumMath) {
    'use strict';

    var factorial = CesiumMath.factorial;

    function calculateCoefficientTerm(x, zIndices, xTable, derivOrder, termOrder, reservedIndices) {
        var result = 0;
        var reserved;
        var i;
        var j;

        if (derivOrder > 0) {
            for (i = 0; i < termOrder; i++) {
                reserved = false;
                for (j = 0; j < reservedIndices.length && !reserved; j++) {
                    if (i === reservedIndices[j]) {
                        reserved = true;
                    }
                }

                if (!reserved) {
                    reservedIndices.push(i);
                    result += calculateCoefficientTerm(x, zIndices, xTable, derivOrder - 1, termOrder, reservedIndices);
                    reservedIndices.splice(reservedIndices.length - 1, 1);
                }
            }

            return result;
        }

        result = 1;
        for (i = 0; i < termOrder; i++) {
            reserved = false;
            for (j = 0; j < reservedIndices.length && !reserved; j++) {
                if (i === reservedIndices[j]) {
                    reserved = true;
                }
            }

            if (!reserved) {
                result *= x - xTable[zIndices[i]];
            }
        }

        return result;
    }

    /**
     * An {@link InterpolationAlgorithm} for performing Hermite interpolation.
     *
     * @exports HermitePolynomialApproximation
     */
    var HermitePolynomialApproximation = {
        type : 'Hermite'
    };

    /**
     * Given the desired degree, returns the number of data points required for interpolation.
     *
     * @param {Number} degree The desired degree of interpolation.
     * @param {Number} [inputOrder=0]  The order of the inputs (0 means just the data, 1 means the data and its derivative, etc).
     * @returns {Number} The number of required data points needed for the desired degree of interpolation.
     * @exception {DeveloperError} degree must be 0 or greater.
     * @exception {DeveloperError} inputOrder must be 0 or greater.
     */
    HermitePolynomialApproximation.getRequiredDataPoints = function(degree, inputOrder) {
        inputOrder = defaultValue(inputOrder, 0);

        //>>includeStart('debug', pragmas.debug);
        if (!defined(degree)) {
            throw new DeveloperError('degree is required.');
        }
        if (degree < 0) {
            throw new DeveloperError('degree must be 0 or greater.');
        }
        if (inputOrder < 0) {
            throw new DeveloperError('inputOrder must be 0 or greater.');
        }
        //>>includeEnd('debug');

        return Math.max(Math.floor((degree + 1) / (inputOrder + 1)), 2);
    };

    /**
     * Interpolates values using Hermite Polynomial Approximation.
     *
     * @param {Number} x The independent variable for which the dependent variables will be interpolated.
     * @param {Number[]} xTable The array of independent variables to use to interpolate.  The values
     * in this array must be in increasing order and the same value must not occur twice in the array.
     * @param {Number[]} yTable The array of dependent variables to use to interpolate.  For a set of three
     * dependent values (p,q,w) at time 1 and time 2 this should be as follows: {p1, q1, w1, p2, q2, w2}.
     * @param {Number} yStride The number of dependent variable values in yTable corresponding to
     * each independent variable value in xTable.
     * @param {Number[]} [result] An existing array into which to store the result.
     * @returns {Number[]} The array of interpolated values, or the result parameter if one was provided.
     */
    HermitePolynomialApproximation.interpolateOrderZero = function(x, xTable, yTable, yStride, result) {
        if (!defined(result)) {
            result = new Array(yStride);
        }

        var i;
        var j;
        var d;
        var s;
        var len;
        var index;
        var length = xTable.length;
        var coefficients = new Array(yStride);

        for (i = 0; i < yStride; i++) {
            result[i] = 0;

            var l = new Array(length);
            coefficients[i] = l;
            for (j = 0; j < length; j++) {
                l[j] = [];
            }
        }

        var zIndicesLength = length, zIndices = new Array(zIndicesLength);

        for (i = 0; i < zIndicesLength; i++) {
            zIndices[i] = i;
        }

        var highestNonZeroCoef = length - 1;
        for (s = 0; s < yStride; s++) {
            for (j = 0; j < zIndicesLength; j++) {
                index = zIndices[j] * yStride + s;
                coefficients[s][0].push(yTable[index]);
            }

            for (i = 1; i < zIndicesLength; i++) {
                var nonZeroCoefficients = false;
                for (j = 0; j < zIndicesLength - i; j++) {
                    var zj = xTable[zIndices[j]];
                    var zn = xTable[zIndices[j + i]];

                    var numerator;
                    if (zn - zj <= 0) {
                        index = zIndices[j] * yStride + yStride * i + s;
                        numerator = yTable[index];
                        coefficients[s][i].push(numerator / factorial(i));
                    } else {
                        numerator = (coefficients[s][i - 1][j + 1] - coefficients[s][i - 1][j]);
                        coefficients[s][i].push(numerator / (zn - zj));
                    }
                    nonZeroCoefficients = nonZeroCoefficients || (numerator !== 0);
                }

                if (!nonZeroCoefficients) {
                    highestNonZeroCoef = i - 1;
                }
            }
        }

        for (d = 0, len = 0; d <= len; d++) {
            for (i = d; i <= highestNonZeroCoef; i++) {
                var tempTerm = calculateCoefficientTerm(x, zIndices, xTable, d, i, []);
                for (s = 0; s < yStride; s++) {
                    var coeff = coefficients[s][i][0];
                    result[s + d * yStride] += coeff * tempTerm;
                }
            }
        }

        return result;
    };

    var arrayScratch = [];

    /**
     * Interpolates values using Hermite Polynomial Approximation.
     *
     * @param {Number} x The independent variable for which the dependent variables will be interpolated.
     * @param {Number[]} xTable The array of independent variables to use to interpolate.  The values
     * in this array must be in increasing order and the same value must not occur twice in the array.
     * @param {Number[]} yTable The array of dependent variables to use to interpolate.  For a set of three
     * dependent values (p,q,w) at time 1 and time 2 this should be as follows: {p1, q1, w1, p2, q2, w2}.
     * @param {Number} yStride The number of dependent variable values in yTable corresponding to
     * each independent variable value in xTable.
     * @param {Number} inputOrder The number of derivatives supplied for input.
     * @param {Number} outputOrder The number of derivatives desired for output.
     * @param {Number[]} [result] An existing array into which to store the result.
     *
     * @returns {Number[]} The array of interpolated values, or the result parameter if one was provided.
     */
    HermitePolynomialApproximation.interpolate = function(x, xTable, yTable, yStride, inputOrder, outputOrder, result) {
        var resultLength = yStride * (outputOrder + 1);
        if (!defined(result)) {
            result = new Array(resultLength);
        }
        for (var r = 0; r < resultLength; r++) {
            result[r] = 0;
        }

        var length = xTable.length;
        // The zIndices array holds copies of the addresses of the xTable values
        // in the range we're looking at. Even though this just holds information already
        // available in xTable this is a much more convenient format.
        var zIndices = new Array(length * (inputOrder + 1));
        for (var i = 0; i < length; i++) {
            for (var j = 0; j < (inputOrder + 1); j++) {
                zIndices[i * (inputOrder + 1) + j] = i;
            }
        }

        var zIndiceslength = zIndices.length;
        var coefficients = arrayScratch;
        var highestNonZeroCoef = fillCoefficientList(coefficients, zIndices, xTable, yTable, yStride, inputOrder);
        var reservedIndices = [];

        var tmp = zIndiceslength * (zIndiceslength + 1) / 2;
        var loopStop = Math.min(highestNonZeroCoef, outputOrder);
        for (var d = 0; d <= loopStop; d++) {
            for (i = d; i <= highestNonZeroCoef; i++) {
                reservedIndices.length = 0;
                var tempTerm = calculateCoefficientTerm(x, zIndices, xTable, d, i, reservedIndices);
                var dimTwo = Math.floor(i * (1 - i) / 2) + (zIndiceslength * i);

                for (var s = 0; s < yStride; s++) {
                    var dimOne = Math.floor(s * tmp);
                    var coef = coefficients[dimOne + dimTwo];
                    result[s + d * yStride] += coef * tempTerm;
                }
            }
        }

        return result;
    };

    function fillCoefficientList(coefficients, zIndices, xTable, yTable, yStride, inputOrder) {
        var j;
        var index;
        var highestNonZero = -1;
        var zIndiceslength = zIndices.length;
        var tmp = zIndiceslength * (zIndiceslength + 1) / 2;

        for (var s = 0; s < yStride; s++) {
            var dimOne = Math.floor(s * tmp);

            for (j = 0; j < zIndiceslength; j++) {
                index = zIndices[j] * yStride * (inputOrder + 1) + s;
                coefficients[dimOne + j] = yTable[index];
            }

            for (var i = 1; i < zIndiceslength; i++) {
                var coefIndex = 0;
                var dimTwo = Math.floor(i * (1 - i) / 2) + (zIndiceslength * i);
                var nonZeroCoefficients = false;

                for (j = 0; j < zIndiceslength - i; j++) {
                    var zj = xTable[zIndices[j]];
                    var zn = xTable[zIndices[j + i]];

                    var numerator;
                    var coefficient;
                    if (zn - zj <= 0) {
                        index = zIndices[j] * yStride * (inputOrder + 1) + yStride * i + s;
                        numerator = yTable[index];
                        coefficient = (numerator / CesiumMath.factorial(i));
                        coefficients[dimOne + dimTwo + coefIndex] = coefficient;
                        coefIndex++;
                    } else {
                        var dimTwoMinusOne = Math.floor((i - 1) * (2 - i) / 2) + (zIndiceslength * (i - 1));
                        numerator = coefficients[dimOne + dimTwoMinusOne + j + 1] - coefficients[dimOne + dimTwoMinusOne + j];
                        coefficient = (numerator / (zn - zj));
                        coefficients[dimOne + dimTwo + coefIndex] = coefficient;
                        coefIndex++;
                    }
                    nonZeroCoefficients = nonZeroCoefficients || (numerator !== 0.0);
                }

                if (nonZeroCoefficients) {
                    highestNonZero = Math.max(highestNonZero, i);
                }
            }
        }

        return highestNonZero;
    }

    return HermitePolynomialApproximation;
});