Skip to content
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

Jump to bottom

Compute gradients with respect to grid spacing #116

Merged
merged 6 commits into from
Aug 20, 2024
Merged

Compute gradients with respect to grid spacing #116

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
86c731b
Compute gradients with respect to grid spacing
ctlee Nov 20, 2023
4aa80be
Merge remote-tracking branch 'upstream/main'
ctlee Dec 19, 2023
1287919
Fixing doc style and flake8 issues
ctlee Dec 20, 2023
965530a
fix: Fix gradient dx, dy order
ctlee Feb 21, 2024
bea5fdd
Merge branch 'MDAnalysis:main' into main
ctlee Jun 17, 2024
35ceed9
Merge branch 'MDAnalysis:main' into main
ctlee Aug 20, 2024
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
10 changes: 8 additions & 2 deletions membrane_curvature/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,8 @@ 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:
Expand Down Expand Up @@ -183,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.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)
Expand Down
28 changes: 16 additions & 12 deletions membrane_curvature/curvature.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@
import numpy as np


def gaussian_curvature(Z):
def gaussian_curvature(Z, *varargs):
ojeda-e marked this conversation as resolved.
Show resolved Hide resolved
"""
Calculate Gaussian curvature from Z cloud points.

Expand All @@ -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
-------
Expand All @@ -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.

Expand All @@ -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
-------
Expand All @@ -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
116 changes: 109 additions & 7 deletions membrane_curvature/tests/test_membrane_curvature.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,104 @@ 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.

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)
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)
ojeda-e marked this conversation as resolved.
Show resolved Hide resolved
# 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.

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)
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)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can also check that has 0 gaussian curvature at origin?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The 0 gaussian curvature at origin condition is implicitly checked by conditions where nx and ny are odd. If the discretization doesn't include (0, 0) then we'd have to interpolate or skip the test. I think comparison against the analytical formula at corresponding grid points should be robust and complete enough.

# 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):
Expand Down Expand Up @@ -401,23 +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)
expected_mean = mean_curvature(dummy_surface, curvature_unwrapped_universe.dy, curvature_unwrapped_universe.dx)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Parametrize over these arguments, to check that the original implementation still functions as intended.

@pytest.mark.parametrize("varargs", [, [], [dy, dx]])
``
Or if the parametrization is too complicated split into two tests.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure what to do here. The original implementation is fundamentally flawed.

The test fixture universe_dummy_wrap has a specific grid spacing and system dimension. Taking derivatives with respect to the default dx, dy value (=1) is like saying:

$f'(x) = \frac{f(x+h) -f(x)}{1}$.

When the derivatives must be taken with respect to the spacing, $h$:

$f'(x) = \frac{f(x+h) -f(x)}{h}$.

To have dx, dy as parameters, the fixture would have to be changed to produce meshes with varied discretization.

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)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same as before

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
Expand Down Expand Up @@ -524,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([[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()
Expand Down