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 import org.apache.commons.math.exception.util.LocalizedFormats; 26 27 /** 28 * Represents a polynomial spline function. 29 * <p> 30 * A <strong>polynomial spline function</strong> consists of a set of 31 * <i>interpolating polynomials</i> and an ascending array of domain 32 * <i>knot points</i>, determining the intervals over which the spline function 33 * is defined by the constituent polynomials. The polynomials are assumed to 34 * have been computed to match the values of another function at the knot 35 * points. The value consistency constraints are not currently enforced by 36 * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among 37 * the polynomials and knot points passed to the constructor.</p> 38 * <p> 39 * N.B.: The polynomials in the <code>polynomials</code> property must be 40 * centered on the knot points to compute the spline function values. 41 * See below.</p> 42 * <p> 43 * The domain of the polynomial spline function is 44 * <code>[smallest knot, largest knot]</code>. Attempts to evaluate the 45 * function at values outside of this range generate IllegalArgumentExceptions. 46 * </p> 47 * <p> 48 * The value of the polynomial spline function for an argument <code>x</code> 49 * is computed as follows: 50 * <ol> 51 * <li>The knot array is searched to find the segment to which <code>x</code> 52 * belongs. If <code>x</code> is less than the smallest knot point or greater 53 * than the largest one, an <code>IllegalArgumentException</code> 54 * is thrown.</li> 55 * <li> Let <code>j</code> be the index of the largest knot point that is less 56 * than or equal to <code>x</code>. The value returned is <br> 57 * <code>polynomials[j](x - knot[j])</code></li></ol></p> 58 * 59 * @version $Revision: 1037327 $ $Date: 2010-11-20 21:57:37 +0100 (sam. 20 nov. 2010) $ 60 */ 61 public class PolynomialSplineFunction 62 implements DifferentiableUnivariateRealFunction { 63 64 /** Spline segment interval delimiters (knots). Size is n+1 for n segments. */ 65 private final 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 final PolynomialFunction polynomials[]; 75 76 /** 77 * Number of spline segments = number of polynomials 78 * = number of partition points - 1 79 */ 80 private final int n; 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 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 } 225
