View Javadoc

1   /*
2    * Copyright 2003-2005 The Apache Software Foundation.
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *      http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
15   */
16  package org.apache.commons.math.analysis;
17  
18  import java.io.Serializable;
19  import java.util.Arrays;
20  
21  import org.apache.commons.math.FunctionEvaluationException;
22  
23  /**
24   * Represents a polynomial spline function.
25   * <p>
26   * A <strong>polynomial spline function</strong> consists of a set of 
27   * <i>interpolating polynomials</i> and an ascending array of domain 
28   * <i>knot points</i>, determining the intervals over which the spline function
29   * is defined by the constituent polynomials.  The polynomials are assumed to
30   * have been computed to match the values of another function at the knot
31   * points.  The value consistency constraints are not currently enforced by 
32   * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
33   * the polynomials and knot points passed to the constructor.
34   * <p>
35   * N.B.:  The polynomials in the <code>polynomials</code> property must be
36   * centered on the knot points to compute the spline function values.  See below.
37   * <p>
38   * The domain of the polynomial spline function is 
39   * <code>[smallest knot, largest knot]</code>.  Attempts to evaluate the
40   * function at values outside of this range generate IllegalArgumentExceptions.
41   * <p>
42   * The value of the polynomial spline function for an argument <code>x</code>
43   * is computed as follows:
44   * <ol>
45   * <li>The knot array is searched to find the segment to which <code>x</code>
46   * belongs.  If <code>x</code> is less than the smallest knot point or greater
47   * than the largest one, an <code>IllegalArgumentException</code>
48   * is thrown.</li>
49   * <li> Let <code>j</code> be the index of the largest knot point that is less
50   * than or equal to <code>x</code>.  The value returned is <br>
51   * <code>polynomials[j](x - knot[j])</code></li></ol>
52   *
53   * @version $Revision: 348761 $ $Date: 2005-11-24 09:04:20 -0700 (Thu, 24 Nov 2005) $
54   */
55  public class PolynomialSplineFunction 
56      implements DifferentiableUnivariateRealFunction, Serializable {
57     
58      /** Serializable version identifier */
59      private static final long serialVersionUID = 7011031166416885789L;
60      
61      /** Spline segment interval delimiters (knots).   Size is n+1 for n segments. */
62      private double knots[];
63  
64      /**
65       * The polynomial functions that make up the spline.  The first element
66       * determines the value of the spline over the first subinterval, the
67       * second over the second, etc.   Spline function values are determined by
68       * evaluating these functions at <code>(x - knot[i])</code> where i is the
69       * knot segment to which x belongs.
70       */
71      private PolynomialFunction polynomials[] = null;
72      
73      /** 
74       * Number of spline segments = number of polynomials
75       *  = number of partition points - 1 
76       */
77      private int n = 0;
78      
79  
80      /**
81       * Construct a polynomial spline function with the given segment delimiters
82       * and interpolating polynomials.
83       * <p>
84       * The constructor copies both arrays and assigns the copies to the knots
85       * and polynomials properties, respectively.
86       * 
87       * @param knots spline segment interval delimiters
88       * @param polynomials polynomial functions that make up the spline
89       * @throws NullPointerException if either of the input arrays is null
90       * @throws IllegalArgumentException if knots has length less than 2,  
91       * <code>polynomials.length != knots.length - 1 </code>, or the knots array
92       * is not strictly increasing.
93       * 
94       */
95      public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
96          if (knots.length < 2) {
97              throw new IllegalArgumentException
98                  ("Not enough knot values -- spline partition must have at least 2 points.");
99          }
100         if (knots.length - 1 != polynomials.length) {
101             throw new IllegalArgumentException 
102             ("Number of polynomial interpolants must match the number of segments.");
103         }
104         if (!isStrictlyIncreasing(knots)) {
105             throw new IllegalArgumentException 
106                 ("Knot values must be strictly increasing.");
107         }
108         
109         this.n = knots.length -1;
110         this.knots = new double[n + 1];
111         System.arraycopy(knots, 0, this.knots, 0, n + 1);
112         this.polynomials = new PolynomialFunction[n];
113         System.arraycopy(polynomials, 0, this.polynomials, 0, n);
114     }
115 
116     /**
117      * Compute the value for the function.
118      * <p>
119      * Throws FunctionEvaluationException if v is outside of the domain of the
120      * function.  The domain is [smallest knot, largest knot].
121      * <p>
122      * See {@link PolynomialSplineFunction} for details on the algorithm for
123      * computing the value of the function.
124      * 
125      * @param v the point for which the function value should be computed
126      * @return the value
127      * @throws FunctionEvaluationException if v is outside of the domain of
128      * of the spline function (less than the smallest knot point or greater
129      * than the largest knot point)
130      */
131     public double value(double v) throws FunctionEvaluationException {
132         if (v < knots[0] || v > knots[n]) {
133             throw new FunctionEvaluationException(v,"Argument outside domain");
134         }
135         int i = Arrays.binarySearch(knots, v);
136         if (i < 0) {
137             i = -i - 2;
138         }
139         //This will handle the case where v is the last knot value
140         //There are only n-1 polynomials, so if v is the last knot
141         //then we will use the last polynomial to calculate the value.
142         if ( i >= polynomials.length ) {
143             i--;
144         }
145         return polynomials[i].value(v - knots[i]);
146     }
147     
148     /**
149      * Returns the derivative of the polynomial spline function as a UnivariateRealFunction
150      * @return  the derivative function
151      */
152     public UnivariateRealFunction derivative() {
153         return polynomialSplineDerivative();
154     }
155     
156     /**
157      * Returns the derivative of the polynomial spline function as a PolynomialSplineFunction
158      * 
159      * @return  the derivative function
160      */
161     public PolynomialSplineFunction polynomialSplineDerivative() {
162         PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
163         for (int i = 0; i < n; i++) {
164             derivativePolynomials[i] = polynomials[i].polynomialDerivative();
165         }
166         return new PolynomialSplineFunction(knots, derivativePolynomials);
167     }
168 
169     /**
170      * Returns the number of spline segments = the number of polynomials 
171      * = the number of knot points - 1.
172      * 
173      * @return the number of spline segments
174      */
175     public int getN() {
176         return n;
177     }
178 
179     /**
180      * Returns a copy of the interpolating polynomials array.
181      * <p>
182      * Returns a fresh copy of the array. Changes made to the copy will
183      * not affect the polynomials property.
184      * 
185      * @return the interpolating polynomials
186      */
187     public PolynomialFunction[] getPolynomials() {
188         PolynomialFunction p[] = new PolynomialFunction[n];
189         System.arraycopy(polynomials, 0, p, 0, n);
190         return p;
191     }
192 
193     /**
194      * Returns an array copy of the knot points.
195      * <p>
196      * Returns a fresh copy of the array. Changes made to the copy
197      * will not affect the knots property.
198      * 
199      * @return the knot points
200      */
201     public double[] getKnots() {
202         double out[] = new double[n + 1];
203         System.arraycopy(knots, 0, out, 0, n + 1);
204         return out;  
205     }
206 
207     /**
208      * Determines if the given array is ordered in a strictly increasing
209      * fashion.
210      * 
211      * @param x the array to examine.
212      * @return <code>true</code> if the elements in <code>x</code> are ordered
213      * in a stricly increasing manner.  <code>false</code>, otherwise.
214      */
215     private static boolean isStrictlyIncreasing(double[] x) {
216         for (int i = 1; i < x.length; ++i) {
217             if (x[i - 1] >= x[i]) {
218                 return false;
219             }
220         }
221         return true;
222     }
223 }