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 import com.irurueta.algebra.AlgebraException;
19 import com.irurueta.algebra.LUDecomposer;
20 import com.irurueta.algebra.Matrix;
21 import com.irurueta.algebra.WrongSizeException;
22
23 /**
24 * Interpolates sparsely defined points using D.G. Krige method.
25 */
26 public class KrigingInterpolator {
27 /**
28 * Data to compute interpolations from.
29 * Each row corresponds to a point.
30 * The number of columns determines the number of dimensions of provided points.
31 */
32 private final Matrix x;
33
34 /**
35 * Variogram of provided data.
36 */
37 private final Variogram vgram;
38
39 /**
40 * Number of dimensions of each point.
41 */
42 private final int ndim;
43
44 /**
45 * Number of provided points.
46 */
47 private final int npt;
48
49 /**
50 * Most recently computed value.
51 */
52 private double lastval;
53
54 /**
55 * Most recently computed error.
56 */
57 private double lasterr;
58
59 private final Matrix dstar;
60
61 private final Matrix vstar;
62
63 private final Matrix yvi;
64
65 /**
66 * LU decomposer.
67 */
68 private final LUDecomposer vi;
69
70 /**
71 * Values of i-th point.
72 */
73 private final double[] xi;
74
75 /**
76 * Constructor.
77 *
78 * @param xx Data to compute interpolations from. Each row corresponds to a point. The
79 * number of columns determines the number of dimensions of provided points.
80 * @param yy Function values for each provided point.
81 * @param vargram Variogram to use.
82 * @param err array containing estimated error of function values.
83 * @throws InterpolationException if initialization fails for numerical reasons.
84 */
85 public KrigingInterpolator(
86 final Matrix xx, final double[] yy, final Variogram vargram, final double[] err)
87 throws InterpolationException {
88 try {
89 x = xx;
90 vgram = vargram;
91 npt = xx.getRows();
92 ndim = xx.getColumns();
93 final var nptPlus1 = npt + 1;
94 dstar = new Matrix(nptPlus1, 1);
95 vstar = new Matrix(nptPlus1, 1);
96 final var v = new Matrix(nptPlus1, nptPlus1);
97 final var y = new Matrix(nptPlus1, 1);
98 yvi = new Matrix(nptPlus1, 1);
99
100 xi = new double[ndim];
101 final var xj = new double[ndim];
102 final var end = ndim - 1;
103
104 int i;
105 int j;
106 // Fill Y and V
107 for (i = 0; i < npt; i++) {
108 y.setElementAtIndex(i, yy[1]);
109 for (j = i; j < npt; j++) {
110 x.getSubmatrixAsArray(i, 0, i, end, xi);
111 x.getSubmatrixAsArray(j, 0, j, end, xj);
112 final var evaluation = vgram.evaluate(rdist(xi, xj));
113 v.setElementAt(i, j, evaluation);
114 v.setElementAt(j, i, evaluation);
115 }
116
117 v.setElementAt(i, npt, 1.0);
118 v.setElementAt(npt, i, 1.0);
119 }
120 v.setElementAt(npt, npt, 0.0);
121 y.setElementAtIndex(npt, 0.0);
122
123 if (err != null) {
124 for (i = 0; i < npt; i++) {
125 v.setElementAt(i, i, v.getElementAt(i, i) - sqr(err[i]));
126 }
127 }
128
129 vi = new LUDecomposer(v);
130 vi.decompose();
131 vi.solve(y, yvi);
132 } catch (final AlgebraException e) {
133 throw new InterpolationException(e);
134 }
135 }
136
137 /**
138 * Constructor.
139 *
140 * @param xx Data to compute interpolations from. Each row corresponds to a point. The
141 * number of columns determines the number of dimensions of provided points.
142 * @param yy Function values for each provided point.
143 * @param vargram Variogram to use.
144 * @throws InterpolationException if initialization fails for numerical reasons.
145 */
146 public KrigingInterpolator(final Matrix xx, final double[] yy, final Variogram vargram)
147 throws InterpolationException {
148 this(xx, yy, vargram, null);
149 }
150
151 /**
152 * Constructor.
153 *
154 * @param xx Data to compute interpolations from. Each row corresponds to a point. The
155 * number of columns determines the number of dimensions of provided points.
156 * @param yy Function values for each provided point.
157 * @throws InterpolationException if initialization fails for numerical reasons.
158 */
159 public KrigingInterpolator(final Matrix xx, final double[] yy) throws InterpolationException {
160 this(xx, yy, new Variogram(xx, yy));
161 }
162
163 /**
164 * Gets number of dimensions of each point.
165 *
166 * @return number of dimensions of each point.
167 */
168 public int getNdim() {
169 return ndim;
170 }
171
172 /**
173 * Gets number of provided points.
174 *
175 * @return number of provided points.
176 */
177 public int getNpt() {
178 return npt;
179 }
180
181 /**
182 * Gets most recently computed interpolated value.
183 *
184 * @return most recently computed interpolatd value.
185 */
186 public double getLastVal() {
187 return lastval;
188 }
189
190 /**
191 * Gets most recently computed error.
192 *
193 * @return most recently computed error.
194 */
195 public double getLastErr() {
196 return lasterr;
197 }
198
199 /**
200 * Returns interpolated value at provided point.
201 *
202 * @param xstar point where interpolation is computed.
203 * @return computed interpolation.
204 * @throws InterpolationException if interpolation fails.
205 */
206 public double interpolate(final double[] xstar) throws InterpolationException {
207 try {
208 int i;
209 final var end = ndim - 1;
210 for (i = 0; i < npt; i++) {
211 x.getSubmatrixAsArray(i, 0, i, end, xi);
212 vstar.setElementAtIndex(i, vgram.evaluate(rdist(xstar, xi)));
213 }
214
215 vstar.setElementAtIndex(npt, 1.0);
216 lastval = 0.;
217 for (i = 0; i <= npt; i++) {
218 lastval += yvi.getElementAtIndex(i) * vstar.getElementAtIndex(i);
219 }
220 return lastval;
221
222 } catch (WrongSizeException e) {
223 throw new InterpolationException(e);
224 }
225 }
226
227 /**
228 * Returns interpolated value at provided point, and estimated error.
229 *
230 * @param xstar point where interpolation is computed.
231 * @param esterr array where estimated error will be stored on its first position.
232 * @return computed interpolation.
233 * @throws InterpolationException if interpolation fails.
234 */
235 public double interpolate(final double[] xstar, final double[] esterr) throws InterpolationException {
236 try {
237 lastval = interpolate(xstar);
238 vi.solve(vstar, dstar);
239 lasterr = 0;
240 for (var i = 0; i <= npt; i++) {
241 lasterr += dstar.getElementAtIndex(i) * vstar.getElementAtIndex(i);
242 }
243
244 lasterr = Math.sqrt(Math.max(0.0, lasterr));
245 esterr[0] = lasterr;
246 return lastval;
247 } catch (final AlgebraException e) {
248 throw new InterpolationException(e);
249 }
250 }
251
252 /**
253 * Computes euclidean distance between two points.
254 *
255 * @param x1 1st point.
256 * @param x2 2nd point.
257 * @return euclidean distance.
258 */
259 private double rdist(final double[] x1, final double[] x2) {
260 var d = 0.0;
261 for (var i = 0; i < ndim; i++) {
262 d += sqr(x1[i] - x2[i]);
263 }
264 return Math.sqrt(d);
265 }
266
267 /**
268 * Computes squared value.
269 *
270 * @param value a value.
271 * @return computed square value.
272 */
273 private static double sqr(final double value) {
274 return value * value;
275 }
276
277 /**
278 * Variogram function.
279 * Follows expression: v(r) = alpha * r^beta
280 * where beta is considered fixed and alpha is fitted by unweighted least squares over all pairs
281 * of data points i and j.
282 * The value of beta should be in the range 1 <= beta < 2.
283 * A good general choice is 1.5, but for functions with a strong linear trend, you may want to
284 * experiment with values as large as 1.99 (The value 2 gives a degenerate matrix and meaningless
285 * results).
286 */
287 public static class Variogram {
288
289 /**
290 * Default Beta value to be used.
291 */
292 public static final double DEFAULT_BETA = 1.5;
293
294 /**
295 * Default offset to use.
296 */
297 public static final double DEFAULT_NUG = 0.0;
298
299 /**
300 * Estimated alpha of variogram.
301 */
302 private final double alph;
303
304 /**
305 * Beta of variogram.
306 */
307 private final double bet;
308
309 /**
310 * Squared offset.
311 */
312 private final double nugsq;
313
314 /**
315 * Constructor.
316 *
317 * @param x Data to compute interpolations from. Each row corresponds to a point. The
318 * number of columns determines the number of dimensions of provided points.
319 * @param y Function values for each provided point.
320 * @param beta Beta to be used for variogram. The value of beta should be in the range
321 * 1 <= beta < 2.
322 * @param nug Offset of variogram.
323 */
324 public Variogram(final Matrix x, final double[] y, final double beta, final double nug) {
325 bet = beta;
326 nugsq = nug * nug;
327
328 int i;
329 int j;
330 int k;
331 final var npt = x.getRows();
332 final var ndim = x.getColumns();
333 double rb;
334 var num = 0.0;
335 var denom = 0.0;
336 for (i = 0; i < npt; i++) {
337 for (j = i + 1; j < npt; j++) {
338 rb = 0.;
339 for (k = 0; k < ndim; k++) {
340 rb += sqr(x.getElementAt(i, k) - x.getElementAt(j, k));
341 }
342 rb = Math.pow(rb, 0.5 * beta);
343 num += rb * (0.5 * sqr(y[i] - y[j]) - nugsq);
344 denom += sqr(rb);
345 }
346 }
347 alph = num / denom;
348 }
349
350 /**
351 * Constructor using default zero offset.
352 *
353 * @param x Data to compute interpolations from. Each row corresponds to a point. The
354 * number of columns determines the number of dimensions of provided points.
355 * @param y Function values for each provided point.
356 * @param beta Beta to be used for variogram. The value of beta should be in the range
357 * 1 <= beta < 2.
358 */
359 public Variogram(final Matrix x, final double[] y, final double beta) {
360 this(x, y, beta, DEFAULT_NUG);
361 }
362
363 /**
364 * Constructor using default beta = 1.5 and zero offset.
365 *
366 * @param x Data to compute interpolations from. Each row corresponds to a point. The
367 * number of columns determines the number of dimensions of provided points.
368 * @param y Function values for each provided point.
369 */
370 public Variogram(final Matrix x, final double[] y) {
371 this(x, y, DEFAULT_BETA);
372 }
373
374 /**
375 * Evaluates variogram at provided distance.
376 *
377 * @param r a distance.
378 * @return variogram evaluation.
379 */
380 public double evaluate(final double r) {
381 return nugsq + alph * Math.pow(r, bet);
382 }
383 }
384 }