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 */
017package org.apache.commons.math4.analysis.polynomials;
018
019import org.apache.commons.math4.analysis.UnivariateFunction;
020import org.apache.commons.math4.exception.DimensionMismatchException;
021import org.apache.commons.math4.exception.NonMonotonicSequenceException;
022import org.apache.commons.math4.exception.NumberIsTooSmallException;
023import org.apache.commons.math4.exception.util.LocalizedFormats;
024import org.apache.commons.math4.util.FastMath;
025import org.apache.commons.math4.util.MathArrays;
026
027/**
028 * Implements the representation of a real polynomial function in
029 * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
030 * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
031 * Analysis</b>, ISBN 038795452X, chapter 2.
032 * <p>
033 * The approximated function should be smooth enough for Lagrange polynomial
034 * to work well. Otherwise, consider using splines instead.</p>
035 *
036 * @since 1.2
037 */
038public 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 NonMonotonicSequenceException
069     * if two abscissae have the same value.
070     */
071    public PolynomialFunctionLagrangeForm(double x[], double y[])
072        throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException {
073        this.x = new double[x.length];
074        this.y = new double[y.length];
075        System.arraycopy(x, 0, this.x, 0, x.length);
076        System.arraycopy(y, 0, this.y, 0, y.length);
077        coefficientsComputed = false;
078
079        if (!verifyInterpolationArray(x, y, false)) {
080            MathArrays.sortInPlace(this.x, this.y);
081            // Second check in case some abscissa is duplicated.
082            verifyInterpolationArray(this.x, this.y, true);
083        }
084    }
085
086    /**
087     * Calculate the function value at the given point.
088     *
089     * @param z Point at which the function value is to be computed.
090     * @return the function value.
091     * @throws DimensionMismatchException if {@code x} and {@code y} have
092     * different lengths.
093     * @throws org.apache.commons.math4.exception.NonMonotonicSequenceException
094     * if {@code x} is not sorted in strictly increasing order.
095     * @throws NumberIsTooSmallException if the size of {@code x} is less
096     * than 2.
097     */
098    @Override
099    public double value(double z) {
100        return evaluateInternal(x, y, z);
101    }
102
103    /**
104     * Returns the degree of the polynomial.
105     *
106     * @return the degree of the polynomial
107     */
108    public int degree() {
109        return x.length - 1;
110    }
111
112    /**
113     * Returns a copy of the interpolating points array.
114     * <p>
115     * Changes made to the returned copy will not affect the polynomial.</p>
116     *
117     * @return a fresh copy of the interpolating points array
118     */
119    public double[] getInterpolatingPoints() {
120        double[] out = new double[x.length];
121        System.arraycopy(x, 0, out, 0, x.length);
122        return out;
123    }
124
125    /**
126     * Returns a copy of the interpolating values array.
127     * <p>
128     * Changes made to the returned copy will not affect the polynomial.</p>
129     *
130     * @return a fresh copy of the interpolating values array
131     */
132    public double[] getInterpolatingValues() {
133        double[] out = new double[y.length];
134        System.arraycopy(y, 0, out, 0, y.length);
135        return out;
136    }
137
138    /**
139     * Returns a copy of the coefficients array.
140     * <p>
141     * Changes made to the returned copy will not affect the polynomial.</p>
142     * <p>
143     * Note that coefficients computation can be ill-conditioned. Use with caution
144     * and only when it is necessary.</p>
145     *
146     * @return a fresh copy of the coefficients array
147     */
148    public double[] getCoefficients() {
149        if (!coefficientsComputed) {
150            computeCoefficients();
151        }
152        double[] out = new double[coefficients.length];
153        System.arraycopy(coefficients, 0, out, 0, coefficients.length);
154        return out;
155    }
156
157    /**
158     * Evaluate the Lagrange polynomial using
159     * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
160     * Neville's Algorithm</a>. It takes O(n^2) time.
161     *
162     * @param x Interpolating points array.
163     * @param y Interpolating values array.
164     * @param z Point at which the function value is to be computed.
165     * @return the function value.
166     * @throws DimensionMismatchException if {@code x} and {@code y} have
167     * different lengths.
168     * @throws NonMonotonicSequenceException
169     * if {@code x} is not sorted in strictly increasing order.
170     * @throws NumberIsTooSmallException if the size of {@code x} is less
171     * than 2.
172     */
173    public static double evaluate(double x[], double y[], double z)
174        throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException {
175        if (verifyInterpolationArray(x, y, false)) {
176            return evaluateInternal(x, y, z);
177        }
178
179        // Array is not sorted.
180        final double[] xNew = new double[x.length];
181        final double[] yNew = new double[y.length];
182        System.arraycopy(x, 0, xNew, 0, x.length);
183        System.arraycopy(y, 0, yNew, 0, y.length);
184
185        MathArrays.sortInPlace(xNew, yNew);
186        // Second check in case some abscissa is duplicated.
187        verifyInterpolationArray(xNew, yNew, true);
188        return evaluateInternal(xNew, yNew, z);
189    }
190
191    /**
192     * Evaluate the Lagrange polynomial using
193     * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
194     * Neville's Algorithm</a>. It takes O(n^2) time.
195     *
196     * @param x Interpolating points array.
197     * @param y Interpolating values array.
198     * @param z Point at which the function value is to be computed.
199     * @return the function value.
200     * @throws DimensionMismatchException if {@code x} and {@code y} have
201     * different lengths.
202     * @throws org.apache.commons.math4.exception.NonMonotonicSequenceException
203     * if {@code x} is not sorted in strictly increasing order.
204     * @throws NumberIsTooSmallException if the size of {@code x} is less
205     * than 2.
206     */
207    private static double evaluateInternal(double x[], double y[], double z) {
208        int nearest = 0;
209        final int n = x.length;
210        final double[] c = new double[n];
211        final double[] d = new double[n];
212        double min_dist = Double.POSITIVE_INFINITY;
213        for (int i = 0; i < n; i++) {
214            // initialize the difference arrays
215            c[i] = y[i];
216            d[i] = y[i];
217            // find out the abscissa closest to z
218            final double dist = FastMath.abs(z - x[i]);
219            if (dist < min_dist) {
220                nearest = i;
221                min_dist = dist;
222            }
223        }
224
225        // initial approximation to the function value at z
226        double value = y[nearest];
227
228        for (int i = 1; i < n; i++) {
229            for (int j = 0; j < n-i; j++) {
230                final double tc = x[j] - z;
231                final double td = x[i+j] - z;
232                final double divider = x[j] - x[i+j];
233                // update the difference arrays
234                final double w = (c[j+1] - d[j]) / divider;
235                c[j] = tc * w;
236                d[j] = td * w;
237            }
238            // sum up the difference terms to get the final value
239            if (nearest < 0.5*(n-i+1)) {
240                value += c[nearest];    // fork down
241            } else {
242                nearest--;
243                value += d[nearest];    // fork up
244            }
245        }
246
247        return value;
248    }
249
250    /**
251     * Calculate the coefficients of Lagrange polynomial from the
252     * interpolation data. It takes O(n^2) time.
253     * Note that this computation can be ill-conditioned: Use with caution
254     * and only when it is necessary.
255     */
256    protected void computeCoefficients() {
257        final int n = degree() + 1;
258        coefficients = new double[n];
259        for (int i = 0; i < n; i++) {
260            coefficients[i] = 0.0;
261        }
262
263        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
264        final double[] c = new double[n+1];
265        c[0] = 1.0;
266        for (int i = 0; i < n; i++) {
267            for (int j = i; j > 0; j--) {
268                c[j] = c[j-1] - c[j] * x[i];
269            }
270            c[0] *= -x[i];
271            c[i+1] = 1;
272        }
273
274        final double[] tc = new double[n];
275        for (int i = 0; i < n; i++) {
276            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
277            double d = 1;
278            for (int j = 0; j < n; j++) {
279                if (i != j) {
280                    d *= x[i] - x[j];
281                }
282            }
283            final double t = y[i] / d;
284            // Lagrange polynomial is the sum of n terms, each of which is a
285            // polynomial of degree n-1. tc[] are the coefficients of the i-th
286            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
287            tc[n-1] = c[n];     // actually c[n] = 1
288            coefficients[n-1] += t * tc[n-1];
289            for (int j = n-2; j >= 0; j--) {
290                tc[j] = c[j+1] + tc[j+1] * x[i];
291                coefficients[j] += t * tc[j];
292            }
293        }
294
295        coefficientsComputed = true;
296    }
297
298    /**
299     * Check that the interpolation arrays are valid.
300     * The arrays features checked by this method are that both arrays have the
301     * same length and this length is at least 2.
302     *
303     * @param x Interpolating points array.
304     * @param y Interpolating values array.
305     * @param abort Whether to throw an exception if {@code x} is not sorted.
306     * @throws DimensionMismatchException if the array lengths are different.
307     * @throws NumberIsTooSmallException if the number of points is less than 2.
308     * @throws org.apache.commons.math4.exception.NonMonotonicSequenceException
309     * if {@code x} is not sorted in strictly increasing order and {@code abort}
310     * is {@code true}.
311     * @return {@code false} if the {@code x} is not sorted in increasing order,
312     * {@code true} otherwise.
313     * @see #evaluate(double[], double[], double)
314     * @see #computeCoefficients()
315     */
316    public static boolean verifyInterpolationArray(double x[], double y[], boolean abort)
317        throws DimensionMismatchException, NumberIsTooSmallException, NonMonotonicSequenceException {
318        if (x.length != y.length) {
319            throw new DimensionMismatchException(x.length, y.length);
320        }
321        if (x.length < 2) {
322            throw new NumberIsTooSmallException(LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length, true);
323        }
324
325        return MathArrays.checkOrder(x, MathArrays.OrderDirection.INCREASING, true, abort);
326    }
327}