From 86c731b20a8909446167da203cf713c3c513e02d Mon Sep 17 00:00:00 2001 From: "Christopher T. Lee" Date: Mon, 20 Nov 2023 15:49:20 -0800 Subject: [PATCH 1/3] Compute gradients with respect to grid spacing Signed-off-by: Christopher T. Lee --- membrane_curvature/base.py | 7 +- membrane_curvature/curvature.py | 28 +++--- .../tests/test_membrane_curvature.py | 96 +++++++++++++++++-- 3 files changed, 109 insertions(+), 22 deletions(-) diff --git a/membrane_curvature/base.py b/membrane_curvature/base.py index 8029dc9..0458b87 100644 --- a/membrane_curvature/base.py +++ b/membrane_curvature/base.py @@ -135,6 +135,9 @@ def __init__(self, universe, select='all', self.n_y_bins = n_y_bins self.x_range = x_range if x_range else (0, universe.dimensions[0]) self.y_range = y_range if y_range else (0, universe.dimensions[1]) + self.dx = (self.x_range[1] - self.x_range[0]) / n_x_bins + self.dy = (self.y_range[1] - self.y_range[0]) / n_y_bins + # Raise if selection doesn't exist if len(self.ag) == 0: @@ -183,8 +186,8 @@ def _single_frame(self): n_y_bins=self.n_y_bins, x_range=self.x_range, y_range=self.y_range) - self.results.mean[self._frame_index] = mean_curvature(self.results.z_surface[self._frame_index]) - self.results.gaussian[self._frame_index] = gaussian_curvature(self.results.z_surface[self._frame_index]) + self.results.mean[self._frame_index] = mean_curvature(self.results.z_surface[self._frame_index], self.dy, self.dx) + self.results.gaussian[self._frame_index] = gaussian_curvature(self.results.z_surface[self._frame_index], self.dy, self.dx) def _conclude(self): self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0) diff --git a/membrane_curvature/curvature.py b/membrane_curvature/curvature.py index 6f6f1ad..d535cfa 100644 --- a/membrane_curvature/curvature.py +++ b/membrane_curvature/curvature.py @@ -46,7 +46,7 @@ import numpy as np -def gaussian_curvature(Z): +def gaussian_curvature(Z, *varargs): """ Calculate Gaussian curvature from Z cloud points. @@ -55,7 +55,9 @@ def gaussian_curvature(Z): ---------- Z: np.ndarray. Multidimensional array of shape (n,n). - + varargs : list of scalar or array, optional + Spacing between f values. Default unitary spacing for all dimensions. + See np.gradient docs for more information. Returns ------- @@ -65,16 +67,16 @@ def gaussian_curvature(Z): """ - Zy, Zx = np.gradient(Z) - Zxy, Zxx = np.gradient(Zx) - Zyy, _ = np.gradient(Zy) + Zy, Zx = np.gradient(Z, *varargs) + Zxy, Zxx = np.gradient(Zx, *varargs) + Zyy, _ = np.gradient(Zy, *varargs) - K = (Zxx * Zyy - (Zxy ** 2)) / (1 + (Zx ** 2) + (Zy ** 2)) ** 2 + K = (Zxx * Zyy - (Zxy**2)) / (1 + (Zx**2) + (Zy**2)) ** 2 return K -def mean_curvature(Z): +def mean_curvature(Z, *varargs): """ Calculates mean curvature from Z cloud points. @@ -83,7 +85,9 @@ def mean_curvature(Z): ---------- Z: np.ndarray. Multidimensional array of shape (n,n). - + varargs : list of scalar or array, optional + Spacing between f values. Default unitary spacing for all dimensions. + See np.gradient docs for more information. Returns ------- @@ -93,11 +97,11 @@ def mean_curvature(Z): """ - Zy, Zx = np.gradient(Z) - Zxy, Zxx = np.gradient(Zx) - Zyy, _ = np.gradient(Zy) + Zy, Zx = np.gradient(Z, *varargs) + Zxy, Zxx = np.gradient(Zx, *varargs) + Zyy, _ = np.gradient(Zy, *varargs) H = (1 + Zx**2) * Zyy + (1 + Zy**2) * Zxx - 2 * Zx * Zy * Zxy - H = H / (2 * (1 + Zx**2 + Zy**2)**(1.5)) + H = H / (2 * (1 + Zx**2 + Zy**2) ** (1.5)) return H diff --git a/membrane_curvature/tests/test_membrane_curvature.py b/membrane_curvature/tests/test_membrane_curvature.py index 5c9f0ee..f4f65e3 100644 --- a/membrane_curvature/tests/test_membrane_curvature.py +++ b/membrane_curvature/tests/test_membrane_curvature.py @@ -121,7 +121,6 @@ } - @pytest.fixture def small_grofile(): u = mda.Universe(GRO_PO4_SMALL) @@ -129,6 +128,87 @@ def small_grofile(): return sel +def parametric_hemisphere( + radius: float, nx: int, ny: int, scale: float = 5 +) -> tuple[np.ndarray, float, float]: + """Generate the height field and analytical mean/gaussian curvature of a hemisphere given a sampling range. The range of the patch is defined by the hemisphere radius and a scale. + + Args: + radius (float): Radius of hemisphere + nx (int): Number of x gridpoints + ny (int): Number of y grid points + patch_scale (float, optional): Scale factor for range of grid. Defaults to 5. + + Returns: + Tuple[np.ndarray, float, float, float, float]: height field, mean curvature, gaussian curvature, dx, dy + """ + x, dx = np.linspace(-radius / scale, radius / scale, nx, retstep=True) + y, dy = np.linspace(-radius / scale, radius / scale, ny, retstep=True) + xx, yy = np.meshgrid(x, y, indexing="ij") + zz = np.sqrt(radius**2 - xx**2 - yy**2) + H = -1 / radius + K = 1 / radius**2 + return zz, H, K, dx, dy + + +def test_parametric_hemisphere_curvature(): + for radius in np.arange(5, 100, 5): + for nx in [7, 23, 47, 67]: + for ny in [7, 23, 47, 67]: + # Ensure that patch_scale is sufficiently large such that points are closely spaced + height_field, H, K, dx, dy = parametric_hemisphere( + radius=radius, nx=nx, ny=ny, scale=200 + ) + mean = mean_curvature(height_field, dx, dy) + # Check only central points to avoid boundary errors from np.gradient + assert np.all(np.isclose(mean[2:-2, 2:-2], H)) + gauss = gaussian_curvature(height_field, dx, dy) + assert np.all(np.isclose(gauss[2:-2, 2:-2], K)) + + +def parametric_saddle( + length: float, nx: int, ny: int, scale: float = 1000 +) -> tuple[np.ndarray, float, float]: + """Generate the height field and analytical mean/gaussian curvature corresponding to a monkey saddle. + + Args: + length (float): Sample from - to + length. + nx (int): Number of x grid points + ny (int): Number of y grid points + scale (float, optional): Scale factor for range of grid. Defaults to 500. + + Returns: + tuple[np.ndarray, np.ndarray, np.ndarray, float, float]: height field, mean curvature, gaussian curvature, dx, dy + """ + x, dx = np.linspace(-length / scale, length / scale, nx, retstep=True) + y, dy = np.linspace(-length / scale, length / scale, ny, retstep=True) + xx, yy = np.meshgrid(x, y, indexing="ij") + zz = xx * (xx**2 - 3 * yy**2) # Monkey saddle + + K = -36 * (xx**2 + yy**2) / (1 + 9 * (xx**2 + yy**2) ** 2) ** 2 + H = ( + 27 + * xx + * (-(xx**4) + 2 * xx**2 * yy**2 + 3 * yy**4) + / (1 + 9 * (xx**2 + yy**2) ** 2) ** (3 / 2) + ) + return zz, K, H, dx, dy + + +def test_parametric_saddle_curvature(): + for nx in [7, 23, 47, 67]: + for ny in [7, 23, 47, 67]: + height_field, K, H, dx, dy = parametric_saddle( + length=10, nx=nx, ny=ny, scale=800 + ) + mean = mean_curvature(height_field, dx, dy) + gauss = gaussian_curvature(height_field, dx, dy) + + # Check only central points to avoid boundary errors from np.gradient + assert np.all(np.isclose(mean[2:-2, 2:-2], H[2:-2, 2:-2])) + assert np.all(np.isclose(gauss[2:-2, 2:-2], K[2:-2, 2:-2])) + + def test_gaussian_curvature_small(): K_test = gaussian_curvature(MEMBRANE_CURVATURE_DATA['z_ref']['small']) for k, k_test in zip(MEMBRANE_CURVATURE_DATA['gaussian_curvature']['small'], K_test): @@ -401,23 +481,23 @@ def test_analysis_get_z_surface_wrap_xy(self, universe_dummy_wrap_xy, dummy_surf # test mean curvature def test_analysis_mean_wrap(self, curvature_unwrapped_universe, dummy_surface): avg_mean = curvature_unwrapped_universe.results.average_mean - expected_mean = mean_curvature(dummy_surface) + expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) assert_almost_equal(avg_mean, expected_mean) def test_analysis_mean_wrap_xy(self, curvature_unwrapped_universe, dummy_surface): avg_mean = curvature_unwrapped_universe.results.average_mean - expected_mean = mean_curvature(dummy_surface) + expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) assert_almost_equal(avg_mean, expected_mean) # test gaussian def test_analysis_gaussian_wrap(self, curvature_unwrapped_universe, dummy_surface): avg_gaussian = curvature_unwrapped_universe.results.average_gaussian - expected_gaussian = gaussian_curvature(dummy_surface) + expected_gaussian = gaussian_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) assert_almost_equal(avg_gaussian, expected_gaussian) def test_analysis_mean_gaussian_wrap_xy(self, curvature_unwrapped_universe, dummy_surface): avg_gaussian = curvature_unwrapped_universe.results.average_gaussian - expected_gaussian = gaussian_curvature(dummy_surface) + expected_gaussian = gaussian_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) assert_almost_equal(avg_gaussian, expected_gaussian) # test using dummy Universe with atoms out of boounds @@ -524,9 +604,9 @@ def test_analysis_get_z_surface(self, universe_dummy_full, x_bin, y_bin, expecte assert_almost_equal(avg_surface, expected_surface) def test_analysis_mean(self, universe_dummy_full): - expected_mean = np.array([[5.54630914e-04, 1.50000000e+01, -8.80203593e-02], - [2.77315457e-04, 2.20748929e-03, 5.01100068e-01], - [2.77315457e-04, 2.20748929e-03, 1.50000000e+01]]) + expected_mean = np.array([[0.001318109566818, 0.0015 , 0.001702672969121], + [0.000659054783409, 0.000725381187923, 0.001467083452574], + [0.000659054783409, 0.000725381187923, 0.0015 ]]) mc = MembraneCurvature(universe_dummy_full, n_x_bins=3, n_y_bins=3).run() From 12879196b85f29a33c11ee5b11a72979f1afde2f Mon Sep 17 00:00:00 2001 From: "Christopher T. Lee" Date: Tue, 19 Dec 2023 16:28:56 -0800 Subject: [PATCH 2/3] Fixing doc style and flake8 issues Signed-off-by: Christopher T. Lee --- membrane_curvature/base.py | 9 ++- .../tests/test_membrane_curvature.py | 72 ++++++++++++------- 2 files changed, 53 insertions(+), 28 deletions(-) diff --git a/membrane_curvature/base.py b/membrane_curvature/base.py index 0458b87..10fdec3 100644 --- a/membrane_curvature/base.py +++ b/membrane_curvature/base.py @@ -138,7 +138,6 @@ def __init__(self, universe, select='all', self.dx = (self.x_range[1] - self.x_range[0]) / n_x_bins self.dy = (self.y_range[1] - self.y_range[0]) / n_y_bins - # Raise if selection doesn't exist if len(self.ag) == 0: raise ValueError("Invalid selection. Empty AtomGroup.") @@ -186,8 +185,12 @@ def _single_frame(self): n_y_bins=self.n_y_bins, x_range=self.x_range, y_range=self.y_range) - self.results.mean[self._frame_index] = mean_curvature(self.results.z_surface[self._frame_index], self.dy, self.dx) - self.results.gaussian[self._frame_index] = gaussian_curvature(self.results.z_surface[self._frame_index], self.dy, self.dx) + self.results.mean[self._frame_index] = mean_curvature( + self.results.z_surface[self._frame_index], self.dy, self.dx + ) + self.results.gaussian[self._frame_index] = gaussian_curvature( + self.results.z_surface[self._frame_index], self.dy, self.dx + ) def _conclude(self): self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0) diff --git a/membrane_curvature/tests/test_membrane_curvature.py b/membrane_curvature/tests/test_membrane_curvature.py index f4f65e3..d56c722 100644 --- a/membrane_curvature/tests/test_membrane_curvature.py +++ b/membrane_curvature/tests/test_membrane_curvature.py @@ -121,6 +121,7 @@ } + @pytest.fixture def small_grofile(): u = mda.Universe(GRO_PO4_SMALL) @@ -131,16 +132,25 @@ def small_grofile(): def parametric_hemisphere( radius: float, nx: int, ny: int, scale: float = 5 ) -> tuple[np.ndarray, float, float]: - """Generate the height field and analytical mean/gaussian curvature of a hemisphere given a sampling range. The range of the patch is defined by the hemisphere radius and a scale. - - Args: - radius (float): Radius of hemisphere - nx (int): Number of x gridpoints - ny (int): Number of y grid points - patch_scale (float, optional): Scale factor for range of grid. Defaults to 5. - - Returns: - Tuple[np.ndarray, float, float, float, float]: height field, mean curvature, gaussian curvature, dx, dy + """Generate the height field and analytical mean/gaussian curvature of a + hemisphere given a sampling range. The range of the patch is defined by + the hemisphere radius and a scale. + + Parameters + ---------- + radius: float + Radius of hemisphere + nx: int + Number of x gridpoints + ny: int + Number of y grid points + patch_scale: float, optional + Scale factor for range of grid. Defaults to 5. + + Returns + ------- + Tuple[np.ndarray, float, float, float, float] + height field, mean curvature, gaussian curvature, dx, dy """ x, dx = np.linspace(-radius / scale, radius / scale, nx, retstep=True) y, dy = np.linspace(-radius / scale, radius / scale, ny, retstep=True) @@ -169,16 +179,24 @@ def test_parametric_hemisphere_curvature(): def parametric_saddle( length: float, nx: int, ny: int, scale: float = 1000 ) -> tuple[np.ndarray, float, float]: - """Generate the height field and analytical mean/gaussian curvature corresponding to a monkey saddle. - - Args: - length (float): Sample from - to + length. - nx (int): Number of x grid points - ny (int): Number of y grid points - scale (float, optional): Scale factor for range of grid. Defaults to 500. - - Returns: - tuple[np.ndarray, np.ndarray, np.ndarray, float, float]: height field, mean curvature, gaussian curvature, dx, dy + """Generate the height field and analytical mean/gaussian curvature + corresponding to a monkey saddle. + + Parameters + ---------- + length: float + Sample from - to + length. + nx: int + Number of x grid points + ny: int + Number of y grid points + scale: float, optional + Scale factor for range of grid. Defaults to 500. + + Returns + ------- + tuple[np.ndarray, np.ndarray, np.ndarray, float, float] + height field, mean curvature, gaussian curvature, dx, dy """ x, dx = np.linspace(-length / scale, length / scale, nx, retstep=True) y, dy = np.linspace(-length / scale, length / scale, ny, retstep=True) @@ -492,12 +510,16 @@ def test_analysis_mean_wrap_xy(self, curvature_unwrapped_universe, dummy_surface # test gaussian def test_analysis_gaussian_wrap(self, curvature_unwrapped_universe, dummy_surface): avg_gaussian = curvature_unwrapped_universe.results.average_gaussian - expected_gaussian = gaussian_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) + expected_gaussian = gaussian_curvature(dummy_surface, + curvature_unwrapped_universe.dy, + curvature_unwrapped_universe.dx) assert_almost_equal(avg_gaussian, expected_gaussian) def test_analysis_mean_gaussian_wrap_xy(self, curvature_unwrapped_universe, dummy_surface): avg_gaussian = curvature_unwrapped_universe.results.average_gaussian - expected_gaussian = gaussian_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) + expected_gaussian = gaussian_curvature(dummy_surface, + curvature_unwrapped_universe.dy, + curvature_unwrapped_universe.dx) assert_almost_equal(avg_gaussian, expected_gaussian) # test using dummy Universe with atoms out of boounds @@ -604,9 +626,9 @@ def test_analysis_get_z_surface(self, universe_dummy_full, x_bin, y_bin, expecte assert_almost_equal(avg_surface, expected_surface) def test_analysis_mean(self, universe_dummy_full): - expected_mean = np.array([[0.001318109566818, 0.0015 , 0.001702672969121], - [0.000659054783409, 0.000725381187923, 0.001467083452574], - [0.000659054783409, 0.000725381187923, 0.0015 ]]) + expected_mean = np.array([[0.001318109566818, 0.0015, 0.001702672969121], + [0.000659054783409, 0.000725381187923, 0.001467083452574], + [0.000659054783409, 0.000725381187923, 0.0015]]) mc = MembraneCurvature(universe_dummy_full, n_x_bins=3, n_y_bins=3).run() From 965530ac7459c373c702c269e42b8f1ce639f738 Mon Sep 17 00:00:00 2001 From: "Christopher T. Lee" Date: Tue, 20 Feb 2024 20:42:26 -0800 Subject: [PATCH 3/3] fix: Fix gradient dx, dy order Signed-off-by: Christopher T. Lee --- membrane_curvature/base.py | 4 ++-- membrane_curvature/curvature.py | 12 ++++++------ membrane_curvature/tests/test_membrane_curvature.py | 12 ++++++------ 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/membrane_curvature/base.py b/membrane_curvature/base.py index 10fdec3..364e9d5 100644 --- a/membrane_curvature/base.py +++ b/membrane_curvature/base.py @@ -186,10 +186,10 @@ def _single_frame(self): x_range=self.x_range, y_range=self.y_range) self.results.mean[self._frame_index] = mean_curvature( - self.results.z_surface[self._frame_index], self.dy, self.dx + self.results.z_surface[self._frame_index], self.dx, self.dy ) self.results.gaussian[self._frame_index] = gaussian_curvature( - self.results.z_surface[self._frame_index], self.dy, self.dx + self.results.z_surface[self._frame_index], self.dx, self.dy ) def _conclude(self): diff --git a/membrane_curvature/curvature.py b/membrane_curvature/curvature.py index d535cfa..fcde91b 100644 --- a/membrane_curvature/curvature.py +++ b/membrane_curvature/curvature.py @@ -67,9 +67,9 @@ def gaussian_curvature(Z, *varargs): """ - Zy, Zx = np.gradient(Z, *varargs) - Zxy, Zxx = np.gradient(Zx, *varargs) - Zyy, _ = np.gradient(Zy, *varargs) + Zx, Zy = np.gradient(Z, *varargs) + Zxx, Zxy = np.gradient(Zx, *varargs) + _, Zyy = np.gradient(Zy, *varargs) K = (Zxx * Zyy - (Zxy**2)) / (1 + (Zx**2) + (Zy**2)) ** 2 @@ -97,9 +97,9 @@ def mean_curvature(Z, *varargs): """ - Zy, Zx = np.gradient(Z, *varargs) - Zxy, Zxx = np.gradient(Zx, *varargs) - Zyy, _ = np.gradient(Zy, *varargs) + Zx, Zy, = np.gradient(Z, *varargs) + Zxx, Zxy = np.gradient(Zx, *varargs) + _, Zyy = np.gradient(Zy, *varargs) H = (1 + Zx**2) * Zyy + (1 + Zy**2) * Zxx - 2 * Zx * Zy * Zxy H = H / (2 * (1 + Zx**2 + Zy**2) ** (1.5)) diff --git a/membrane_curvature/tests/test_membrane_curvature.py b/membrane_curvature/tests/test_membrane_curvature.py index d56c722..03c6ef8 100644 --- a/membrane_curvature/tests/test_membrane_curvature.py +++ b/membrane_curvature/tests/test_membrane_curvature.py @@ -499,27 +499,27 @@ def test_analysis_get_z_surface_wrap_xy(self, universe_dummy_wrap_xy, dummy_surf # test mean curvature def test_analysis_mean_wrap(self, curvature_unwrapped_universe, dummy_surface): avg_mean = curvature_unwrapped_universe.results.average_mean - expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) + expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dx, curvature_unwrapped_universe.dy) assert_almost_equal(avg_mean, expected_mean) def test_analysis_mean_wrap_xy(self, curvature_unwrapped_universe, dummy_surface): avg_mean = curvature_unwrapped_universe.results.average_mean - expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx) + expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dx, curvature_unwrapped_universe.dy) assert_almost_equal(avg_mean, expected_mean) # test gaussian def test_analysis_gaussian_wrap(self, curvature_unwrapped_universe, dummy_surface): avg_gaussian = curvature_unwrapped_universe.results.average_gaussian expected_gaussian = gaussian_curvature(dummy_surface, - curvature_unwrapped_universe.dy, - curvature_unwrapped_universe.dx) + curvature_unwrapped_universe.dx, + curvature_unwrapped_universe.dy) assert_almost_equal(avg_gaussian, expected_gaussian) def test_analysis_mean_gaussian_wrap_xy(self, curvature_unwrapped_universe, dummy_surface): avg_gaussian = curvature_unwrapped_universe.results.average_gaussian expected_gaussian = gaussian_curvature(dummy_surface, - curvature_unwrapped_universe.dy, - curvature_unwrapped_universe.dx) + curvature_unwrapped_universe.dx, + curvature_unwrapped_universe.dy) assert_almost_equal(avg_gaussian, expected_gaussian) # test using dummy Universe with atoms out of boounds