Coverage Report - org.apache.commons.math.analysis.SplineInterpolator

Classes in this File	Line Coverage	Branch Coverage	Complexity
SplineInterpolator	97%	100%	12.0;12

1		/*
2		* Copyright 2003-2004 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		/**
19		* Computes a natural (a.k.a. "free", "unclamped") cubic spline interpolation for the data set.
20		* <p>
21		* The {@link #interpolate(double[], double[])} method returns a {@link PolynomialSplineFunction}
22		* consisting of n cubic polynomials, defined over the subintervals determined by the x values,
23		* x[0] < x[i] ... < x[n]. The x values are referred to as "knot points."
24		* <p>
25		* The value of the PolynomialSplineFunction at a point x that is greater than or equal to the smallest
26		* knot point and strictly less than the largest knot point is computed by finding the subinterval to which
27		* x belongs and computing the value of the corresponding polynomial at <code>x - x[i] </code> where
28		* <code>i</code> is the index of the subinterval. See {@link PolynomialSplineFunction} for more details.
29		* <p>
30		* The interpolating polynomials satisfy: <ol>
31		* <li>The value of the PolynomialSplineFunction at each of the input x values equals the
32		* corresponding y value.</li>
33		* <li>Adjacent polynomials are equal through two derivatives at the knot points (i.e., adjacent polynomials
34		* "match up" at the knot points, as do their first and second derivatives).</li>
35		* </ol>
36		* <p>
37		* The cubic spline interpolation algorithm implemented is as described in R.L. Burden, J.D. Faires,
38		* <u>Numerical Analysis</u>, 4th Ed., 1989, PWS-Kent, ISBN 0-53491-585-X, pp 126-131.
39		*
40		* @version $Revision$ $Date: 2005-02-26 05:11:52 -0800 (Sat, 26 Feb 2005) $
41		*
42		*/
43	10	public class SplineInterpolator implements UnivariateRealInterpolator {
44
45		/**
46		* Computes an interpolating function for the data set.
47		* @param x the arguments for the interpolation points
48		* @param y the values for the interpolation points
49		* @return a function which interpolates the data set
50		*/
51		public UnivariateRealFunction interpolate(double x[], double y[]) {
52	12	if (x.length != y.length) {
53	2	throw new IllegalArgumentException("Dataset arrays must have same length.");
54		}
55
56	10	if (x.length < 3) {
57	0	throw new IllegalArgumentException
58		("At least 3 datapoints are required to compute a spline interpolant");
59		}
60
61		// Number of intervals. The number of data points is n + 1.
62	10	int n = x.length - 1;
63
64	42	for (int i = 0; i < n; i++) {
65	34	if (x[i] >= x[i + 1]) {
66	2	throw new IllegalArgumentException("Dataset x values must be strictly increasing.");
67		}
68		}
69
70		// Differences between knot points
71	8	double h[] = new double[n];
72	38	for (int i = 0; i < n; i++) {
73	30	h[i] = x[i + 1] - x[i];
74		}
75
76	8	double mu[] = new double[n];
77	8	double z[] = new double[n + 1];
78	8	mu[0] = 0d;
79	8	z[0] = 0d;
80	8	double g = 0;
81	30	for (int i = 1; i < n; i++) {
82	22	g = 2d * (x[i+1] - x[i - 1]) - h[i - 1] * mu[i -1];
83	22	mu[i] = h[i] / g;
84	22	z[i] = (3d * (y[i + 1] * h[i - 1] - y[i] * (x[i + 1] - x[i - 1])+ y[i - 1] * h[i]) /
85		(h[i - 1] * h[i]) - h[i - 1] * z[i - 1]) / g;
86		}
87
88		// cubic spline coefficients -- b is linear, c quadratic, d is cubic (original y's are constants)
89	8	double b[] = new double[n];
90	8	double c[] = new double[n + 1];
91	8	double d[] = new double[n];
92
93	8	z[n] = 0d;
94	8	c[n] = 0d;
95
96	38	for (int j = n -1; j >=0; j--) {
97	30	c[j] = z[j] - mu[j] * c[j + 1];
98	30	b[j] = (y[j + 1] - y[j]) / h[j] - h[j] * (c[j + 1] + 2d * c[j]) / 3d;
99	30	d[j] = (c[j + 1] - c[j]) / (3d * h[j]);
100		}
101
102	8	PolynomialFunction polynomials[] = new PolynomialFunction[n];
103	8	double coefficients[] = new double[4];
104	38	for (int i = 0; i < n; i++) {
105	30	coefficients[0] = y[i];
106	30	coefficients[1] = b[i];
107	30	coefficients[2] = c[i];
108	30	coefficients[3] = d[i];
109	30	polynomials[i] = new PolynomialFunction(coefficients);
110		}
111
112	8	return new PolynomialSplineFunction(x, polynomials);
113		}
114
115		}