    001    /*
    002     * Licensed to the Apache Software Foundation (ASF) under one or more
    003     * contributor license agreements.  See the NOTICE file distributed with
    004     * this work for additional information regarding copyright ownership.
    005     * The ASF licenses this file to You under the Apache License, Version 2.0
    006     * (the "License"); you may not use this file except in compliance with
    007     * the License.  You may obtain a copy of the License at
    008     *
    009     *      http://www.apache.org/licenses/LICENSE-2.0
    010     *
    011     * Unless required by applicable law or agreed to in writing, software
    012     * distributed under the License is distributed on an "AS IS" BASIS,
    013     * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
    014     * See the License for the specific language governing permissions and
    015     * limitations under the License.
    016     */
    017    package org.apache.commons.math3.analysis.polynomials;
    018
    019    import org.apache.commons.math3.analysis.UnivariateFunction;
    020    import org.apache.commons.math3.util.FastMath;
    021    import org.apache.commons.math3.util.MathArrays;
    022    import org.apache.commons.math3.exception.DimensionMismatchException;
    023    import org.apache.commons.math3.exception.NumberIsTooSmallException;
    024    import org.apache.commons.math3.exception.util.LocalizedFormats;
    025
    026    /**
    027     * Implements the representation of a real polynomial function in
    028     * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
    029     * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
    030     * Analysis</b>, ISBN 038795452X, chapter 2.
    031     * <p>
    032     * The approximated function should be smooth enough for Lagrange polynomial
    033     * to work well. Otherwise, consider using splines instead.</p>
    034     *
    035     * @version $Id: PolynomialFunctionLagrangeForm.java 1364387 2012-07-22 18:14:11Z tn $
    036     * @since 1.2
    037     */
    038    public class PolynomialFunctionLagrangeForm implements UnivariateFunction {
    039        /**
    040         * The coefficients of the polynomial, ordered by degree -- i.e.
    041         * coefficients[0] is the constant term and coefficients[n] is the
    042         * coefficient of x^n where n is the degree of the polynomial.
    043         */
    044        private double coefficients[];
    045        /**
    046         * Interpolating points (abscissas).
    047         */
    048        private final double x[];
    049        /**
    050         * Function values at interpolating points.
    051         */
    052        private final double y[];
    053        /**
    054         * Whether the polynomial coefficients are available.
    055         */
    056        private boolean coefficientsComputed;
    057
    058        /**
    059         * Construct a Lagrange polynomial with the given abscissas and function
    060         * values. The order of interpolating points are not important.
    061         * <p>
    062         * The constructor makes copy of the input arrays and assigns them.</p>
    063         *
    064         * @param x interpolating points
    065         * @param y function values at interpolating points
    066         * @throws DimensionMismatchException if the array lengths are different.
    067         * @throws NumberIsTooSmallException if the number of points is less than 2.
    068         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    069         * if two abscissae have the same value.
    070         */
    071        public PolynomialFunctionLagrangeForm(double x[], double y[]) {
    072            this.x = new double[x.length];
    073            this.y = new double[y.length];
    074            System.arraycopy(x, 0, this.x, 0, x.length);
    075            System.arraycopy(y, 0, this.y, 0, y.length);
    076            coefficientsComputed = false;
    077
    078            if (!verifyInterpolationArray(x, y, false)) {
    079                MathArrays.sortInPlace(this.x, this.y);
    080                // Second check in case some abscissa is duplicated.
    081                verifyInterpolationArray(this.x, this.y, true);
    082            }
    083        }
    084
    085        /**
    086         * Calculate the function value at the given point.
    087         *
    088         * @param z Point at which the function value is to be computed.
    089         * @return the function value.
    090         * @throws DimensionMismatchException if {@code x} and {@code y} have
    091         * different lengths.
    092         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    093         * if {@code x} is not sorted in strictly increasing order.
    094         * @throws NumberIsTooSmallException if the size of {@code x} is less
    095         * than 2.
    096         */
    097        public double value(double z) {
    098            return evaluateInternal(x, y, z);
    099        }
    100
    101        /**
    102         * Returns the degree of the polynomial.
    103         *
    104         * @return the degree of the polynomial
    105         */
    106        public int degree() {
    107            return x.length - 1;
    108        }
    109
    110        /**
    111         * Returns a copy of the interpolating points array.
    112         * <p>
    113         * Changes made to the returned copy will not affect the polynomial.</p>
    114         *
    115         * @return a fresh copy of the interpolating points array
    116         */
    117        public double[] getInterpolatingPoints() {
    118            double[] out = new double[x.length];
    119            System.arraycopy(x, 0, out, 0, x.length);
    120            return out;
    121        }
    122
    123        /**
    124         * Returns a copy of the interpolating values array.
    125         * <p>
    126         * Changes made to the returned copy will not affect the polynomial.</p>
    127         *
    128         * @return a fresh copy of the interpolating values array
    129         */
    130        public double[] getInterpolatingValues() {
    131            double[] out = new double[y.length];
    132            System.arraycopy(y, 0, out, 0, y.length);
    133            return out;
    134        }
    135
    136        /**
    137         * Returns a copy of the coefficients array.
    138         * <p>
    139         * Changes made to the returned copy will not affect the polynomial.</p>
    140         * <p>
    141         * Note that coefficients computation can be ill-conditioned. Use with caution
    142         * and only when it is necessary.</p>
    143         *
    144         * @return a fresh copy of the coefficients array
    145         */
    146        public double[] getCoefficients() {
    147            if (!coefficientsComputed) {
    148                computeCoefficients();
    149            }
    150            double[] out = new double[coefficients.length];
    151            System.arraycopy(coefficients, 0, out, 0, coefficients.length);
    152            return out;
    153        }
    154
    155        /**
    156         * Evaluate the Lagrange polynomial using
    157         * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
    158         * Neville's Algorithm</a>. It takes O(n^2) time.
    159         *
    160         * @param x Interpolating points array.
    161         * @param y Interpolating values array.
    162         * @param z Point at which the function value is to be computed.
    163         * @return the function value.
    164         * @throws DimensionMismatchException if {@code x} and {@code y} have
    165         * different lengths.
    166         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    167         * if {@code x} is not sorted in strictly increasing order.
    168         * @throws NumberIsTooSmallException if the size of {@code x} is less
    169         * than 2.
    170         */
    171        public static double evaluate(double x[], double y[], double z) {
    172            if (verifyInterpolationArray(x, y, false)) {
    173                return evaluateInternal(x, y, z);
    174            }
    175
    176            // Array is not sorted.
    177            final double[] xNew = new double[x.length];
    178            final double[] yNew = new double[y.length];
    179            System.arraycopy(x, 0, xNew, 0, x.length);
    180            System.arraycopy(y, 0, yNew, 0, y.length);
    181
    182            MathArrays.sortInPlace(xNew, yNew);
    183            // Second check in case some abscissa is duplicated.
    184            verifyInterpolationArray(xNew, yNew, true);
    185            return evaluateInternal(xNew, yNew, z);
    186        }
    187
    188        /**
    189         * Evaluate the Lagrange polynomial using
    190         * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
    191         * Neville's Algorithm</a>. It takes O(n^2) time.
    192         *
    193         * @param x Interpolating points array.
    194         * @param y Interpolating values array.
    195         * @param z Point at which the function value is to be computed.
    196         * @return the function value.
    197         * @throws DimensionMismatchException if {@code x} and {@code y} have
    198         * different lengths.
    199         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    200         * if {@code x} is not sorted in strictly increasing order.
    201         * @throws NumberIsTooSmallException if the size of {@code x} is less
    202         * than 2.
    203         */
    204        private static double evaluateInternal(double x[], double y[], double z) {
    205            int nearest = 0;
    206            final int n = x.length;
    207            final double[] c = new double[n];
    208            final double[] d = new double[n];
    209            double min_dist = Double.POSITIVE_INFINITY;
    210            for (int i = 0; i < n; i++) {
    211                // initialize the difference arrays
    212                c[i] = y[i];
    213                d[i] = y[i];
    214                // find out the abscissa closest to z
    215                final double dist = FastMath.abs(z - x[i]);
    216                if (dist < min_dist) {
    217                    nearest = i;
    218                    min_dist = dist;
    219                }
    220            }
    221
    222            // initial approximation to the function value at z
    223            double value = y[nearest];
    224
    225            for (int i = 1; i < n; i++) {
    226                for (int j = 0; j < n-i; j++) {
    227                    final double tc = x[j] - z;
    228                    final double td = x[i+j] - z;
    229                    final double divider = x[j] - x[i+j];
    230                    // update the difference arrays
    231                    final double w = (c[j+1] - d[j]) / divider;
    232                    c[j] = tc * w;
    233                    d[j] = td * w;
    234                }
    235                // sum up the difference terms to get the final value
    236                if (nearest < 0.5*(n-i+1)) {
    237                    value += c[nearest];    // fork down
    238                } else {
    239                    nearest--;
    240                    value += d[nearest];    // fork up
    241                }
    242            }
    243
    244            return value;
    245        }
    246
    247        /**
    248         * Calculate the coefficients of Lagrange polynomial from the
    249         * interpolation data. It takes O(n^2) time.
    250         * Note that this computation can be ill-conditioned: Use with caution
    251         * and only when it is necessary.
    252         */
    253        protected void computeCoefficients() {
    254            final int n = degree() + 1;
    255            coefficients = new double[n];
    256            for (int i = 0; i < n; i++) {
    257                coefficients[i] = 0.0;
    258            }
    259
    260            // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
    261            final double[] c = new double[n+1];
    262            c[0] = 1.0;
    263            for (int i = 0; i < n; i++) {
    264                for (int j = i; j > 0; j--) {
    265                    c[j] = c[j-1] - c[j] * x[i];
    266                }
    267                c[0] *= -x[i];
    268                c[i+1] = 1;
    269            }
    270
    271            final double[] tc = new double[n];
    272            for (int i = 0; i < n; i++) {
    273                // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
    274                double d = 1;
    275                for (int j = 0; j < n; j++) {
    276                    if (i != j) {
    277                        d *= x[i] - x[j];
    278                    }
    279                }
    280                final double t = y[i] / d;
    281                // Lagrange polynomial is the sum of n terms, each of which is a
    282                // polynomial of degree n-1. tc[] are the coefficients of the i-th
    283                // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
    284                tc[n-1] = c[n];     // actually c[n] = 1
    285                coefficients[n-1] += t * tc[n-1];
    286                for (int j = n-2; j >= 0; j--) {
    287                    tc[j] = c[j+1] + tc[j+1] * x[i];
    288                    coefficients[j] += t * tc[j];
    289                }
    290            }
    291
    292            coefficientsComputed = true;
    293        }
    294
    295        /**
    296         * Check that the interpolation arrays are valid.
    297         * The arrays features checked by this method are that both arrays have the
    298         * same length and this length is at least 2.
    299         *
    300         * @param x Interpolating points array.
    301         * @param y Interpolating values array.
    302         * @param abort Whether to throw an exception if {@code x} is not sorted.
    303         * @throws DimensionMismatchException if the array lengths are different.
    304         * @throws NumberIsTooSmallException if the number of points is less than 2.
    305         * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
    306         * if {@code x} is not sorted in strictly increasing order and {@code abort}
    307         * is {@code true}.
    308         * @return {@code false} if the {@code x} is not sorted in increasing order,
    309         * {@code true} otherwise.
    310         * @see #evaluate(double[], double[], double)
    311         * @see #computeCoefficients()
    312         */
    313        public static boolean verifyInterpolationArray(double x[], double y[], boolean abort) {
    314            if (x.length != y.length) {
    315                throw new DimensionMismatchException(x.length, y.length);
    316            }
    317            if (x.length < 2) {
    318                throw new NumberIsTooSmallException(LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
    319            }
    320
    321            return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort);
    322        }
    323    }