Skip to content

Commit 73ce32d

Browse files
author
Gilles Sadowski
committed
MATH-1650: Clamped spline interpolation (WIP).
1 parent e016130 commit 73ce32d

File tree

4 files changed

+343
-0
lines changed

4 files changed

+343
-0
lines changed

‎commons-math-legacy/src/main/java/org/apache/commons/math4/legacy/analysis/interpolation/ClampedSplineInterpolator.java‎

Lines changed: 161 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,161 @@
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.math4.legacy.analysis.interpolation;
18+
19+
import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
20+
import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialSplineFunction;
21+
import org.apache.commons.math4.legacy.core.MathArrays;
22+
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
23+
import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
24+
import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
25+
import org.apache.commons.math4.legacy.exception.util.LocalizedFormats;
26+
27+
/**
28+
* Clamped cubic spline interpolator.
29+
* The interpolating function consists in cubic polynomial functions defined over the
30+
* subintervals determined by the "knot points".
31+
*
32+
* The interpolating polynomials satisfy:
33+
* <ol>
34+
* <li>The value of the interpolating function at each of the input {@code x} values
35+
* equals the corresponding input {@code y} value.</li>
36+
* <li>Adjacent polynomials are equal through two derivatives at the knot points
37+
* (i.e., adjacent polynomials "match up" at the knot points, as do their first and
38+
* second derivatives).</li>
39+
* <li>The clamped boundary condition, i.e. the interpolating function takes "a
40+
* specific direction" at both its start point and its end point by providing the
41+
* desired first derivative values (slopes) as function parameters to
42+
* {@link #interpolate(double[], double[], double, double)}.</li>
43+
* </ol>
44+
*
45+
* The algorithm is implemented as described in
46+
* <blockquote>
47+
* R.L. Burden, J.D. Faires,
48+
* <em>Numerical Analysis</em>, 9th Ed., 2010, Cengage Learning, ISBN 0-538-73351-9, pp 153-156.
49+
* </blockquote>
50+
*
51+
*/
52+
public class ClampedSplineInterpolator implements UnivariateInterpolator {
53+
/**
54+
*
55+
* The first derivatives evaluated at the first and last knot points are
56+
* approximated from a natural/unclamped spline that passes through the same
57+
* set of points.
58+
*
59+
* @param x Arguments for the interpolation points.
60+
* @param y Values for the interpolation points.
61+
* @return the interpolating function.
62+
* @throws DimensionMismatchException if {@code x} and {@code y} have different sizes.
63+
* @throws NumberIsTooSmallException if the size of {@code x < 3}.
64+
* @throws NonMonotonicSequenceException if {@code x} is not sorted in strict increasing order.
65+
*/
66+
@Override
67+
public PolynomialSplineFunction interpolate(final double[] x,
68+
final double[] y) {
69+
final SplineInterpolator spliner = new SplineInterpolator();
70+
final PolynomialSplineFunction spline = spliner.interpolate(x, y);
71+
final PolynomialSplineFunction derivativeSpline = spline.polynomialSplineDerivative();
72+
final double fpStart = derivativeSpline.value(x[0]);
73+
final double fpEnd = derivativeSpline.value(x[x.length - 1]);
74+
75+
return this.interpolate(x, y, fpStart, fpEnd);
76+
}
77+
78+
/**
79+
* Computes an interpolating function for the data set with defined
80+
* boundary conditions.
81+
*
82+
* @param x Arguments for the interpolation points.
83+
* @param y Values for the interpolation points.
84+
* @param fpStart First derivative at the starting point of the returned
85+
* spline function (starting slope).
86+
* @param fpEnd First derivative at the ending point of the returned
87+
* spline function (ending slope).
88+
* @return the interpolating function.
89+
* @throws DimensionMismatchException if {@code x} and {@code y} have different sizes.
90+
* @throws NumberIsTooSmallException if the size of {@code x < 3}.
91+
* @throws NonMonotonicSequenceException if {@code x} is not sorted in strict increasing order.
92+
*/
93+
public PolynomialSplineFunction interpolate(final double[] x,
94+
final double[] y,
95+
final double fpStart,
96+
final double fpEnd) {
97+
if (x.length != y.length) {
98+
throw new DimensionMismatchException(x.length, y.length);
99+
}
100+
101+
if (x.length < 3) {
102+
throw new NumberIsTooSmallException(LocalizedFormats.NUMBER_OF_POINTS,
103+
x.length, 3, true);
104+
}
105+
106+
// Number of intervals. The number of data points is n + 1.
107+
final int n = x.length - 1;
108+
109+
MathArrays.checkOrder(x);
110+
111+
// Differences between knot points.
112+
final double[] h = new double[n];
113+
for (int i = 0; i < n; i++) {
114+
h[i] = x[i + 1] - x[i];
115+
}
116+
117+
final double[] mu = new double[n];
118+
final double[] z = new double[n + 1];
119+
120+
final double alpha0 = 3d * (y[1] - y[0]) / h[0] - 3d * fpStart;
121+
final double alphaN = 3d * fpEnd - 3d * (y[n] - y[n - 1]) / h[n - 1];
122+
123+
mu[0] = 0.5d;
124+
final double ell0 = 2d * h[0];
125+
z[0] = alpha0 / ell0;
126+
127+
for (int i = 1; i < n; i++) {
128+
final double alpha = (3d / h[i]) * (y[i + 1] - y[i]) - (3d / h[i - 1]) * (y[i] - y[i - 1]);
129+
final double ell = 2d * (x[i + 1] - x[i - 1]) - h[i - 1] * mu[i - 1];
130+
mu[i] = h[i] / ell;
131+
z[i] = (alpha - h[i - 1] * z[i - 1]) / ell;
132+
}
133+
134+
// Cubic spline coefficients.
135+
final double[] b = new double[n]; // Linear.
136+
final double[] c = new double[n + 1]; // Quadratic.
137+
final double[] d = new double[n]; // Cubic.
138+
139+
final double ellN = h[n - 1] * (2d - mu[n - 1]);
140+
z[n] = (alphaN - h[n - 1] * z[n - 1]) / ellN;
141+
c[n] = z[n];
142+
143+
for (int j = n - 1; j >= 0; j--) {
144+
c[j] = z[j] - mu[j] * c[j + 1];
145+
b[j] = ((y[j + 1] - y[j]) / h[j]) - h[j] * (c[j + 1] + 2d * c[j]) / 3d;
146+
d[j] = (c[j + 1] - c[j]) / (3d * h[j]);
147+
}
148+
149+
final PolynomialFunction[] polynomials = new PolynomialFunction[n];
150+
final double[] coefficients = new double[4];
151+
for (int i = 0; i < n; i++) {
152+
coefficients[0] = y[i];
153+
coefficients[1] = b[i];
154+
coefficients[2] = c[i];
155+
coefficients[3] = d[i];
156+
polynomials[i] = new PolynomialFunction(coefficients);
157+
}
158+
159+
return new PolynomialSplineFunction(x, polynomials);
160+
}
161+
}

