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 }