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    }