1/*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements.  See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License.  You may obtain a copy of the License at
8 *
9 *      http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17package org.apache.commons.math.analysis.polynomials;
18
19import org.apache.commons.math.DuplicateSampleAbscissaException;
20import org.apache.commons.math.MathRuntimeException;
21import org.apache.commons.math.analysis.UnivariateRealFunction;
22import org.apache.commons.math.FunctionEvaluationException;
23import org.apache.commons.math.exception.util.LocalizedFormats;
24import org.apache.commons.math.util.FastMath;
25
26/**
27 * Implements the representation of a real polynomial function in
28 * <a href="http://mathworld.wolfram.com/LagrangeInterpolatingPolynomial.html">
29 * Lagrange Form</a>. For reference, see <b>Introduction to Numerical
30 * Analysis</b>, ISBN 038795452X, chapter 2.
31 * <p>
32 * The approximated function should be smooth enough for Lagrange polynomial
33 * to work well. Otherwise, consider using splines instead.</p>
34 *
35 * @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 févr. 2011) $
36 * @since 1.2
37 */
38public class PolynomialFunctionLagrangeForm implements UnivariateRealFunction {
39
40    /**
41     * The coefficients of the polynomial, ordered by degree -- i.e.
42     * coefficients[0] is the constant term and coefficients[n] is the
43     * coefficient of x^n where n is the degree of the polynomial.
44     */
45    private double coefficients[];
46
47    /**
48     * Interpolating points (abscissas).
49     */
50    private final double x[];
51
52    /**
53     * Function values at interpolating points.
54     */
55    private final double y[];
56
57    /**
58     * Whether the polynomial coefficients are available.
59     */
60    private boolean coefficientsComputed;
61
62    /**
63     * Construct a Lagrange polynomial with the given abscissas and function
64     * values. The order of interpolating points are not important.
65     * <p>
66     * The constructor makes copy of the input arrays and assigns them.</p>
67     *
68     * @param x interpolating points
69     * @param y function values at interpolating points
70     * @throws IllegalArgumentException if input arrays are not valid
71     */
72    public PolynomialFunctionLagrangeForm(double x[], double y[])
73        throws IllegalArgumentException {
74
75        verifyInterpolationArray(x, y);
76        this.x = new double[x.length];
77        this.y = new double[y.length];
78        System.arraycopy(x, 0, this.x, 0, x.length);
79        System.arraycopy(y, 0, this.y, 0, y.length);
80        coefficientsComputed = false;
81    }
82
83    /** {@inheritDoc} */
84    public double value(double z) throws FunctionEvaluationException {
85        try {
86            return evaluate(x, y, z);
87        } catch (DuplicateSampleAbscissaException e) {
88            throw new FunctionEvaluationException(z, e.getSpecificPattern(), e.getGeneralPattern(), e.getArguments());
89        }
90    }
91
92    /**
93     * Returns the degree of the polynomial.
94     *
95     * @return the degree of the polynomial
96     */
97    public int degree() {
98        return x.length - 1;
99    }
100
101    /**
102     * Returns a copy of the interpolating points array.
103     * <p>
104     * Changes made to the returned copy will not affect the polynomial.</p>
105     *
106     * @return a fresh copy of the interpolating points array
107     */
108    public double[] getInterpolatingPoints() {
109        double[] out = new double[x.length];
110        System.arraycopy(x, 0, out, 0, x.length);
111        return out;
112    }
113
114    /**
115     * Returns a copy of the interpolating values array.
116     * <p>
117     * Changes made to the returned copy will not affect the polynomial.</p>
118     *
119     * @return a fresh copy of the interpolating values array
120     */
121    public double[] getInterpolatingValues() {
122        double[] out = new double[y.length];
123        System.arraycopy(y, 0, out, 0, y.length);
124        return out;
125    }
126
127    /**
128     * Returns a copy of the coefficients array.
129     * <p>
130     * Changes made to the returned copy will not affect the polynomial.</p>
131     * <p>
132     * Note that coefficients computation can be ill-conditioned. Use with caution
133     * and only when it is necessary.</p>
134     *
135     * @return a fresh copy of the coefficients array
136     */
137    public double[] getCoefficients() {
138        if (!coefficientsComputed) {
139            computeCoefficients();
140        }
141        double[] out = new double[coefficients.length];
142        System.arraycopy(coefficients, 0, out, 0, coefficients.length);
143        return out;
144    }
145
146    /**
147     * Evaluate the Lagrange polynomial using
148     * <a href="http://mathworld.wolfram.com/NevillesAlgorithm.html">
149     * Neville's Algorithm</a>. It takes O(N^2) time.
150     * <p>
151     * This function is made public static so that users can call it directly
152     * without instantiating PolynomialFunctionLagrangeForm object.</p>
153     *
154     * @param x the interpolating points array
155     * @param y the interpolating values array
156     * @param z the point at which the function value is to be computed
157     * @return the function value
158     * @throws DuplicateSampleAbscissaException if the sample has duplicate abscissas
159     * @throws IllegalArgumentException if inputs are not valid
160     */
161    public static double evaluate(double x[], double y[], double z) throws
162        DuplicateSampleAbscissaException, IllegalArgumentException {
163
164        verifyInterpolationArray(x, y);
165
166        int nearest = 0;
167        final int n = x.length;
168        final double[] c = new double[n];
169        final double[] d = new double[n];
170        double min_dist = Double.POSITIVE_INFINITY;
171        for (int i = 0; i < n; i++) {
172            // initialize the difference arrays
173            c[i] = y[i];
174            d[i] = y[i];
175            // find out the abscissa closest to z
176            final double dist = FastMath.abs(z - x[i]);
177            if (dist < min_dist) {
178                nearest = i;
179                min_dist = dist;
180            }
181        }
182
183        // initial approximation to the function value at z
184        double value = y[nearest];
185
186        for (int i = 1; i < n; i++) {
187            for (int j = 0; j < n-i; j++) {
188                final double tc = x[j] - z;
189                final double td = x[i+j] - z;
190                final double divider = x[j] - x[i+j];
191                if (divider == 0.0) {
192                    // This happens only when two abscissas are identical.
193                    throw new DuplicateSampleAbscissaException(x[i], i, i+j);
194                }
195                // update the difference arrays
196                final double w = (c[j+1] - d[j]) / divider;
197                c[j] = tc * w;
198                d[j] = td * w;
199            }
200            // sum up the difference terms to get the final value
201            if (nearest < 0.5*(n-i+1)) {
202                value += c[nearest];    // fork down
203            } else {
204                nearest--;
205                value += d[nearest];    // fork up
206            }
207        }
208
209        return value;
210    }
211
212    /**
213     * Calculate the coefficients of Lagrange polynomial from the
214     * interpolation data. It takes O(N^2) time.
215     * <p>
216     * Note this computation can be ill-conditioned. Use with caution
217     * and only when it is necessary.</p>
218     *
219     * @throws ArithmeticException if any abscissas coincide
220     */
221    protected void computeCoefficients() throws ArithmeticException {
222
223        final int n = degree() + 1;
224        coefficients = new double[n];
225        for (int i = 0; i < n; i++) {
226            coefficients[i] = 0.0;
227        }
228
229        // c[] are the coefficients of P(x) = (x-x[0])(x-x[1])...(x-x[n-1])
230        final double[] c = new double[n+1];
231        c[0] = 1.0;
232        for (int i = 0; i < n; i++) {
233            for (int j = i; j > 0; j--) {
234                c[j] = c[j-1] - c[j] * x[i];
235            }
236            c[0] *= -x[i];
237            c[i+1] = 1;
238        }
239
240        final double[] tc = new double[n];
241        for (int i = 0; i < n; i++) {
242            // d = (x[i]-x[0])...(x[i]-x[i-1])(x[i]-x[i+1])...(x[i]-x[n-1])
243            double d = 1;
244            for (int j = 0; j < n; j++) {
245                if (i != j) {
246                    d *= x[i] - x[j];
247                }
248            }
249            if (d == 0.0) {
250                // This happens only when two abscissas are identical.
251                for (int k = 0; k < n; ++k) {
252                    if ((i != k) && (x[i] == x[k])) {
253                        throw MathRuntimeException.createArithmeticException(
254                              LocalizedFormats.IDENTICAL_ABSCISSAS_DIVISION_BY_ZERO,
255                              i, k, x[i]);
256                    }
257                }
258            }
259            final double t = y[i] / d;
260            // Lagrange polynomial is the sum of n terms, each of which is a
261            // polynomial of degree n-1. tc[] are the coefficients of the i-th
262            // numerator Pi(x) = (x-x[0])...(x-x[i-1])(x-x[i+1])...(x-x[n-1]).
263            tc[n-1] = c[n];     // actually c[n] = 1
264            coefficients[n-1] += t * tc[n-1];
265            for (int j = n-2; j >= 0; j--) {
266                tc[j] = c[j+1] + tc[j+1] * x[i];
267                coefficients[j] += t * tc[j];
268            }
269        }
270
271        coefficientsComputed = true;
272    }
273
274    /**
275     * Verifies that the interpolation arrays are valid.
276     * <p>
277     * The arrays features checked by this method are that both arrays have the
278     * same length and this length is at least 2.
279     * </p>
280     * <p>
281     * The interpolating points must be distinct. However it is not
282     * verified here, it is checked in evaluate() and computeCoefficients().
283     * </p>
284     *
285     * @param x the interpolating points array
286     * @param y the interpolating values array
287     * @throws IllegalArgumentException if not valid
288     * @see #evaluate(double[], double[], double)
289     * @see #computeCoefficients()
290     */
291    public static void verifyInterpolationArray(double x[], double y[])
292        throws IllegalArgumentException {
293
294        if (x.length != y.length) {
295            throw MathRuntimeException.createIllegalArgumentException(
296                  LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, x.length, y.length);
297        }
298
299        if (x.length < 2) {
300            throw MathRuntimeException.createIllegalArgumentException(
301                  LocalizedFormats.WRONG_NUMBER_OF_POINTS, 2, x.length);
302        }
303
304    }
305}
306