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