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 */
066public 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 NonMonotonicSequenceException if the {@code knots} array is not strictly increasing.
099     *
100     */
101    public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[])
102        throws NullArgumentException, NumberIsTooSmallException,
103               DimensionMismatchException, NonMonotonicSequenceException{
104        if (knots == null ||
105            polynomials == null) {
106            throw new NullArgumentException();
107        }
108        if (knots.length < 2) {
109            throw new NumberIsTooSmallException(LocalizedFormats.NOT_ENOUGH_POINTS_IN_SPLINE_PARTITION,
110                                                2, knots.length, false);
111        }
112        if (knots.length - 1 != polynomials.length) {
113            throw new DimensionMismatchException(polynomials.length, knots.length);
114        }
115        MathArrays.checkOrder(knots);
116
117        this.n = knots.length -1;
118        this.knots = new double[n + 1];
119        System.arraycopy(knots, 0, this.knots, 0, n + 1);
120        this.polynomials = new PolynomialFunction[n];
121        System.arraycopy(polynomials, 0, this.polynomials, 0, n);
122    }
123
124    /**
125     * Compute the value for the function.
126     * See {@link PolynomialSplineFunction} for details on the algorithm for
127     * computing the value of the function.
128     *
129     * @param v Point for which the function value should be computed.
130     * @return the value.
131     * @throws OutOfRangeException if {@code v} is outside of the domain of the
132     * spline function (smaller than the smallest knot point or larger than the
133     * largest knot point).
134     */
135    public double value(double v) {
136        if (v < knots[0] || v > knots[n]) {
137            throw new OutOfRangeException(v, knots[0], knots[n]);
138        }
139        int i = Arrays.binarySearch(knots, v);
140        if (i < 0) {
141            i = -i - 2;
142        }
143        // This will handle the case where v is the last knot value
144        // There are only n-1 polynomials, so if v is the last knot
145        // then we will use the last polynomial to calculate the value.
146        if ( i >= polynomials.length ) {
147            i--;
148        }
149        return polynomials[i].value(v - knots[i]);
150    }
151
152    /**
153     * Get the derivative of the polynomial spline function.
154     *
155     * @return the derivative function.
156     */
157    public UnivariateFunction derivative() {
158        return polynomialSplineDerivative();
159    }
160
161    /**
162     * Get the derivative of the polynomial spline function.
163     *
164     * @return the derivative function.
165     */
166    public PolynomialSplineFunction polynomialSplineDerivative() {
167        PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
168        for (int i = 0; i < n; i++) {
169            derivativePolynomials[i] = polynomials[i].polynomialDerivative();
170        }
171        return new PolynomialSplineFunction(knots, derivativePolynomials);
172    }
173
174
175    /** {@inheritDoc}
176     * @since 3.1
177     */
178    public DerivativeStructure value(final DerivativeStructure t) {
179        final double t0 = t.getValue();
180        if (t0 < knots[0] || t0 > knots[n]) {
181            throw new OutOfRangeException(t0, knots[0], knots[n]);
182        }
183        int i = Arrays.binarySearch(knots, t0);
184        if (i < 0) {
185            i = -i - 2;
186        }
187        // This will handle the case where t is the last knot value
188        // There are only n-1 polynomials, so if t is the last knot
189        // then we will use the last polynomial to calculate the value.
190        if ( i >= polynomials.length ) {
191            i--;
192        }
193        return polynomials[i].value(t.subtract(knots[i]));
194    }
195
196    /**
197     * Get the number of spline segments.
198     * It is also the number of polynomials and the number of knot points - 1.
199     *
200     * @return the number of spline segments.
201     */
202    public int getN() {
203        return n;
204    }
205
206    /**
207     * Get a copy of the interpolating polynomials array.
208     * It returns a fresh copy of the array. Changes made to the copy will
209     * not affect the polynomials property.
210     *
211     * @return the interpolating polynomials.
212     */
213    public PolynomialFunction[] getPolynomials() {
214        PolynomialFunction p[] = new PolynomialFunction[n];
215        System.arraycopy(polynomials, 0, p, 0, n);
216        return p;
217    }
218
219    /**
220     * Get an array copy of the knot points.
221     * It returns a fresh copy of the array. Changes made to the copy
222     * will not affect the knots property.
223     *
224     * @return the knot points.
225     */
226    public double[] getKnots() {
227        double out[] = new double[n + 1];
228        System.arraycopy(knots, 0, out, 0, n + 1);
229        return out;
230    }
231
232    /**
233     * Indicates whether a point is within the interpolation range.
234     *
235     * @param x Point.
236     * @return {@code true} if {@code x} is a valid point.
237     */
238    public boolean isValidPoint(double x) {
239        if (x < knots[0] ||
240            x > knots[n]) {
241            return false;
242        } else {
243            return true;
244        }
245    }
246}