1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
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 com.irurueta.numerical.interpolation;
17
18 /**
19 * Computes polynomial interpolation.
20 * For large sets of ata, this interpolator might return inaccurate results.
21 * Additionally, accuracy worsens as the polynomial degree to interpolate increases.
22 */
23 public class PolynomialInterpolator extends BaseInterpolator {
24
25 /**
26 * An indication of interpolation error on the y values of the last call to
27 * {@link #interpolate(double)}.
28 */
29 private double dy;
30
31 /**
32 * Constructor.
33 *
34 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
35 * decreasing)
36 * @param y y values to interpolate to.
37 * @param m length of x's and y's to take into account. Must be less or equal than x or y
38 * length.
39 * @throws IllegalArgumentException if x or y have invalid length or m exceeds length of x or y.
40 */
41 public PolynomialInterpolator(final double[] x, final double[] y, final int m) {
42 this(x, y, m, true);
43 }
44
45 /**
46 * Constructor.
47 *
48 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
49 * decreasing)
50 * @param y y values to interpolate to.
51 * @param m length of x's and y's to take into account. Must be less or equal than x or y
52 * length.
53 * @param check true to make validations, false otherwise.
54 * @throws IllegalArgumentException if x or y have invalid length or m exceeds length of x or y.
55 */
56 public PolynomialInterpolator(final double[] x, final double[] y, final int m, final boolean check) {
57 super(x, y, m, check);
58 dy = 0.0;
59 }
60
61 /**
62 * Constructor.
63 *
64 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
65 * decreasing)
66 * @param y y values to interpolate to.
67 * @throws IllegalArgumentException if x or y have invalid length or m exceeds length of x or y.
68 */
69 public PolynomialInterpolator(final double[] x, final double[] y) {
70 this(x, y, x.length);
71 }
72
73 /**
74 * Gets an indication of the error of interpolation on the y values.
75 *
76 * @return indication of error of interpolation.
77 */
78 public double getDy() {
79 return dy;
80 }
81
82 /**
83 * Actual interpolation method.
84 *
85 * @param jl index where value x to be interpolated in located in the array of xx.
86 * @param x value to obtain interpolation for.
87 * @return interpolated value.
88 * @throws InterpolationException if interpolation fails.
89 */
90 @SuppressWarnings("Duplicates")
91 @Override
92 public double rawinterp(final int jl, final double x) throws InterpolationException {
93 int i;
94 int m;
95 int ns = 0;
96 double y;
97 double den;
98 double dif;
99 double dift;
100 double ho;
101 double hp;
102 double w;
103 final var xa = xx;
104 final var ya = yy;
105 final var c = new double[mm];
106 final var d = new double[mm];
107 dif = Math.abs(x - xa[jl]);
108 for (i = 0; i < mm; i++) {
109 // Here we find the index ns of the closest table entry
110 if ((dift = Math.abs(x - xa[jl + i])) < dif) {
111 ns = i;
112 dif = dift;
113 }
114 // and initialize the tableau of c's and d's
115 c[i] = ya[jl + i];
116 d[i] = ya[jl + i];
117 }
118 // This is the initial approximation to y
119 y = ya[jl + ns--];
120 for (m = 1; m < mm; m++) {
121 // For each column of the tableau
122 for (i = 0; i < mm - m; i++) {
123 // we loop over the current c's and d's and update them
124 ho = xa[jl + i] - x;
125 hp = xa[jl + i + m] - x;
126 w = c[i + 1] - d[i];
127 den = ho - hp;
128 if (den == 0.0) {
129 // This error can occur only if two input xa's are (to within rounding error)
130 // identical
131 throw new InterpolationException();
132 }
133 den = w / den;
134 // Here the c’s and d’s are updated.
135 d[i] = hp * den;
136 c[i] = ho * den;
137 }
138 dy = 2 * (ns + 1) < (mm - m) ? c[ns + 1] : d[ns--];
139 y += dy;
140 // After each column in the tableau is completed, we decide which correction, c or d,
141 // we want to add to our accumulating value of y, i.e., which path to take through the
142 // tableau — forking up or down. We do this in such a way as to take the most “straight
143 // line” route through the tableau to its apex, updating ns accordingly to keep track
144 // of where we are. This route keeps the partial approximations centered (insofar as
145 // possible) on the target x. The last dy added is thus the error indication.
146 }
147 return y;
148 }
149 }