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.math.analysis.polynomials;
018    
019    import java.util.Arrays;
020    
021    import org.apache.commons.math.ArgumentOutsideDomainException;
022    import org.apache.commons.math.MathRuntimeException;
023    import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
024    import org.apache.commons.math.analysis.UnivariateRealFunction;
025    
026    /**
027     * Represents a polynomial spline function.
028     * <p>
029     * A <strong>polynomial spline function</strong> consists of a set of 
030     * <i>interpolating polynomials</i> and an ascending array of domain 
031     * <i>knot points</i>, determining the intervals over which the spline function
032     * is defined by the constituent polynomials.  The polynomials are assumed to
033     * have been computed to match the values of another function at the knot
034     * points.  The value consistency constraints are not currently enforced by 
035     * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
036     * the polynomials and knot points passed to the constructor.</p>
037     * <p>
038     * N.B.:  The polynomials in the <code>polynomials</code> property must be
039     * centered on the knot points to compute the spline function values.  
040     * See below.</p>
041     * <p>
042     * The domain of the polynomial spline function is 
043     * <code>[smallest knot, largest knot]</code>.  Attempts to evaluate the
044     * function at values outside of this range generate IllegalArgumentExceptions.
045     * </p>
046     * <p>
047     * The value of the polynomial spline function for an argument <code>x</code>
048     * is computed as follows:
049     * <ol>
050     * <li>The knot array is searched to find the segment to which <code>x</code>
051     * belongs.  If <code>x</code> is less than the smallest knot point or greater
052     * than the largest one, an <code>IllegalArgumentException</code>
053     * is thrown.</li>
054     * <li> Let <code>j</code> be the index of the largest knot point that is less
055     * than or equal to <code>x</code>.  The value returned is <br>
056     * <code>polynomials[j](x - knot[j])</code></li></ol></p>
057     *
058     * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
059     */
060    public class PolynomialSplineFunction 
061        implements DifferentiableUnivariateRealFunction {
062    
063        /** Spline segment interval delimiters (knots).   Size is n+1 for n segments. */
064        private double knots[];
065    
066        /**
067         * The polynomial functions that make up the spline.  The first element
068         * determines the value of the spline over the first subinterval, the
069         * second over the second, etc.   Spline function values are determined by
070         * evaluating these functions at <code>(x - knot[i])</code> where i is the
071         * knot segment to which x belongs.
072         */
073        private PolynomialFunction polynomials[] = null;
074        
075        /** 
076         * Number of spline segments = number of polynomials
077         *  = number of partition points - 1 
078         */
079        private int n = 0;
080        
081    
082        /**
083         * Construct a polynomial spline function with the given segment delimiters
084         * and interpolating polynomials.
085         * <p>
086         * The constructor copies both arrays and assigns the copies to the knots
087         * and polynomials properties, respectively.</p>
088         * 
089         * @param knots spline segment interval delimiters
090         * @param polynomials polynomial functions that make up the spline
091         * @throws NullPointerException if either of the input arrays is null
092         * @throws IllegalArgumentException if knots has length less than 2,  
093         * <code>polynomials.length != knots.length - 1 </code>, or the knots array
094         * is not strictly increasing.
095         * 
096         */
097        public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
098            if (knots.length < 2) {
099                throw MathRuntimeException.createIllegalArgumentException(
100                      "spline partition must have at least {0} points, got {1}",
101                      2, knots.length);
102            }
103            if (knots.length - 1 != polynomials.length) {
104                throw MathRuntimeException.createIllegalArgumentException(
105                      "number of polynomial interpolants must match the number of segments ({0} != {1} - 1)",
106                      polynomials.length, knots.length);
107            }
108            if (!isStrictlyIncreasing(knots)) {
109                throw MathRuntimeException.createIllegalArgumentException(
110                      "knot values must be strictly increasing");
111            }
112            
113            this.n = knots.length -1;
114            this.knots = new double[n + 1];
115            System.arraycopy(knots, 0, this.knots, 0, n + 1);
116            this.polynomials = new PolynomialFunction[n];
117            System.arraycopy(polynomials, 0, this.polynomials, 0, n);
118        }
119    
120        /**
121         * Compute the value for the function.
122         * <p>
123         * Throws FunctionEvaluationException if v is outside of the domain of the
124         * function.  The domain is [smallest knot, largest knot].</p>
125         * <p>
126         * See {@link PolynomialSplineFunction} for details on the algorithm for
127         * computing the value of the function.</p>
128         * 
129         * @param v the point for which the function value should be computed
130         * @return the value
131         * @throws ArgumentOutsideDomainException if v is outside of the domain of
132         * of the spline function (less than the smallest knot point or greater
133         * than the largest knot point)
134         */
135        public double value(double v) throws ArgumentOutsideDomainException {
136            if (v < knots[0] || v > knots[n]) {
137                throw new ArgumentOutsideDomainException(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         * Returns the derivative of the polynomial spline function as a UnivariateRealFunction
154         * @return  the derivative function
155         */
156        public UnivariateRealFunction derivative() {
157            return polynomialSplineDerivative();
158        }
159        
160        /**
161         * Returns the derivative of the polynomial spline function as a PolynomialSplineFunction
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         * Returns the number of spline segments = the number of polynomials 
175         * = the number of knot points - 1.
176         * 
177         * @return the number of spline segments
178         */
179        public int getN() {
180            return n;
181        }
182    
183        /**
184         * Returns a copy of the interpolating polynomials array.
185         * <p>
186         * Returns a fresh copy of the array. Changes made to the copy will
187         * not affect the polynomials property.</p>
188         * 
189         * @return the interpolating polynomials
190         */
191        public PolynomialFunction[] getPolynomials() {
192            PolynomialFunction p[] = new PolynomialFunction[n];
193            System.arraycopy(polynomials, 0, p, 0, n);
194            return p;
195        }
196    
197        /**
198         * Returns an array copy of the knot points.
199         * <p>
200         * Returns a fresh copy of the array. Changes made to the copy
201         * will not affect the knots property.</p>
202         * 
203         * @return the knot points
204         */
205        public double[] getKnots() {
206            double out[] = new double[n + 1];
207            System.arraycopy(knots, 0, out, 0, n + 1);
208            return out;  
209        }
210    
211        /**
212         * Determines if the given array is ordered in a strictly increasing
213         * fashion.
214         * 
215         * @param x the array to examine.
216         * @return <code>true</code> if the elements in <code>x</code> are ordered
217         * in a stricly increasing manner.  <code>false</code>, otherwise.
218         */
219        private static boolean isStrictlyIncreasing(double[] x) {
220            for (int i = 1; i < x.length; ++i) {
221                if (x[i - 1] >= x[i]) {
222                    return false;
223                }
224            }
225            return true;
226        }
227    }