View Javadoc

1   /*
2    * Licensed to the Apache Software Foundation (ASF) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * The ASF licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *      http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.apache.commons.math.analysis.polynomials;
18  
19  import java.util.Arrays;
20  
21  import org.apache.commons.math.ArgumentOutsideDomainException;
22  import org.apache.commons.math.MathRuntimeException;
23  import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
24  import org.apache.commons.math.analysis.UnivariateRealFunction;
25  
26  /**
27   * Represents a polynomial spline function.
28   * <p>
29   * A <strong>polynomial spline function</strong> consists of a set of 
30   * <i>interpolating polynomials</i> and an ascending array of domain 
31   * <i>knot points</i>, determining the intervals over which the spline function
32   * is defined by the constituent polynomials.  The polynomials are assumed to
33   * have been computed to match the values of another function at the knot
34   * points.  The value consistency constraints are not currently enforced by 
35   * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
36   * the polynomials and knot points passed to the constructor.</p>
37   * <p>
38   * N.B.:  The polynomials in the <code>polynomials</code> property must be
39   * centered on the knot points to compute the spline function values.  
40   * See below.</p>
41   * <p>
42   * The domain of the polynomial spline function is 
43   * <code>[smallest knot, largest knot]</code>.  Attempts to evaluate the
44   * function at values outside of this range generate IllegalArgumentExceptions.
45   * </p>
46   * <p>
47   * The value of the polynomial spline function for an argument <code>x</code>
48   * is computed as follows:
49   * <ol>
50   * <li>The knot array is searched to find the segment to which <code>x</code>
51   * belongs.  If <code>x</code> is less than the smallest knot point or greater
52   * than the largest one, an <code>IllegalArgumentException</code>
53   * is thrown.</li>
54   * <li> Let <code>j</code> be the index of the largest knot point that is less
55   * than or equal to <code>x</code>.  The value returned is <br>
56   * <code>polynomials[j](x - knot[j])</code></li></ol></p>
57   *
58   * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $
59   */
60  public class PolynomialSplineFunction 
61      implements DifferentiableUnivariateRealFunction {
62  
63      /** Spline segment interval delimiters (knots).   Size is n+1 for n segments. */
64      private double knots[];
65  
66      /**
67       * The polynomial functions that make up the spline.  The first element
68       * determines the value of the spline over the first subinterval, the
69       * second over the second, etc.   Spline function values are determined by
70       * evaluating these functions at <code>(x - knot[i])</code> where i is the
71       * knot segment to which x belongs.
72       */
73      private PolynomialFunction polynomials[] = null;
74      
75      /** 
76       * Number of spline segments = number of polynomials
77       *  = number of partition points - 1 
78       */
79      private int n = 0;
80      
81  
82      /**
83       * Construct a polynomial spline function with the given segment delimiters
84       * and interpolating polynomials.
85       * <p>
86       * The constructor copies both arrays and assigns the copies to the knots
87       * and polynomials properties, respectively.</p>
88       * 
89       * @param knots spline segment interval delimiters
90       * @param polynomials polynomial functions that make up the spline
91       * @throws NullPointerException if either of the input arrays is null
92       * @throws IllegalArgumentException if knots has length less than 2,  
93       * <code>polynomials.length != knots.length - 1 </code>, or the knots array
94       * is not strictly increasing.
95       * 
96       */
97      public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
98          if (knots.length < 2) {
99              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 }