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