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 }