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.Matrix;
19 import com.irurueta.algebra.WrongSizeException;
20
21 /**
22 * Interpolates sparsely defined points of dimension "dim" using Shepard interpolation, which is
23 * a simplification of Radial Basis Function interpolation that achieves less accurate results but
24 * having less computational cost.
25 */
26 public class ShepardInterpolator extends BaseRadialBasisFunctionInterpolator {
27 /**
28 * Values of function at provided points. Must have the same length as the number of rows of
29 * provided points matrix.
30 */
31 private final double[] vals;
32
33 /**
34 * Negative value of p parameter controlling Shepard power-law function phi(r) = r^-p.
35 * Typically, p lies between 1 and 3.
36 */
37 private final double pneg;
38
39 /**
40 * Constructor.
41 *
42 * @param ptss Matrix containing points to interpolate from. Each row contains one point.
43 * Matrix will have n points (rows) having a dimension (columns) equal to dim.
44 * @param valss values of function at provided points. Must have the same length as the number
45 * of rows of provided points matrix.
46 * @param p p parameter controlling Shepard power-law function phi(r) = r^-p. Typically, p
47 * lies between 1 and 3.
48 * @throws IllegalArgumentException if provided values array does not match the number of
49 * points (rows) in provided matrix.
50 */
51 public ShepardInterpolator(final Matrix ptss, final double[] valss, final double p) {
52 super(ptss);
53
54 if (valss.length != n) {
55 throw new IllegalArgumentException("wrong length of values");
56 }
57
58 vals = valss;
59 pneg = -p;
60 }
61
62 /**
63 * Constructor using default p parameter, which is equal to 2.0.
64 *
65 * @param ptss Matrix containing points to interpolate from. Each row contains one point.
66 * Matrix will have n points (rows) having a dimension (columns) equal to dim.
67 * @param valss values of function at provided points. Must have the same length as the number
68 * of rows of provided points matrix.
69 * @throws IllegalArgumentException if provided values array does not match the number of
70 * points (rows) in provided matrix.
71 */
72 public ShepardInterpolator(final Matrix ptss, final double[] valss) {
73 this(ptss, valss, 2.0);
74 }
75
76 /**
77 * Returns the interpolated function value at a dim-dimensional point pt.
78 *
79 * @param pt dim-dimensional point where interpolation must be computed.
80 * @return result of interpolation.
81 * @throws IllegalArgumentException if provided point has an invalid length.
82 */
83 public double interpolate(final double[] pt) {
84 double r;
85 double w;
86 var sum = 0.0;
87 var sumw = 0.0;
88 if (pt.length != dim) {
89 throw new IllegalArgumentException("Wrong point length");
90 }
91
92 try {
93 for (var i = 0; i < n; i++) {
94 final var endCol = dim - 1;
95 pts.getSubmatrixAsArray(i, 0, i, endCol, pi);
96 if ((r = rad(pt, pi)) == 0.0) return vals[i];
97 sum += (w = Math.pow(r, pneg));
98 sumw += w * vals[i];
99 }
100 } catch (final WrongSizeException ignore) {
101 // never happens
102 }
103
104 return sumw / sum;
105 }
106 }