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    }