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.

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

Add coordinate wrapping #48

Merged
merged 12 commits into from
Aug 2, 2021
8 changes: 3 additions & 5 deletions membrane_curvature/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,6 @@ 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])
print(universe.dimensions)

# Raise if selection doesn't exist
if len(self.ag) == 0:
Expand All @@ -146,10 +145,8 @@ def __init__(self, universe, select='all',
logger.warn(msg)

# Apply wrapping coordinates
if self.wrap:
self.ag.wrap()
# Warning
else:
if not self.wrap:
# Warning
msg = (" `wrap == False` may result in inaccurate calculation "
"of membrane curvature. Surfaces will be derived from "
"a reduced number of atoms. \n "
Expand All @@ -171,6 +168,7 @@ def _prepare(self):
self.n_y_bins), np.nan)

def _single_frame(self):
# Apply wrapping coordinates
if self.wrap:
self.ag.wrap()
# Populate a slice with np.arrays of surface, mean, and gaussian per frame
Expand Down
7 changes: 5 additions & 2 deletions membrane_curvature/surface.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
"""

import numpy as np
import warnings


def derive_surface(atoms, n_cells_x, n_cells_y, max_width_x, max_width_y):
Expand Down Expand Up @@ -88,14 +89,16 @@ def get_z_surface(coordinates, n_x_bins=10, n_y_bins=10, x_range=(0, 100), y_ran

try:
if l < 0 or m < 0:
lilyminium marked this conversation as resolved.
Show resolved Hide resolved
print("Atom outside grid boundaries. Skipping atom.")
msg = ("Atom outside grid boundaries. Skipping atom.")
warnings.warn(msg)
continue

grid_z_coordinates[l, m] += z
grid_norm_unit[l, m] += 1

except IndexError:
print("Atom outside grid boundaries. Skipping atom.")
msg = ("Atom outside grid boundaries. Skipping atom.")
warnings.warn(msg)

z_surface = normalized_grid(grid_z_coordinates, grid_norm_unit)

Expand Down
100 changes: 46 additions & 54 deletions membrane_curvature/tests/test_membrane_curvature.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
"""


from numpy.core.fromnumeric import mean
import pytest
from membrane_curvature.surface import normalized_grid, derive_surface, get_z_surface
from membrane_curvature.curvature import mean_curvature, gaussian_curvature
Expand Down Expand Up @@ -239,20 +238,20 @@ def universe_dummy_full(self):

@pytest.fixture()
def universe_dummy_wrap(self):
# Atoms out of bounds in x
# +-----------+
# | | 7 | 8 | 9
# | | 8 | 9 | 7
# +-----------+
# 4 | 5 | 6 | |
# 6 | 4 | 5 | |
# +-----------+
# | | 1 | 2 | 3
# | | 2 | 3 | 1
# +-----------+
a = np.array([[300., 0., 110.], [100., 0., 150.], [200., 0., 150.],
[0., 100., 150.], [100., 100., 150.], [-100., 100., 150.],
[300., 200., 110.], [100., 200., 150.], [200., 200., 150.]])

u = mda.Universe(a, n_atoms=9)
u.dimensions = [300, 300, 300, 90., 90., 90.]

return u

# Equivalent to universe_dummy_wrap when wrapping is applied.
Expand All @@ -275,6 +274,18 @@ def dummy_surface(self):
[150., 150., 150.]])
return surface

@pytest.fixture()
def curvature_unwrapped_universe(self, universe_dummy_wrap):
return MembraneCurvature(universe_dummy_wrap,
n_x_bins=3,
n_y_bins=3).run()

@pytest.fixture()
def curvature_unwrapped_universe_xy(self, universe_dummy_wrap_xy):
return MembraneCurvature(universe_dummy_wrap_xy,
n_x_bins=3,
n_y_bins=3).run()

def test_invalid_selection(self, universe):
with pytest.raises(ValueError, match=r'Invalid selection'):
MembraneCurvature(universe, select='name P')
Expand Down Expand Up @@ -366,20 +377,16 @@ def test_analysis_mean_no_wrap(self, universe):
# test using dummy Universe with atoms out of boounds
# with wrap=True (default)
# +-----------+ +-----------+
# | | 7 | 8 | 9 | 9 | 7 | 8 |
# | | 8 | 9 | 7 | 7 | 8 | 9 |
# +-----------+ +-----------+
# 4 | 5 | 6 | | ---> | 5 | 6 | 4 |
# 6 | 4 | 5 | | ---> | 4 | 5 | 6 |
# +-----------+ +-----------+
# | | 1 | 2 | 3 | 3 | 1 | 2 |
# | | 2 | 3 | 1 | 1 | 2 | 3 |
# +-----------+ +-----------+
#
# test surface in universe with atoms out of bounds in x
def test_analysis_get_z_surface_wrap(self, universe_dummy_wrap, dummy_surface):
x_bin = y_bin = 3
mc = MembraneCurvature(universe_dummy_wrap,
n_x_bins=x_bin,
n_y_bins=y_bin).run()
avg_surface = mc.results.average_z_surface
def test_analysis_get_z_surface_wrap(self, curvature_unwrapped_universe, dummy_surface):
avg_surface = curvature_unwrapped_universe.results.average_z_surface
assert_almost_equal(avg_surface, dummy_surface)

