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