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.util.MathArrays;
022 import org.apache.commons.math.analysis.DifferentiableUnivariateRealFunction;
023 import org.apache.commons.math.analysis.UnivariateRealFunction;
024 import org.apache.commons.math.exception.OutOfRangeException;
025 import org.apache.commons.math.exception.NumberIsTooSmallException;
026 import org.apache.commons.math.exception.DimensionMismatchException;
027 import org.apache.commons.math.exception.NullArgumentException;
028 import org.apache.commons.math.exception.util.LocalizedFormats;
029
030 /**
031 * Represents a polynomial spline function.
032 * <p>
033 * A <strong>polynomial spline function</strong> consists of a set of
034 * <i>interpolating polynomials</i> and an ascending array of domain
035 * <i>knot points</i>, determining the intervals over which the spline function
036 * is defined by the constituent polynomials. The polynomials are assumed to
037 * have been computed to match the values of another function at the knot
038 * points. The value consistency constraints are not currently enforced by
039 * <code>PolynomialSplineFunction</code> itself, but are assumed to hold among
040 * the polynomials and knot points passed to the constructor.</p>
041 * <p>
042 * N.B.: The polynomials in the <code>polynomials</code> property must be
043 * centered on the knot points to compute the spline function values.
044 * See below.</p>
045 * <p>
046 * The domain of the polynomial spline function is
047 * <code>[smallest knot, largest knot]</code>. Attempts to evaluate the
048 * function at values outside of this range generate IllegalArgumentExceptions.
049 * </p>
050 * <p>
051 * The value of the polynomial spline function for an argument <code>x</code>
052 * is computed as follows:
053 * <ol>
054 * <li>The knot array is searched to find the segment to which <code>x</code>
055 * belongs. If <code>x</code> is less than the smallest knot point or greater
056 * than the largest one, an <code>IllegalArgumentException</code>
057 * is thrown.</li>
058 * <li> Let <code>j</code> be the index of the largest knot point that is less
059 * than or equal to <code>x</code>. The value returned is <br>
060 * <code>polynomials[j](x - knot[j])</code></li></ol></p>
061 *
062 * @version $Id: PolynomialSplineFunction.java 1182134 2011-10-11 22:55:08Z erans $
063 */
064 public class PolynomialSplineFunction implements DifferentiableUnivariateRealFunction {
065 /**
066 * Spline segment interval delimiters (knots).
067 * Size is n + 1 for n segments.
068 */
069 private final double knots[];
070 /**
071 * The polynomial functions that make up the spline. The first element
072 * determines the value of the spline over the first subinterval, the
073 * second over the second, etc. Spline function values are determined by
074 * evaluating these functions at {@code (x - knot[i])} where i is the
075 * knot segment to which x belongs.
076 */
077 private final PolynomialFunction polynomials[];
078 /**
079 * Number of spline segments. It is equal to the number of polynomials and
080 * to the number of partition points - 1.
081 */
082 private final int n;
083
084
085 /**
086 * Construct a polynomial spline function with the given segment delimiters
087 * and interpolating polynomials.
088 * The constructor copies both arrays and assigns the copies to the knots
089 * and polynomials properties, respectively.
090 *
091 * @param knots Spline segment interval delimiters.
092 * @param polynomials Polynomial functions that make up the spline.
093 * @throws NullArgumentException if either of the input arrays is {@code null}.
094 * @throws NumberIsTooSmallException if knots has length less than 2.
095 * @throws DimensionMismatchException if {@code polynomials.length != knots.length - 1}.
096 * @throws org.apache.commons.math.exception.NonMonotonicSequenceException if
097 * the {@code knots} array is not strictly increasing.
098 *
099 */
100 public PolynomialSplineFunction(double knots[], PolynomialFunction polynomials[]) {
101 if (knots == null ||
102 polynomials == null) {
103 throw new NullArgumentException();
104 }
105 if (knots.length < 2) {
106 throw new NumberIsTooSmallException(LocalizedFormats.NOT_ENOUGH_POINTS_IN_SPLINE_PARTITION,
107 2, knots.length, false);
108 }
109 if (knots.length - 1 != polynomials.length) {
110 throw new DimensionMismatchException(polynomials.length, knots.length);
111 }
112 MathArrays.checkOrder(knots);
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.
125 *
126 * @param v Point for which the function value should be computed.
127 * @return the value.
128 * @throws OutOfRangeException if {@code v} is outside of the domain of the
129 * spline function (smaller than the smallest knot point or larger than the
130 * largest knot point).
131 */
132 public double value(double v) {
133 if (v < knots[0] || v > knots[n]) {
134 throw new OutOfRangeException(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 * Get the derivative of the polynomial spline function.
151 *
152 * @return the derivative function.
153 */
154 public UnivariateRealFunction derivative() {
155 return polynomialSplineDerivative();
156 }
157
158 /**
159 * Get the derivative of the polynomial spline function.
160 *
161 * @return the derivative function.
162 */
163 public PolynomialSplineFunction polynomialSplineDerivative() {
164 PolynomialFunction derivativePolynomials[] = new PolynomialFunction[n];
165 for (int i = 0; i < n; i++) {
166 derivativePolynomials[i] = polynomials[i].polynomialDerivative();
167 }
168 return new PolynomialSplineFunction(knots, derivativePolynomials);
169 }
170
171 /**
172 * Get the number of spline segments.
173 * It is also the number of polynomials and the number of knot points - 1.
174 *
175 * @return the number of spline segments.
176 */
177 public int getN() {
178 return n;
179 }
180
181 /**
182 * Get a copy of the interpolating polynomials array.
183 * It returns a fresh copy of the array. Changes made to the copy will
184 * not affect the polynomials property.
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 * Get an array copy of the knot points.
196 * It 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 }