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 }