‎commons-math-legacy/src/test/java/org/apache/commons/math4/legacy/analysis/interpolation/ClampedSplineInterpolatorTest.java‎

Lines changed: 176 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,176 @@
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.math4.legacy.analysis.interpolation;
18+
19+
import org.apache.commons.math4.legacy.TestUtils;
20+
import org.apache.commons.math4.legacy.analysis.UnivariateFunction;
21+
import org.apache.commons.math4.legacy.analysis.integration.SimpsonIntegrator;
22+
import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialFunction;
23+
import org.apache.commons.math4.legacy.analysis.polynomials.PolynomialSplineFunction;
24+
import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
25+
import org.apache.commons.math4.legacy.exception.NonMonotonicSequenceException;
26+
import org.apache.commons.math4.legacy.exception.NumberIsTooSmallException;
27+
import org.junit.Assert;
28+
import org.junit.Test;
29+
30+
/**
31+
* Test the ClampedSplineInterpolator.
32+
*/
33+
public class ClampedSplineInterpolatorTest {
34+
/** Error tolerance for spline interpolator value at knot points. */
35+
private static final double KNOT_TOL = 1e-14;
36+
/** Error tolerance for interpolating polynomial coefficients. */
37+
private static final double COEF_TOL = 1e-14;
38+
39+
@Test
40+
public void testInterpolateLinearDegenerateTwoSegment() {
41+
final double[] x = {0, 0.5, 1};
42+
final double[] y = {1, Math.exp(0.5), Math.exp(1)};
43+
final double fpo = 1;
44+
final double fpn = Math.exp(1);
45+
final ClampedSplineInterpolator i = new ClampedSplineInterpolator();
46+
final PolynomialSplineFunction f = i.interpolate(x, y, fpo, fpn);
47+
verifyInterpolation(f, x, y);
48+
verifyConsistency(f, x);
49+
50+
// Verify coefficients using analytical values
51+
final PolynomialFunction[] polynomials = f.getPolynomials();
52+
53+
final double[] target0 = {1, 1, 0.4889506772539256, 0.21186881109317435};
54+
final double[] target1 = {1.6487212707001282, 1.6478522855738063, 0.8067538938936871, 0.35156753198873575};
55+
TestUtils.assertEquals(polynomials[0].getCoefficients(), target0, COEF_TOL);
56+
TestUtils.assertEquals(polynomials[1].getCoefficients(), target1, COEF_TOL);
57+
}
58+
59+
@Test
60+
public void testInterpolateLinearDegenerateThreeSegment() {
61+
final double[] x = {0, 1, 2, 3};
62+
final double[] y = {1, Math.exp(1), Math.exp(2), Math.exp(3)};
63+
final double fpo = 1;
64+
final double fpn = Math.exp(3);
65+
final ClampedSplineInterpolator i = new ClampedSplineInterpolator();
66+
final PolynomialSplineFunction f = i.interpolate(x, y, fpo, fpn);
67+
verifyInterpolation(f, x, y);
68+
verifyConsistency(f, x);
69+
70+
// Verify coefficients using analytical values
71+
final PolynomialFunction[] polynomials = f.getPolynomials();
72+
73+
final double[] target0 = {1, 0.9999999999999999, 0.4446824969658283, 0.27359933149321697};
74+
final double[] target1 = {2.718281828459045, 2.710162988411307, 1.2654804914454791, 0.6951307906148195};
75+
final double[] target2 = {7.38905609893065, 7.326516343146723, 3.3508728632899376, 2.019091617820356};
76+
TestUtils.assertEquals(polynomials[0].getCoefficients(), target0, COEF_TOL);
77+
TestUtils.assertEquals(polynomials[1].getCoefficients(), target1, COEF_TOL);
78+
TestUtils.assertEquals(polynomials[2].getCoefficients(), target2, COEF_TOL);
79+
}
80+
81+
@Test(expected=DimensionMismatchException.class)
82+
public void testArrayLengthMismatch() {
83+
// Data set arrays of different size.
84+
new ClampedSplineInterpolator().interpolate(new double[] {1, 2, 3, 4},
85+
new double[] {2, 3, 5},
86+
2, 1);
87+
}
88+
@Test(expected=NonMonotonicSequenceException.class)
89+
public void testUnsortedArray() {
90+
// Knot values not sorted.
91+
new ClampedSplineInterpolator().interpolate(new double[] {1, 3, 2 },
92+
new double[] {2, 3, 5},
93+
2, 1);
94+
}
95+
@Test(expected=NumberIsTooSmallException.class)
96+
public void testInsufficientData() {
97+
// Not enough data to interpolate.
98+
new ClampedSplineInterpolator().interpolate(new double[] {1, 2 },
99+
new double[] {2, 3},
100+
2, 1);
101+
}
102+
103+
/**
104+
* Verifies that a clamped spline <i>without</i> specified boundary conditions behaves similar to a natural
105+
* ("unclamped") spline.
106+
*
107+
* <p>
108+
* Using the exponential function <code>e<sup>x</sup></code> over the interval <code>[0, 3]</code>, the test
109+
* evaluates:
110+
* <ol>
111+
* <li>The integral of a clamped spline with specified boundary conditions (endpoint slopes/derivatives).</li>
112+
* <li>The integral of a clamped spline without specified boundary conditions.</li>
113+
* <li>The integral of a natural spline (without boundary conditions by default).</li>
114+
* </ol>
115+
*
116+
* These integrals are compared against the direct integral of the function <code>e<sup>x</sup></code> over the same
117+
* interval.</p>
118+
*
119+
* <p>
120+
* This test is based on "Example 4" in R.L. Burden, J.D. Faires,
121+
* <u>Numerical Analysis</u>, 9th Ed., 2010, Cengage Learning, ISBN 0-538-73351-9, pp 156-157.
122+
* </p>
123+
*/
124+
@Test
125+
public void testIntegral() {
126+
final double[] x = {0, 1, 2, 3};
127+
final double[] y = {1, Math.exp(1), Math.exp(2), Math.exp(3)};
128+
final double fpo = 1;
129+
final double fpn = Math.exp(3);
130+
131+
final ClampedSplineInterpolator clampedSplineInterpolator = new ClampedSplineInterpolator();
132+
final PolynomialSplineFunction clampedSpline = clampedSplineInterpolator.interpolate(x, y, fpo, fpn);
133+
final PolynomialSplineFunction clampedSplineAsNaturalSpline = clampedSplineInterpolator.interpolate(x, y);
134+
135+
final SplineInterpolator naturalSplineInterpolator = new SplineInterpolator();
136+
final PolynomialSplineFunction naturalSpline = naturalSplineInterpolator.interpolate(x, y);
137+
138+
final SimpsonIntegrator integrator = new SimpsonIntegrator();
139+
140+
final double clampedSplineIntegral = integrator.integrate(1000, clampedSpline, 0, 3);
141+
final double clampedSplineAsNaturalSplineIntegral = integrator.integrate(1000, clampedSplineAsNaturalSpline, 0, 3);
142+
final double naturalSplineIntegral = integrator.integrate(1000, naturalSpline, 0, 3);
143+
final double exponentialFunctionIntegral = integrator.integrate(1000, (arg) -> Math.exp(arg), 0, 3);
144+
145+
Assert.assertEquals(Math.abs(clampedSplineAsNaturalSplineIntegral - naturalSplineIntegral), 0, 0);
146+
Assert.assertEquals(Math.abs(exponentialFunctionIntegral - clampedSplineIntegral), 0.02589, 0.1);
147+
Assert.assertEquals(Math.abs(exponentialFunctionIntegral - naturalSplineIntegral), 0.46675, 0.1);
148+
}
149+
150+
/**
151+
* Verifies that f(x[i]) = y[i] for i = 0, ..., n-1 (where n is common length).
152+
*/
153+
private void verifyInterpolation(PolynomialSplineFunction f,
154+
double[] x, double[] y) {
155+
for (int i = 0; i < x.length; i++) {
156+
Assert.assertEquals(f.value(x[i]), y[i], KNOT_TOL);
157+
}
158+
}
159+
160+
/**
161+
* Verifies that interpolating polynomials satisfy consistency requirement: adjacent polynomials must agree through
162+
* two derivatives at knot points.
163+
*/
164+
private void verifyConsistency(PolynomialSplineFunction f,
165+
double[] x) {
166+
PolynomialFunction polynomials[] = f.getPolynomials();
167+
for (int i = 1; i < x.length - 2; i++) {
168+
// evaluate polynomials and derivatives at x[i + 1]
169+
Assert.assertEquals(polynomials[i].value(x[i + 1] - x[i]), polynomials[i + 1].value(0), 0.1);
170+
Assert.assertEquals(polynomials[i].polynomialDerivative().value(x[i + 1] - x[i]),
171+
polynomials[i + 1].polynomialDerivative().value(0), 0.5);
172+
Assert.assertEquals(polynomials[i].polynomialDerivative().polynomialDerivative().value(x[i + 1] - x[i]),
173+
polynomials[i + 1].polynomialDerivative().polynomialDerivative().value(0), 0.5);
174+
}
175+
}
176+
}

‎pom.xml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1054,6 +1054,9 @@ This is avoided by creating an empty directory when svn is not available.
10541054
<contributor>
10551055
<name>Samy Badjoudj</name>
10561056
</contributor>
1057+
<contributor>
1058+
<name>Michael Scholz</name>
1059+
</contributor>
10571060
</contributors>
10581061

10591062
</project>

‎src/changes/changes.xml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,9 @@ Caveat:
9696
to support the whole codebase (it was one of the main reasons for
9797
creating more focused components).
9898
">
99+
<action dev="erans" type="add" issue="MATH-1650" due-to="Michael Scholz">
100+
New feature: "ClampedSplineInterpolator".
101+
</action>
99102
<action dev="erans" type="update" issue="MATH-1669" due-to="Wolff Bock von Wuelfingen">
100103
Javadoc: Fix broken link.
101104
</action>

0 commit comments

Comments
 (0)