Skip to content

Improve RBF performance #148

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Apr 20, 2018
Merged

Improve RBF performance #148

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 0 additions & 6 deletions docs/source/_summaries/pygem.radial.RBF._distance_matrix.rst
View file
Open in desktop

This file was deleted.

1 change: 0 additions & 1 deletion docs/source/radial.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,6 @@ Radial
:toctree: _summaries
:nosignatures:

RBF._distance_matrix
RBF._get_weights
RBF.beckert_wendland_c2_basis
RBF.gaussian_spline
Expand Down
129 changes: 53 additions & 76 deletions pygem/radial.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,12 +32,12 @@
:math:`\\gamma_i` is the weight, corresponding to the a-priori selected
:math:`\\mathcal{N}_C` control points, associated to the :math:`i`-th basis
function, and :math:`\\varphi(\\| \\boldsymbol{x} - \\boldsymbol{x_{C_i}}
\\|)` a radial function based on the Euclidean distance between the
control points position :math:`\\boldsymbol{x_{C_i}}` and
:math:`\\boldsymbol{x}`. A radial basis function, generally, is a
real-valued function whose value depends only on the distance from the
origin, so that :math:`\\varphi(\\boldsymbol{x}) = \\tilde{\\varphi}(\\|
\\boldsymbol{x} \\|)`.
\\|)` a radial function based on the Euclidean distance between the control
points position :math:`\\boldsymbol{x_{C_i}}` and :math:`\\boldsymbol{x}`.
A radial basis function, generally, is a real-valued function whose value
depends only on the distance from the origin, so that
:math:`\\varphi(\\boldsymbol{x}) = \\tilde{\\varphi}(\\| \\boldsymbol{x}
\\|)`.

The matrix version of the formula above is:

Expand Down Expand Up @@ -142,17 +142,15 @@ def gaussian_spline(X, r):
It implements the following formula:

.. math::
\\varphi(\\| \\boldsymbol{x} \\|) =
e^{-\\frac{\\| \\boldsymbol{x} \\|^2}{r^2}}
\\varphi(\\boldsymbol{x}) = e^{-\\frac{\\boldsymbol{x}^2}{r^2}}

:param numpy.ndarray X: the vector x in the formula above.
:param float r: the parameter r in the formula above.

:return: result: the result of the formula above.
:rtype: float
"""
norm = np.linalg.norm(X)
result = np.exp(-(norm * norm) / (r * r))
result = np.exp(-(X * X) / (r * r))
return result

@staticmethod
Expand All @@ -161,17 +159,15 @@ def multi_quadratic_biharmonic_spline(X, r):
It implements the following formula:

.. math::
\\varphi(\\| \\boldsymbol{x} \\|) =
\\sqrt{\\| \\boldsymbol{x} \\|^2 + r^2}
\\varphi(\\boldsymbol{x}) = \\sqrt{\\boldsymbol{x}^2 + r^2}

:param numpy.ndarray X: the vector x in the formula above.
:param float r: the parameter r in the formula above.

:return: result: the result of the formula above.
:rtype: float
"""
norm = np.linalg.norm(X)
result = np.sqrt((norm * norm) + (r * r))
result = np.sqrt((X * X) + (r * r))
return result

@staticmethod
Expand All @@ -180,17 +176,16 @@ def inv_multi_quadratic_biharmonic_spline(X, r):
It implements the following formula:

.. math::
\\varphi(\\| \\boldsymbol{x} \\|) =
(\\| \\boldsymbol{x} \\|^2 + r^2 )^{-\\frac{1}{2}}
\\varphi(\\boldsymbol{x}) =
(\\boldsymbol{x}^2 + r^2 )^{-\\frac{1}{2}}

:param numpy.ndarray X: the vector x in the formula above.
:param float r: the parameter r in the formula above.

:return: result: the result of the formula above.
:rtype: float
"""
norm = np.linalg.norm(X)
result = 1.0 / (np.sqrt((norm * norm) + (r * r)))
result = 1.0 / (np.sqrt((X * X) + (r * r)))
return result

@staticmethod
Expand All @@ -199,9 +194,9 @@ def thin_plate_spline(X, r):
It implements the following formula:

.. math::
\\varphi(\\| \\boldsymbol{x} \\|) =
\\left\\| \\frac{\\boldsymbol{x} }{r} \\right\\|^2
\\ln \\left\\| \\frac{\\boldsymbol{x} }{r} \\right\\|
\\varphi(\\boldsymbol{x}) =
\\left(\\frac{\\boldsymbol{x}}{r}\\right)^2
\\ln\\frac{\\boldsymbol{x}}{r}

:param numpy.ndarray X: the vector x in the formula above.
:param float r: the parameter r in the formula above.
Expand All @@ -210,10 +205,8 @@ def thin_plate_spline(X, r):
:rtype: float
"""
arg = X / r
norm = np.linalg.norm(arg)
result = norm * norm
if norm > 0:
result *= np.log(norm)
result = arg * arg
result = np.where(arg > 0, result * np.log(arg), result)
return result

@staticmethod
Expand All @@ -222,21 +215,18 @@ def beckert_wendland_c2_basis(X, r):
It implements the following formula:

.. math::
\\varphi(\\| \\boldsymbol{x} \\|) =
\\left( 1 - \\frac{\\| \\boldsymbol{x} \\|}{r} \\right)^4_+
\\left( 4 \\frac{\\| \\boldsymbol{x} \\|}{r} + 1 \\right)
\\varphi(\\boldsymbol{x}) =
\\left( 1 - \\frac{\\boldsymbol{x}}{r}\\right)^4 +
\\left( 4 \\frac{ \\boldsymbol{x} }{r} + 1 \\right)

:param numpy.ndarray X: the vector x in the formula above.
:param float r: the parameter r in the formula above.

:return: result: the result of the formula above.
:rtype: float
"""
norm = np.linalg.norm(X)
arg = norm / r
first = 0
if (1 - arg) > 0:
first = np.power((1 - arg), 4)
arg = X / r
first = np.where((1 - arg) > 0, np.power((1 - arg), 4), 0)
second = (4 * arg) + 1
result = first * second
return result
Expand All @@ -247,18 +237,18 @@ def polyharmonic_spline(self, X, r):

.. math::

\\varphi(\\| \\boldsymbol{x} \\|) =
\\varphi(\\boldsymbol{x}) =
\\begin{cases}
\\|\\frac{\\boldsymbol{x}}{r}\\|^k
\\frac{\\boldsymbol{x}}{r}^k
\\quad & \\text{if}~k = 1,3,5,...\\\\
\\|\\frac{\\boldsymbol{x}}{r}\\|^{k-1}
\\ln(\\|\\frac{\\boldsymbol{x}}{r}\\|^
{\\|\\frac{\\boldsymbol{x}}{r}\\|})
\\quad & \\text{if}~\\|\\frac{\\boldsymbol{x}}{r}\\| < 1,
\\frac{\\boldsymbol{x}}{r}^{k-1}
\\ln(\\frac{\\boldsymbol{x}}{r}^
{\\frac{\\boldsymbol{x}}{r}})
\\quad & \\text{if}~\\frac{\\boldsymbol{x}}{r} < 1,
~k = 2,4,6,...\\\\
\\|\\frac{\\boldsymbol{x}}{r}\\|^k
\\ln(\\|\\frac{\\boldsymbol{x}}{r}\\|)
\\quad & \\text{if}~\\|\\frac{\\boldsymbol{x}}{r}\\| \\ge 1,
\\frac{\\boldsymbol{x}}{r}^k
\\ln(\\frac{\\boldsymbol{x}}{r})
\\quad & \\text{if}~\\frac{\\boldsymbol{x}}{r} \\ge 1,
~k = 2,4,6,...\\\\
\\end{cases}

Expand All @@ -270,33 +260,18 @@ def polyharmonic_spline(self, X, r):
"""

k = self.parameters.power
r_sc = np.linalg.norm(X) / r
r_sc = X / r

# k odd
if k & 1:
return np.power(r_sc, k)

print(r_sc)
# k even
if r_sc < 1:
return np.power(r_sc, k - 1) * np.log(np.power(r_sc, r_sc))
else:
return np.power(r_sc, k) * np.log(r_sc)

def _distance_matrix(self, X1, X2):
"""
This private method returns the following matrix:
:math:`\\boldsymbol{D_{ij}} = \\varphi(\\| \\boldsymbol{x_i} -
\\boldsymbol{y_j} \\|)`

:param numpy.ndarray X1: the vector x in the formula above.
:param numpy.ndarray X2: the vector y in the formula above.

:return: matrix: the matrix D.
:rtype: numpy.ndarray
"""
matrix = cdist(X1, X2,
lambda x, y: self.basis(x - y, self.parameters.radius))
return matrix
result = np.where(r_sc < 1,
np.power(r_sc, k - 1) * np.log(np.power(r_sc, r_sc)),
np.power(r_sc, k) * np.log(r_sc))
return result

def _get_weights(self, X, Y):
"""
Expand All @@ -315,17 +290,17 @@ def _get_weights(self, X, Y):
:return: weights: the matrix with the weights and the polynomial terms.
:rtype: numpy.matrix
"""
n_points = X.shape[0]
dim = X.shape[1]
identity = np.ones((n_points, 1))
dist = self._distance_matrix(X, X)
H = np.bmat([[dist, identity,
X], [identity.T,
np.zeros((1, 1)),
np.zeros((1, dim))],
[X.T, np.zeros((dim, 1)),
np.zeros((dim, dim))]])
rhs = np.bmat([[Y], [np.zeros((1, dim))], [np.zeros((dim, dim))]])
n_points, dim = X.shape
H = np.zeros((n_points + 3 + 1, n_points + 3 + 1))
H[:n_points, :n_points] = self.basis(
cdist(X, X), self.parameters.radius)
H[n_points, :n_points] = 1.0
H[:n_points, n_points] = 1.0
H[:n_points, -3:] = X
H[-3:, :n_points] = X.T

rhs = np.zeros((n_points + 3 + 1, dim))
rhs[:n_points, :] = Y
weights = np.linalg.solve(H, rhs)
return weights

Expand All @@ -335,8 +310,10 @@ def perform(self):
execution it sets `self.modified_mesh_points`.
"""
n_points = self.original_mesh_points.shape[0]
dist = self._distance_matrix(self.original_mesh_points,
self.parameters.original_control_points)
dist = self.basis(
cdist(self.original_mesh_points,
self.parameters.original_control_points),
self.parameters.radius)
identity = np.ones((n_points, 1))
H = np.bmat([[dist, identity, self.original_mesh_points]])
self.modified_mesh_points = np.asarray(np.dot(H, self.weights))
18 changes: 9 additions & 9 deletions tests/test_radial.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ def test_gaussian_spline(self):
filename='tests/test_datasets/parameters_rbf_default.prm')
rbf = rad.RBF(params, self.get_cube_mesh_points())
value = rbf.gaussian_spline(
np.array([0.5, 1, 2, 0.2]).reshape(4, 1), 0.2)
np.linalg.norm(np.array([0.5, 1, 2, 0.2])), 0.2)
np.testing.assert_almost_equal(value, 0.0)

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

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

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

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

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

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

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

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