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 of dimension "dim" using a Radial Basis Function.
25 */
26 public class RadialBasisFunctionInterpolator extends BaseRadialBasisFunctionInterpolator {
27 /**
28 * Computed weights to compute interpolation from provided points.
29 */
30 private final Matrix w;
31
32 /**
33 * Radial basis function defining a value based on the distance of two points.
34 */
35 private final RadialBasisFunction fn;
36
37 /**
38 * Indicates whether normalized Radial Basis Function (RBF) must be used or not.
39 */
40 private final boolean norm;
41
42 /**
43 * Constructor.
44 *
45 * @param ptss Matrix containing points to interpolate from. Each row contains one point.
46 * Matrix will have n points (rows) having a dimension (columns) equal to dim.
47 * @param valss values of function at provided points. Must have the same length as the number
48 * of rows of provided points matrix.
49 * @param func function to be used as Radial Basis Function (RBF).
50 * @param nrbf true to normalize RBF, false otherwise.
51 * @throws InterpolationException if provided points are redundant and result in a degenerate
52 * solution.
53 * @throws IllegalArgumentException if provided values array does not match the number of
54 * points (rows) in provided matrix.
55 */
56 public RadialBasisFunctionInterpolator(
57 final Matrix ptss, final double[] valss, final RadialBasisFunction func, final boolean nrbf)
58 throws InterpolationException {
59 super(ptss);
60
61 if (valss.length != n) {
62 throw new IllegalArgumentException("wrong length of values");
63 }
64
65 try {
66 w = new Matrix(n, 1);
67 fn = func;
68 norm = nrbf;
69
70 final var pj = new double[dim];
71
72 int i;
73 int j;
74 double sum;
75 final var rbf = new Matrix(n, n);
76 final var rhs = new Matrix(n, 1);
77 for (i = 0; i < n; i++) {
78 // Fill the matrix phi(|ri - rj|) and the right hand sisde (rhs) vector
79 sum = 0.;
80 for (j = 0; j < n; j++) {
81 final var endCol = dim - 1;
82 pts.getSubmatrixAsArray(i, 0, i, endCol, pi);
83 pts.getSubmatrixAsArray(j, 0, j, endCol, pj);
84 final var value = fn.evaluate(rad(pi, pj));
85 rbf.setElementAt(i, j, value);
86 sum += value;
87 }
88
89 if (norm) {
90 rhs.setElementAtIndex(i, sum * valss[i]);
91 } else {
92 rhs.setElementAtIndex(i, valss[i]);
93 }
94 }
95
96 // Solve the set of linear equations
97 final var lu = new LUDecomposer(rbf);
98 lu.decompose();
99 lu.solve(rhs, w);
100 } catch (final AlgebraException e) {
101 throw new InterpolationException(e);
102 }
103 }
104
105 /**
106 * Constructor.
107 *
108 * @param ptss Matrix containing points to interpolate from. Each row contains one point.
109 * Matrix will have n points (rows) having a dimension (columns) equal to dim.
110 * @param valss values of function at provided points. Must have the same length as the number
111 * of rows of provided points matrix.
112 * @param func function to be used as Radial Basis Function (RBF).
113 * @throws InterpolationException if provided points are redundant and result in a degenerate
114 * solution.
115 * @throws IllegalArgumentException if provided values array does not match the number of
116 * points (rows) in provided matrix.
117 */
118 public RadialBasisFunctionInterpolator(
119 final Matrix ptss, final double[] valss, final RadialBasisFunction func) throws InterpolationException {
120 this(ptss, valss, func, false);
121 }
122
123 /**
124 * Returns the interpolated function value at a dim-dimensional point pt.
125 *
126 * @param pt dim-dimensional point where interpolation must be computed.
127 * @return result of interpolation.
128 * @throws IllegalArgumentException if provided point has an invalid length.
129 */
130 @Override
131 public double interpolate(final double[] pt) {
132 if (pt.length != dim) {
133 throw new IllegalArgumentException("Wrong point length");
134 }
135
136 double fval;
137 var sum = 0.0;
138 var sumw = 0.0;
139
140 try {
141 for (var i = 0; i < n; i++) {
142 final var endCol = dim - 1;
143 pts.getSubmatrixAsArray(i, 0, i, endCol, pi);
144 fval = fn.evaluate(rad(pt, pi));
145 sumw += w.getElementAtIndex(i) * fval;
146 sum += fval;
147 }
148 } catch (final WrongSizeException ignore) {
149 // never happens
150 }
151
152 return norm ? sumw / sum : sumw;
153 }
154 }