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