    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.math.analysis.polynomials;
    018
    019    import java.util.Arrays;
    020
    021    import org.apache.commons.math.ArgumentOutsideDomainException;
    022    import org.apache.commons.math.MathRuntimeException;
    023    import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
    024    import org.apache.commons.math.analysis.UnivariateRealFunction;
    025    import org.apache.commons.math.exception.util.LocalizedFormats;
    026
    027    /**
    028     * Represents a polynomial spline function.
    029     * <p>
    030     * A <strong>polynomial spline function</strong> consists of a set of
    031     * <i>interpolating polynomials</i> and an ascending array of domain
    032     * <i>knot points</i>, determining the intervals over which the spline function
    033     * is defined by the constituent polynomials.  The polynomials are assumed to
    034     * have been computed to match the values of another function at the knot
    035     * points.  The value consistency constraints are not currently enforced by
    036     * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
    037     * the polynomials and knot points passed to the constructor.</p>
    038     * <p>
    039     * N.B.:  The polynomials in the <code>polynomials</code> property must be
    040     * centered on the knot points to compute the spline function values.
    041     * See below.</p>
    042     * <p>
    043     * The domain of the polynomial spline function is
    044     * <code>[smallest knot, largest knot]</code>.  Attempts to evaluate the
    045     * function at values outside of this range generate IllegalArgumentExceptions.
    046     * </p>
    047     * <p>
    048     * The value of the polynomial spline function for an argument <code>x</code>
    049     * is computed as follows:
    050     * <ol>
    051     * <li>The knot array is searched to find the segment to which <code>x</code>
    052     * belongs.  If <code>x</code> is less than the smallest knot point or greater
    053     * than the largest one, an <code>IllegalArgumentException</code>
    054     * is thrown.</li>
    055     * <li> Let <code>j</code> be the index of the largest knot point that is less
    056     * than or equal to <code>x</code>.  The value returned is <br>
    057     * <code>polynomials[j](x - knot[j])</code></li></ol></p>
    058     *
    059     * @version $Revision: 1037327 $ $Date: 2010-11-20 21:57:37 +0100 (sam. 20 nov. 2010) $
    060     */
    061    public class PolynomialSplineFunction
    062        implements DifferentiableUnivariateRealFunction {
    063
    064        /** Spline segment interval delimiters (knots).   Size is n+1 for n segments. */
    065        private final double knots[];
    066
    067        /**
    068         * The polynomial functions that make up the spline.  The first element
    069         * determines the value of the spline over the first subinterval, the
    070         * second over the second, etc.   Spline function values are determined by
    071         * evaluating these functions at <code>(x - knot[i])</code> where i is the
    072         * knot segment to which x belongs.
    073         */
    074        private final PolynomialFunction polynomials[];
    075
    076        /**
    077         * Number of spline segments = number of polynomials
    078         *  = number of partition points - 1
    079         */
    080        private final int n;
    081
    082
    083        /**
    084         * Construct a polynomial spline function with the given segment delimiters
    085         * and interpolating polynomials.
    086         * <p>
    087         * The constructor copies both arrays and assigns the copies to the knots
    088         * and polynomials properties, respectively.</p>
    089         *
    090         * @param knots spline segment interval delimiters
    091         * @param polynomials polynomial functions that make up the spline
    092         * @throws NullPointerException if either of the input arrays is null
    093         * @throws IllegalArgumentException if knots has length less than 2,
    094         * <code>polynomials.length != knots.length - 1 </code>, or the knots array
    095         * is not strictly increasing.
    096         *
    097         */
    098        public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
    099            if (knots.length < 2) {
    100                throw MathRuntimeException.createIllegalArgumentException(
    101                      LocalizedFormats.NOT_ENOUGH_POINTS_IN_SPLINE_PARTITION,
    102                      2, knots.length);
    103            }
    104            if (knots.length - 1 != polynomials.length) {
    105                throw MathRuntimeException.createIllegalArgumentException(
    106                      LocalizedFormats.POLYNOMIAL_INTERPOLANTS_MISMATCH_SEGMENTS,
    107                      polynomials.length, knots.length);
    108            }
    109            if (!isStrictlyIncreasing(knots)) {
    110                throw MathRuntimeException.createIllegalArgumentException(
    111                      LocalizedFormats.NOT_STRICTLY_INCREASING_KNOT_VALUES);
    112            }
    113
    114            this.n = knots.length -1;
    115            this.knots = new double[n + 1];
    116            System.arraycopy(knots, 0, this.knots, 0, n + 1);
    117            this.polynomials = new PolynomialFunction[n];
    118            System.arraycopy(polynomials, 0, this.polynomials, 0, n);
    119        }
    120
    121        /**
    122         * Compute the value for the function.
    123         * See {@link PolynomialSplineFunction} for details on the algorithm for
    124         * computing the value of the function.</p>
    125         *
    126         * @param v the point for which the function value should be computed
    127         * @return the value
    128         * @throws ArgumentOutsideDomainException if v is outside of the domain of
    129         * of the spline function (less than the smallest knot point or greater
    130         * than the largest knot point)
    131         */
    132        public double value(double v) throws ArgumentOutsideDomainException {
    133            if (v < knots[0] || v > knots[n]) {
    134                throw new ArgumentOutsideDomainException(v, knots[0], knots[n]);
    135            }
    136            int i = Arrays.binarySearch(knots, v);
    137            if (i < 0) {
    138                i = -i - 2;
    139            }
    140            //This will handle the case where v is the last knot value
    141            //There are only n-1 polynomials, so if v is the last knot
    142            //then we will use the last polynomial to calculate the value.
    143            if ( i >= polynomials.length ) {
    144                i--;
    145            }
    146            return polynomials[i].value(v - knots[i]);
    147        }
    148
    149        /**
    150         * Returns the derivative of the polynomial spline function as a UnivariateRealFunction
    151         * @return  the derivative function
    152         */
    153        public UnivariateRealFunction derivative() {
    154            return polynomialSplineDerivative();
    155        }
    156
    157        /**
    158         * Returns the derivative of the polynomial spline function as a PolynomialSplineFunction
    159         *
    160         * @return  the derivative function
    161         */
    162        public PolynomialSplineFunction polynomialSplineDerivative() {
    163            PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
    164            for (int i = 0; i < n; i++) {
    165                derivativePolynomials[i] = polynomials[i].polynomialDerivative();
    166            }
    167            return new PolynomialSplineFunction(knots, derivativePolynomials);
    168        }
    169
    170        /**
    171         * Returns the number of spline segments = the number of polynomials
    172         * = the number of knot points - 1.
    173         *
    174         * @return the number of spline segments
    175         */
    176        public int getN() {
    177            return n;
    178        }
    179
    180        /**
    181         * Returns a copy of the interpolating polynomials array.
    182         * <p>
    183         * Returns a fresh copy of the array. Changes made to the copy will
    184         * not affect the polynomials property.</p>
    185         *
    186         * @return the interpolating polynomials
    187         */
    188        public PolynomialFunction[] getPolynomials() {
    189            PolynomialFunction p[] = new PolynomialFunction[n];
    190            System.arraycopy(polynomials, 0, p, 0, n);
    191            return p;
    192        }
    193
    194        /**
    195         * Returns an array copy of the knot points.
    196         * <p>
    197         * Returns a fresh copy of the array. Changes made to the copy
    198         * will not affect the knots property.</p>
    199         *
    200         * @return the knot points
    201         */
    202        public double[] getKnots() {
    203            double out[] = new double[n + 1];
    204            System.arraycopy(knots, 0, out, 0, n + 1);
    205            return out;
    206        }
    207
    208        /**
    209         * Determines if the given array is ordered in a strictly increasing
    210         * fashion.
    211         *
    212         * @param x the array to examine.
    213         * @return <code>true</code> if the elements in <code>x</code> are ordered
    214         * in a stricly increasing manner.  <code>false</code>, otherwise.
    215         */
    216        private static boolean isStrictlyIncreasing(double[] x) {
    217            for (int i = 1; i < x.length; ++i) {
    218                if (x[i - 1] >= x[i]) {
    219                    return false;
    220                }
    221            }
    222            return true;
    223        }
    224    }