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