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 * Abstract base class used by all interpolation implementations.
20 */
21 public abstract class BaseInterpolator {
22 /**
23 * Length of x and y values to be interpolated.
24 */
25 protected final int n;
26
27 /**
28 * Length of data to be taken into account on x's and y's.
29 */
30 protected final int mm;
31
32 protected int cor = 0;
33
34 private int jsav = 0;
35
36 private final int dj;
37
38 /**
39 * X values to be used for interpolation estimation. Must be monotonic (either increasing or
40 * decreasing).
41 */
42 protected final double[] xx;
43
44 /**
45 * Y values to be used for interpolation estimation.
46 */
47 protected final double[] yy;
48
49 /**
50 * Constructor.
51 *
52 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
53 * decreasing)
54 * @param y y values to interpolate to.
55 * @param m length of x's and y's to take into account. Must be less or equal than x or y
56 * length.
57 * @param check true to make validations, false otherwise.
58 * @throws IllegalArgumentException if x or y have invalid length or m exceeds length of x or y.
59 */
60 BaseInterpolator(final double[] x, final double[] y, final int m, final boolean check) {
61 n = x.length;
62 mm = m;
63
64 if (check) {
65 if (n != y.length) {
66 throw new IllegalArgumentException("mismatched x and y length");
67 }
68 if (n < 2 || mm < 2) {
69 throw new IllegalArgumentException("x length is too small");
70 }
71 if (mm > n) {
72 throw new IllegalArgumentException("m exceeds length of x or y");
73 }
74 }
75
76 xx = x;
77 yy = y;
78
79 dj = Math.min(1, (int) Math.pow(n, 0.25));
80 }
81
82 /**
83 * Constructor.
84 *
85 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
86 * decreasing)
87 * @param y y values to interpolate to.
88 * @param m length of x's and y's to take into account. Must be less or equal than x or y
89 * length.
90 * @throws IllegalArgumentException if x or y have invalid length or m exceeds length of x or y.
91 */
92 BaseInterpolator(final double[] x, final double[] y, final int m) {
93 this(x, y, m, true);
94 }
95
96 /**
97 * Given a value x, returns an interpolated value, using data pointed to by {@link #xx} and
98 * {@link #yy}.
99 *
100 * @param x value to obtain interpolation for.
101 * @return interpolated value.
102 * @throws InterpolationException if interpolation fails.
103 */
104 public double interpolate(final double x) throws InterpolationException {
105 final var jlo = cor != 0 ? hunt(x) : locate(x);
106 return rawinterp(jlo, x);
107 }
108
109 /**
110 * Actual interpolation method to be implemented by subclasses.
111 *
112 * @param jlo index where value x to be interpolated in located in the array of xx.
113 * @param x value to obtain interpolation for.
114 * @return interpolated value.
115 * @throws InterpolationException if interpolation fails.
116 */
117 public abstract double rawinterp(final int jlo, final double x) throws InterpolationException;
118
119 /**
120 * Given a value x, returns a value j such that x is (insofar as possible) centered in the
121 * subrange xx[j..j+mm-1], where xx is the stored array. The value in xx must be monotonic,
122 * either increasing or decreasing. The returned value is not less than 0, nor greater than
123 * n - 1.
124 * @param x value to obtain interpolation for.
125 * @return position where value to obtain interpolation lies in the array of {@link #xx}.
126 */
127 @SuppressWarnings("Duplicates")
128 protected int locate(final double x) {
129 int ju;
130 int jm;
131 int jl;
132 // True if ascending order of table, false otherwise.
133 final var ascend = (xx[n - 1] >= xx[0]);
134
135 // Initialize lower and upper limits.
136 jl = 0;
137 ju = n - 1;
138
139 while (ju - jl > 1) {
140 // If we are not yet done, compute a midpoint
141 jm = (ju + jl) >> 1;
142 if (x >= xx[jm] == ascend) {
143 // and replace either the lower limit
144 jl = jm;
145 } else {
146 // or the upper limit, as appropriate
147 ju = jm;
148 }
149
150 // Repeat until the test condition is satisfied
151 }
152
153 // Decide whether to use hunt or locate next time
154 cor = Math.abs(jl - jsav) > dj ? 0 : 1;
155 jsav = jl;
156 return Math.max(0, Math.min(n - mm, jl - ((mm - 2) >> 1)));
157 }
158
159 /**
160 * Given a value x, returns a value j such that x is (insofar as possible) centered in the
161 * subrange xx[j..j+mm-1], where xx is the stored array. The value in xx must be monotonic,
162 * either increasing or decreasing. The returned value is not less than 0, nor greater than
163 * n - 1.
164 * This method starts with a guessed position in the table. It first "hunts" either up or own,
165 * in increments of 1, then 2, then 3, etc. until the desired value is bracketed. It then
166 * bisects in the bracketed interval. At worst, this routine is about a factor of 2 slower than
167 * {@link #locate(double)} (if the hunt phase expands to include the whole table). At best, it
168 * can be a factor of log(n)/log(2) faster than {@link #locate(double)} if the desired point is
169 * usually quite close to the input guess.
170 *
171 * @param x value to obtain interpolation for.
172 * @return position where value to obtain interpolation lies in the array of {@link #xx}.
173 */
174 @SuppressWarnings("Duplicates")
175 protected int hunt(final double x) {
176 var jl = jsav;
177 int jm;
178 int ju;
179 var inc = 1;
180
181 // Ture if ascending order of table, false otherwise
182 final var ascnd = (xx[n - 1] >= xx[0]);
183 if (jl < 0 || jl > n - 1) {
184 // Input guess not useful. Go immediately to bisection
185 jl = 0;
186 ju = n - 1;
187 } else {
188 if (x >= xx[jl] == ascnd) {
189 // Hunt up
190 for (; ; ) {
191 ju = jl + inc;
192 if (ju >= n - 1) {
193 // Off end of table
194 ju = n - 1;
195 break;
196 } else if (x < xx[ju] == ascnd) {
197 // Found bracket
198 break;
199 } else {
200 // Not done, so double the increment and try again
201 jl = ju;
202 inc += inc;
203 }
204 }
205 } else {
206 // Hunt down
207 ju = jl;
208 for (; ; ) {
209 jl = jl - inc;
210 if (jl <= 0) {
211 // Off end of table
212 jl = 0;
213 break;
214 } else if (x >= xx[jl] == ascnd) {
215 // Found bracket
216 break;
217 } else {
218 // Not done, so double the increment and try again
219 ju = jl;
220 inc += inc;
221 }
222 }
223 }
224 }
225
226 // Hunt is done, so begin the final bisection phase
227 while (ju - jl > 1) {
228 jm = (ju + jl) >> 1;
229 if (x >= xx[jm] == ascnd) {
230 jl = jm;
231 } else {
232 ju = jm;
233 }
234 }
235
236 // Decide whether to use hunt or locate next time
237 cor = Math.abs(jl - jsav) > dj ? 0 : 1;
238 jsav = jl;
239 return Math.max(0, Math.min(n - mm, jl - ((mm - 2) >> 1)));
240 }
241 }