Skip to content

Commit f8c867b

Browse files
mtezzelemtezzele
authored
Merge pull request #148 from ndem0/radial_performance
Improve RBF performance
2 parents bbb5147 + 8bc0f63 commit f8c867b

File tree

4 files changed

+62
-92
lines changed

4 files changed

+62
-92
lines changed

‎docs/source/_summaries/pygem.radial.RBF._distance_matrix.rst‎

Lines changed: 0 additions & 6 deletions
This file was deleted.

‎docs/source/radial.rst‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,6 @@ Radial
99
:toctree: _summaries
1010
:nosignatures:
1111

12-
RBF._distance_matrix
1312
RBF._get_weights
1413
RBF.beckert_wendland_c2_basis
1514
RBF.gaussian_spline

‎pygem/radial.py‎

Lines changed: 53 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -32,12 +32,12 @@
3232
:math:`\\gamma_i` is the weight, corresponding to the a-priori selected
3333
:math:`\\mathcal{N}_C` control points, associated to the :math:`i`-th basis
3434
function, and :math:`\\varphi(\\| \\boldsymbol{x} - \\boldsymbol{x_{C_i}}
35-
\\|)` a radial function based on the Euclidean distance between the
36-
control points position :math:`\\boldsymbol{x_{C_i}}` and
37-
:math:`\\boldsymbol{x}`. A radial basis function, generally, is a
38-
real-valued function whose value depends only on the distance from the
39-
origin, so that :math:`\\varphi(\\boldsymbol{x}) = \\tilde{\\varphi}(\\|
40-
\\boldsymbol{x} \\|)`.
35+
\\|)` a radial function based on the Euclidean distance between the control
36+
points position :math:`\\boldsymbol{x_{C_i}}` and :math:`\\boldsymbol{x}`.
37+
A radial basis function, generally, is a real-valued function whose value
38+
depends only on the distance from the origin, so that
39+
:math:`\\varphi(\\boldsymbol{x}) = \\tilde{\\varphi}(\\| \\boldsymbol{x}
40+
\\|)`.
4141
4242
The matrix version of the formula above is:
4343
@@ -142,17 +142,15 @@ def gaussian_spline(X, r):
142142
It implements the following formula:
143143
144144
.. math::
145-
\\varphi(\\| \\boldsymbol{x} \\|) =
146-
e^{-\\frac{\\| \\boldsymbol{x} \\|^2}{r^2}}
145+
\\varphi(\\boldsymbol{x}) = e^{-\\frac{\\boldsymbol{x}^2}{r^2}}
147146
148147
:param numpy.ndarray X: the vector x in the formula above.
149148
:param float r: the parameter r in the formula above.
150149
151150
:return: result: the result of the formula above.
152151
:rtype: float
153152
"""
154-
norm = np.linalg.norm(X)
155-
result = np.exp(-(norm * norm) / (r * r))
153+
result = np.exp(-(X * X) / (r * r))
156154
return result
157155

158156
@staticmethod
@@ -161,17 +159,15 @@ def multi_quadratic_biharmonic_spline(X, r):
161159
It implements the following formula:
162160
163161
.. math::
164-
\\varphi(\\| \\boldsymbol{x} \\|) =
165-
\\sqrt{\\| \\boldsymbol{x} \\|^2 + r^2}
162+
\\varphi(\\boldsymbol{x}) = \\sqrt{\\boldsymbol{x}^2 + r^2}
166163
167164
:param numpy.ndarray X: the vector x in the formula above.
168165
:param float r: the parameter r in the formula above.
169166
170167
:return: result: the result of the formula above.
171168
:rtype: float
172169
"""
173-
norm = np.linalg.norm(X)
174-
result = np.sqrt((norm * norm) + (r * r))
170+
result = np.sqrt((X * X) + (r * r))
175171
return result
176172

177173
@staticmethod
@@ -180,17 +176,16 @@ def inv_multi_quadratic_biharmonic_spline(X, r):
180176
It implements the following formula:
181177
182178
.. math::
183-
\\varphi(\\| \\boldsymbol{x} \\|) =
184-
(\\| \\boldsymbol{x} \\|^2 + r^2 )^{-\\frac{1}{2}}
179+
\\varphi(\\boldsymbol{x}) =
180+
(\\boldsymbol{x}^2 + r^2 )^{-\\frac{1}{2}}
185181
186182
:param numpy.ndarray X: the vector x in the formula above.
187183
:param float r: the parameter r in the formula above.
188184
189185
:return: result: the result of the formula above.
190186
:rtype: float
191187
"""
192-
norm = np.linalg.norm(X)
193-
result = 1.0 / (np.sqrt((norm * norm) + (r * r)))
188+
result = 1.0 / (np.sqrt((X * X) + (r * r)))
194189
return result
195190

196191
@staticmethod
@@ -199,9 +194,9 @@ def thin_plate_spline(X, r):
199194
It implements the following formula:
200195
201196
.. math::
202-
\\varphi(\\| \\boldsymbol{x} \\|) =
203-
\\left\\| \\frac{\\boldsymbol{x} }{r} \\right\\|^2
204-
\\ln \\left\\| \\frac{\\boldsymbol{x} }{r} \\right\\|
197+
\\varphi(\\boldsymbol{x}) =
198+
\\left(\\frac{\\boldsymbol{x}}{r}\\right)^2
199+
\\ln\\frac{\\boldsymbol{x}}{r}
205200
206201
:param numpy.ndarray X: the vector x in the formula above.
207202
:param float r: the parameter r in the formula above.
@@ -210,10 +205,8 @@ def thin_plate_spline(X, r):
210205
:rtype: float
211206
"""
212207
arg = X / r
213-
norm = np.linalg.norm(arg)
214-
result = norm * norm
215-
if norm > 0:
216-
result *= np.log(norm)
208+
result = arg * arg
209+
result = np.where(arg > 0, result * np.log(arg), result)
217210
return result
218211

219212
@staticmethod
@@ -222,21 +215,18 @@ def beckert_wendland_c2_basis(X, r):
222215
It implements the following formula:
223216
224217
.. math::
225-
\\varphi(\\| \\boldsymbol{x} \\|) =
226-
\\left( 1 - \\frac{\\| \\boldsymbol{x} \\|}{r} \\right)^4_+
227-
\\left( 4 \\frac{\\| \\boldsymbol{x} \\|}{r} + 1 \\right)
218+
\\varphi(\\boldsymbol{x}) =
219+
\\left( 1 - \\frac{\\boldsymbol{x}}{r}\\right)^4 +
220+
\\left( 4 \\frac{ \\boldsymbol{x} }{r} + 1 \\right)
228221
229222
:param numpy.ndarray X: the vector x in the formula above.
230223
:param float r: the parameter r in the formula above.
231224
232225
:return: result: the result of the formula above.
233226
:rtype: float
234227
"""
235-
norm = np.linalg.norm(X)
236-
arg = norm / r
237-
first = 0
238-
if (1 - arg) > 0:
239-
first = np.power((1 - arg), 4)
228+
arg = X / r
229+
first = np.where((1 - arg) > 0, np.power((1 - arg), 4), 0)
240230
second = (4 * arg) + 1
241231
result = first * second
242232
return result
@@ -247,18 +237,18 @@ def polyharmonic_spline(self, X, r):
247237
248238
.. math::
249239
250-
\\varphi(\\| \\boldsymbol{x} \\|) =
240+
\\varphi(\\boldsymbol{x}) =
251241
\\begin{cases}
252-
\\|\\frac{\\boldsymbol{x}}{r}\\|^k
242+
\\frac{\\boldsymbol{x}}{r}^k
253243
\\quad & \\text{if}~k = 1,3,5,...\\\\
254-
\\|\\frac{\\boldsymbol{x}}{r}\\|^{k-1}
255-
\\ln(\\|\\frac{\\boldsymbol{x}}{r}\\|^
256-
{\\|\\frac{\\boldsymbol{x}}{r}\\|})
257-
\\quad & \\text{if}~\\|\\frac{\\boldsymbol{x}}{r}\\| < 1,
244+
\\frac{\\boldsymbol{x}}{r}^{k-1}
245+
\\ln(\\frac{\\boldsymbol{x}}{r}^
246+
{\\frac{\\boldsymbol{x}}{r}})
247+
\\quad & \\text{if}~\\frac{\\boldsymbol{x}}{r} < 1,
258248
~k = 2,4,6,...\\\\
259-
\\|\\frac{\\boldsymbol{x}}{r}\\|^k
260-
\\ln(\\|\\frac{\\boldsymbol{x}}{r}\\|)
261-
\\quad & \\text{if}~\\|\\frac{\\boldsymbol{x}}{r}\\| \\ge 1,
249+
\\frac{\\boldsymbol{x}}{r}^k
250+
\\ln(\\frac{\\boldsymbol{x}}{r})
251+
\\quad & \\text{if}~\\frac{\\boldsymbol{x}}{r} \\ge 1,
262252
~k = 2,4,6,...\\\\
263253
\\end{cases}
264254
@@ -270,33 +260,18 @@ def polyharmonic_spline(self, X, r):
270260
"""
271261

272262
k = self.parameters.power
273-
r_sc = np.linalg.norm(X) / r
263+
r_sc = X / r
274264

275265
# k odd
276266
if k & 1:
277267
return np.power(r_sc, k)
278268

269+
print(r_sc)
279270
# k even
280-
if r_sc < 1:
281-
return np.power(r_sc, k - 1) * np.log(np.power(r_sc, r_sc))
282-
else:
283-
return np.power(r_sc, k) * np.log(r_sc)
284-
285-
def _distance_matrix(self, X1, X2):
286-
"""
287-
This private method returns the following matrix:
288-
:math:`\\boldsymbol{D_{ij}} = \\varphi(\\| \\boldsymbol{x_i} -
289-
\\boldsymbol{y_j} \\|)`
290-
291-
:param numpy.ndarray X1: the vector x in the formula above.
292-
:param numpy.ndarray X2: the vector y in the formula above.
293-
294-
:return: matrix: the matrix D.
295-
:rtype: numpy.ndarray
296-
"""
297-
matrix = cdist(X1, X2,
298-
lambda x, y: self.basis(x - y, self.parameters.radius))
299-
return matrix
271+
result = np.where(r_sc < 1,
272+
np.power(r_sc, k - 1) * np.log(np.power(r_sc, r_sc)),
273+
np.power(r_sc, k) * np.log(r_sc))
274+
return result
300275

301276
def _get_weights(self, X, Y):
302277
"""
@@ -315,17 +290,17 @@ def _get_weights(self, X, Y):
315290
:return: weights: the matrix with the weights and the polynomial terms.
316291
:rtype: numpy.matrix
317292
"""
318-
n_points = X.shape[0]
319-
dim = X.shape[1]
320-
identity = np.ones((n_points, 1))
321-
dist = self._distance_matrix(X, X)
322-
H = np.bmat([[dist, identity,
323-
X], [identity.T,
324-
np.zeros((1, 1)),
325-
np.zeros((1, dim))],
326-
[X.T, np.zeros((dim, 1)),
327-
np.zeros((dim, dim))]])
328-
rhs = np.bmat([[Y], [np.zeros((1, dim))], [np.zeros((dim, dim))]])
293+
n_points, dim = X.shape
294+
H = np.zeros((n_points + 3 + 1, n_points + 3 + 1))
295+
H[:n_points, :n_points] = self.basis(
296+
cdist(X, X), self.parameters.radius)
297+
H[n_points, :n_points] = 1.0
298+
H[:n_points, n_points] = 1.0
299+
H[:n_points, -3:] = X
300+
H[-3:, :n_points] = X.T
301+
302+
rhs = np.zeros((n_points + 3 + 1, dim))
303+
rhs[:n_points, :] = Y
329304
weights = np.linalg.solve(H, rhs)
330305
return weights
331306

@@ -335,8 +310,10 @@ def perform(self):
335310
execution it sets `self.modified_mesh_points`.
336311
"""
337312
n_points = self.original_mesh_points.shape[0]
338-
dist = self._distance_matrix(self.original_mesh_points,
339-
self.parameters.original_control_points)
313+
dist = self.basis(
314+
cdist(self.original_mesh_points,
315+
self.parameters.original_control_points),
316+
self.parameters.radius)
340317
identity = np.ones((n_points, 1))
341318
H = np.bmat([[dist, identity, self.original_mesh_points]])
342319
self.modified_mesh_points = np.asarray(np.dot(H, self.weights))

‎tests/test_radial.py‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -84,7 +84,7 @@ def test_gaussian_spline(self):
8484
filename='tests/test_datasets/parameters_rbf_default.prm')
8585
rbf = rad.RBF(params, self.get_cube_mesh_points())
8686
value = rbf.gaussian_spline(
87-
np.array([0.5, 1, 2, 0.2]).reshape(4, 1), 0.2)
87+
np.linalg.norm(np.array([0.5, 1, 2, 0.2])), 0.2)
8888
np.testing.assert_almost_equal(value, 0.0)
8989

9090
def test_multi_quadratic_biharmonic_spline(self):
@@ -93,7 +93,7 @@ def test_multi_quadratic_biharmonic_spline(self):
9393
filename='tests/test_datasets/parameters_rbf_default.prm')
9494
rbf = rad.RBF(params, self.get_cube_mesh_points())
9595
value = rbf.multi_quadratic_biharmonic_spline(
96-
np.array([0.5, 1, 2, 0.2]).reshape(4, 1), 0.2)
96+
np.linalg.norm(np.array([0.5, 1, 2, 0.2])), 0.2)
9797
np.testing.assert_almost_equal(value, 2.30867927612)
9898

9999
def test_inv_multi_quadratic_biharmonic_spline(self):
@@ -102,7 +102,7 @@ def test_inv_multi_quadratic_biharmonic_spline(self):
102102
filename='tests/test_datasets/parameters_rbf_default.prm')
103103
rbf = rad.RBF(params, self.get_cube_mesh_points())
104104
value = rbf.inv_multi_quadratic_biharmonic_spline(
105-
np.array([0.5, 1, 2, 0.2]).reshape(4, 1), 0.2)
105+
np.linalg.norm(np.array([0.5, 1, 2, 0.2])), 0.2)
106106
np.testing.assert_almost_equal(value, 0.433148081824)
107107

108108
def test_thin_plate_spline(self):
@@ -111,7 +111,7 @@ def test_thin_plate_spline(self):
111111
filename='tests/test_datasets/parameters_rbf_default.prm')
112112
rbf = rad.RBF(params, self.get_cube_mesh_points())
113113
value = rbf.thin_plate_spline(
114-
np.array([0.5, 1, 2, 0.2]).reshape(4, 1), 0.2)
114+
np.linalg.norm(np.array([0.5, 1, 2, 0.2])), 0.2)
115115
np.testing.assert_almost_equal(value, 323.000395428)
116116

117117
def test_beckert_wendland_c2_basis_01(self):
@@ -120,7 +120,7 @@ def test_beckert_wendland_c2_basis_01(self):
120120
filename='tests/test_datasets/parameters_rbf_default.prm')
121121
rbf = rad.RBF(params, self.get_cube_mesh_points())
122122
value = rbf.beckert_wendland_c2_basis(
123-
np.array([0.5, 1, 2, 0.2]).reshape(4, 1), 0.2)
123+
np.linalg.norm(np.array([0.5, 1, 2, 0.2])), 0.2)
124124
np.testing.assert_almost_equal(value, 0.0)
125125

126126
def test_beckert_wendland_c2_basis_02(self):
@@ -129,7 +129,7 @@ def test_beckert_wendland_c2_basis_02(self):
129129
filename='tests/test_datasets/parameters_rbf_default.prm')
130130
rbf = rad.RBF(params, self.get_cube_mesh_points())
131131
value = rbf.beckert_wendland_c2_basis(
132-
np.array([0.1, 0.15, -0.2]).reshape(3, 1), 0.9)
132+
np.linalg.norm(np.array([0.1, 0.15, -0.2])), 0.9)
133133
np.testing.assert_almost_equal(value, 0.529916819595)
134134

135135
def test_polyharmonic_spline_k_even(self):
@@ -139,7 +139,7 @@ def test_polyharmonic_spline_k_even(self):
139139
params.power = 3
140140
rbf = rad.RBF(params, self.get_cube_mesh_points())
141141
value = rbf.polyharmonic_spline(
142-
np.array([0.1, 0.15, -0.2]).reshape(3, 1), 0.9)
142+
np.linalg.norm(np.array([0.1, 0.15, -0.2])), 0.9)
143143
np.testing.assert_almost_equal(value, 0.02677808)
144144

145145
def test_polyharmonic_spline_k_odd1(self):
@@ -149,7 +149,7 @@ def test_polyharmonic_spline_k_odd1(self):
149149
params.power = 2
150150
rbf = rad.RBF(params, self.get_cube_mesh_points())
151151
value = rbf.polyharmonic_spline(
152-
np.array([0.1, 0.15, -0.2]).reshape(3, 1), 0.9)
152+
np.linalg.norm(np.array([0.1, 0.15, -0.2])), 0.9)
153153
np.testing.assert_almost_equal(value, -0.1080092)
154154

155155
def test_polyharmonic_spline_k_odd2(self):
@@ -159,5 +159,5 @@ def test_polyharmonic_spline_k_odd2(self):
159159
params.power = 2
160160
rbf = rad.RBF(params, self.get_cube_mesh_points())
161161
value = rbf.polyharmonic_spline(
162-
np.array([0.1, 0.15, -0.2]).reshape(3, 1), 0.2)
162+
np.linalg.norm(np.array([0.1, 0.15, -0.2])), 0.2)
163163
np.testing.assert_almost_equal(value, 0.53895331)

0 commit comments

Comments
 (0)