# test surface in universe with atoms out of bounds in x and y
Copy link
Member

Choose a reason for hiding this comment

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

You could use the curvature_unwrapped_universe_xy fixture here too, right? (the line after)

Copy link
Member Author

Choose a reason for hiding this comment

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

Oh, yeah, I missed that, thanks!

Expand All @@ -392,12 +399,8 @@ def test_analysis_get_z_surface_wrap_xy(self, universe_dummy_wrap_xy, dummy_surf
assert_almost_equal(avg_surface, dummy_surface)

# test mean curvature
def test_analysis_mean_wrap(self, universe_dummy_wrap, dummy_surface):
x_bin = y_bin = 3
mc = MembraneCurvature(universe_dummy_wrap,
n_x_bins=x_bin,
n_y_bins=y_bin).run()
avg_mean = mc.results.average_mean
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)
assert_almost_equal(avg_mean, expected_mean)

Expand All @@ -411,12 +414,8 @@ def test_analysis_mean_wrap_xy(self, universe_dummy_wrap_xy, dummy_surface):
assert_almost_equal(avg_mean, expected_mean)

# test gaussian
def test_analysis_gaussian_wrap(self, universe_dummy_wrap, dummy_surface):
x_bin = y_bin = 3
mc = MembraneCurvature(universe_dummy_wrap,
n_x_bins=x_bin,
n_y_bins=y_bin).run()
avg_gaussian = mc.results.average_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)
assert_almost_equal(avg_gaussian, expected_gaussian)

Expand All @@ -432,11 +431,11 @@ def test_analysis_mean_gaussian_wrap_xy(self, universe_dummy_wrap_xy, dummy_surf
# test using dummy Universe with atoms out of boounds
# with wrap=False
# +-----------+
# | | 7 | 8 | 9
# | | 8 | 9 | 7
# +-----------+
# 4 | 5 | 6 | |
# 6 | 4 | 5 | |
# +-----------+
# | | 1 | 2 | 3
# | | 2 | 3 | 1
# +-----------+
# test surface
# with wrap=False in universe with atoms out of bounds in x
Expand Down Expand Up @@ -542,35 +541,28 @@ def test_analysis_mean(self, universe_dummy_full):
avg_mean = mc.results.average_mean
assert_almost_equal(avg_mean, expected_mean)

@pytest.mark.parametrize('dummy_array, expected_surface', [
# test with negative x coordinates
(np.array([[0., 0., 150.], [-100., 0., 150.], [- 200., 0., 150.],
[0., 100., 150.], [-100., 100., 120.], [-200., 100., 120.],
[0., 200., 120.], [-100., 200., 120.], [-200., 200., 120.]]),
np.array([[150., 150., 120.],
[150., 120., 120.],
[150., 120., 120.]])),
# test with negative y coordinates
(np.array([[0., 0., 150.], [100., 0., 150.], [200., 0., 150.],
[0., -100., 150.], [100., -100., 120.], [200., -100., 120.],
[0., -200., 120.], [100., -200., 120.], [200., -200., 120.]]),
np.array([[150., 120., 150.],
[150., 120., 120.],
[150., 120., 120.]])),
# test with negative z coordinates
(np.array([[0., 0., -150.], [100., 0., -150.], [200., 0., -150.],
[0., 100., -150.], [100., 100., 120.], [200., 100., 120.],
[0., 200., 120.], [100., 200., 120.], [200., 200., 120.]]),
@pytest.mark.parametrize('x_bin, y_bin, box_dim, dummy_array, expected_surface', [
# test with negative z coordinates with 3 bins
(3, 3, 300, np.array([[0., 0., -150.], [100., 0., -150.], [200., 0., -150.],
[0., 100., -150.], [100., 100., 120.], [200., 100., 120.],
[0., 200., 120.], [100., 200., 120.], [200., 200., 120.]]),
np.array([[150., 150., 120.],
[150., 120., 120.],
[150., 120., 120.]]))
[150., 120., 120.]])),
# test with negative z coordinates with 4 bins
(4, 4, 400, np.array([[0., 0., -150.], [100., 0., -150.], [200., 0., -150.], [300., 0., -150.],
[0., 100., -150.], [100., 100., 120.], [200., 100., 120.], [300., 100., -150.],
[0., 200., 120.], [100., 200., 120.], [200., 200., 120.], [300., 200., -150.],
[0., 300., -150.], [100., 300., -150.], [200., 300., -150.], [300., 300., -150.]]),
np.array([[150., 150., 120., 150.],
[150., 120., 120., 150.],
[150., 120., 120., 150.],
[150., 150., 150., 150.]]))
])
def test_analysis_wrapping_coordinates(self, dummy_array, expected_surface):
x_bin, y_bin, = 3, 3
x_range, y_range = (0, 300), (0, 300)
def test_analysis_wrapping_coordinates(self, x_bin, y_bin, box_dim, dummy_array, expected_surface):
x_range, y_range = (0, box_dim), (0, box_dim)
u = mda.Universe(dummy_array, n_atoms=len(dummy_array))
box_dim = 300
u.dimensions = [box_dim, box_dim, box_dim, 90., 90., 90.]
u.dimensions = [box_dim, box_dim, 300, 90., 90., 90.]
# Check with wrapped coords in base
mc = MembraneCurvature(u, select='all',
n_x_bins=x_bin,
Expand Down