From 8bc0f6377099ecac6064daea3d1124b266328e5f Mon Sep 17 00:00:00 2001 From: Nicola Demo Date: Fri, 20 Apr 2018 13:50:37 +0200 Subject: [PATCH] Improve RBF performance --- .../pygem.radial.RBF._distance_matrix.rst | 6 - docs/source/radial.rst | 1 - pygem/radial.py | 129 +++++++----------- tests/test_radial.py | 18 +-- 4 files changed, 62 insertions(+), 92 deletions(-) delete mode 100644 docs/source/_summaries/pygem.radial.RBF._distance_matrix.rst diff --git a/docs/source/_summaries/pygem.radial.RBF._distance_matrix.rst b/docs/source/_summaries/pygem.radial.RBF._distance_matrix.rst deleted file mode 100644 index 40b26f9..0000000 --- a/docs/source/_summaries/pygem.radial.RBF._distance_matrix.rst +++ /dev/null @@ -1,6 +0,0 @@ -pygem.radial.RBF._distance_matrix -================================= - -.. currentmodule:: pygem.radial - -.. automethod:: RBF._distance_matrix \ No newline at end of file diff --git a/docs/source/radial.rst b/docs/source/radial.rst index ef04083..09eb76d 100644 --- a/docs/source/radial.rst +++ b/docs/source/radial.rst @@ -9,7 +9,6 @@ Radial :toctree: _summaries :nosignatures: - RBF._distance_matrix RBF._get_weights RBF.beckert_wendland_c2_basis RBF.gaussian_spline diff --git a/pygem/radial.py b/pygem/radial.py index 51718b1..1c8fb50 100644 --- a/pygem/radial.py +++ b/pygem/radial.py @@ -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: @@ -142,8 +142,7 @@ 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. @@ -151,8 +150,7 @@ def gaussian_spline(X, r): :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 @@ -161,8 +159,7 @@ 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. @@ -170,8 +167,7 @@ def multi_quadratic_biharmonic_spline(X, r): :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 @@ -180,8 +176,8 @@ 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. @@ -189,8 +185,7 @@ def inv_multi_quadratic_biharmonic_spline(X, r): :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 @@ -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. @@ -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 @@ -222,9 +215,9 @@ 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. @@ -232,11 +225,8 @@ def beckert_wendland_c2_basis(X, r): :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 @@ -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} @@ -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): """ @@ -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 @@ -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)) diff --git a/tests/test_radial.py b/tests/test_radial.py index 5b8ffcb..88f9366 100644 --- a/tests/test_radial.py +++ b/tests/test_radial.py @@ -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): @@ -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): @@ -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): @@ -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): @@ -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): @@ -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): @@ -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): @@ -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): @@ -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)