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

Surface Interpolation Function [Fixes Issue45] #52

Closed
wants to merge 8 commits into from
Closed

Surface Interpolation Function [Fixes Issue45] #52

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
33b9748
Added interpolation function with docstrings
Jul 22, 2021
13ea4a9
Added surface_interpolation tests
Jul 22, 2021
36bd38e
Fixed PEP8
Jul 22, 2021
7382cec
Fixed PEP8
Jul 22, 2021
10fd0a7
inverted original mask_nans in interpolation function
Jul 25, 2021
fc4ec95
Replaced assertions by assert_allclose.
Aug 5, 2021
3a2ad5d
Added interpolation to base with tests.
Aug 5, 2021
7d7f4ca
changed atol in assert_allclose
Aug 7, 2021
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
32 changes: 24 additions & 8 deletions membrane_curvature/base.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@

import numpy as np
import warnings
from .surface import get_z_surface
from .surface import get_z_surface, surface_interpolation
from .curvature import mean_curvature, gaussian_curvature

import MDAnalysis
Expand All @@ -32,7 +32,7 @@ class MembraneCurvature(AnalysisBase):
----------
universe : Universe or AtomGroup
An MDAnalysis Universe object.
select : str or iterable of str, optional.
select : str or iterable of str, optional.
The selection string of an atom selection to use as a
reference to derive a surface.
wrap : bool, optional
Expand All @@ -57,13 +57,13 @@ class MembraneCurvature(AnalysisBase):
results.gaussian_curvature : ndarray
Gaussian curvature associated to the surface.
Arrays of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.average_z_surface : ndarray
Average of the array elements in `z_surface`.
results.average_z_surface : ndarray
Average of the array elements in `z_surface`.
Each array has shape (`n_x_bins`, `n_y_bins`)
results.average_mean_curvature : ndarray
results.average_mean_curvature : ndarray
Average of the array elements in `mean_curvature`.
Each array has shape (`n_x_bins`, `n_y_bins`)
results.average_gaussian_curvature: ndarray
results.average_gaussian_curvature: ndarray
Average of the array elements in `gaussian_curvature`.
Each array has shape (`n_x_bins`, `n_y_bins`)

Expand Down Expand Up @@ -120,7 +120,8 @@ def __init__(self, universe, select='all',
n_x_bins=100, n_y_bins=100,
x_range=None,
y_range=None,
wrap=True, **kwargs):
wrap=True,
interpolation=False, **kwargs):

super().__init__(universe.universe.trajectory, **kwargs)
self.ag = universe.select_atoms(select)
Expand All @@ -129,6 +130,7 @@ 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.interpolation = interpolation

# Raise if selection doesn't exist
if len(self.ag) == 0:
Expand Down Expand Up @@ -156,7 +158,6 @@ def __init__(self, universe, select='all',
logger.warn(msg)

def _prepare(self):
# Initialize empty np.array with results
self.results.z_surface = np.full((self.n_frames,
self.n_x_bins,
self.n_y_bins), np.nan)
Expand All @@ -166,6 +167,9 @@ def _prepare(self):
self.results.gaussian = np.full((self.n_frames,
self.n_x_bins,
self.n_y_bins), np.nan)
self.results.interpolated_z_surface = np.full((self.n_frames,
self.n_x_bins,
self.n_y_bins), np.nan)

def _single_frame(self):
# Apply wrapping coordinates
Expand All @@ -179,8 +183,20 @@ def _single_frame(self):
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])
# Perform interpolation fo surface
if self.interpolation:
# Populate a slice with np.arrays of interpolated surface
self.results.interpolated_z_surface[self._frame_index] = surface_interpolation(
self.results.z_surface[self._frame_index])

def _conclude(self):
self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0)
self.results.average_mean = np.nanmean(self.results.mean, axis=0)
self.results.average_gaussian = np.nanmean(self.results.gaussian, axis=0)
if self.interpolation:
self.results.average_interpolated_z_surface = np.nanmean(
self.results.interpolated_z_surface[self._frame_index], axis=0)
Copy link
Member

Choose a reason for hiding this comment

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

Suggested change
self.results.interpolated_z_surface[self._frame_index], axis=0)
self.results.interpolated_z_surface, axis=0)

And with the others too -- if you specify frame_index, you'll wind up only getting the mean of the last frame.

Copy link
Member

Choose a reason for hiding this comment

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

Do you need np.nanmean? Do you expect np.nan values in your interpolated surfaces? If not, using np.mean will reveal errors in case there ever are np.nan values.

self.results.average_interpolated_mean = np.nanmean(mean_curvature(
self.results.interpolated_z_surface[self._frame_index]), axis=0)
self.results.average_interpolated_gaussian = np.nanmean(gaussian_curvature(
self.results.interpolated_z_surface[self._frame_index]), axis=0)
65 changes: 61 additions & 4 deletions membrane_curvature/surface.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ def derive_surface(atoms, n_cells_x, n_cells_y, max_width_x, max_width_y):

Returns
-------
z_coordinates: numpy.ndarray
z_coordinates: np.ndarray
Average z-coordinate values. Return Numpy array of floats of
shape `(n_cells_x, n_cells_y)`.

Expand All @@ -55,8 +55,8 @@ def get_z_surface(coordinates, n_x_bins=10, n_y_bins=10, x_range=(0, 100), y_ran

Parameters
----------
coordinates : numpy.ndarray
Coordinates of AtomGroup. Numpy array of shape=(n_atoms, 3).
coordinates : np.ndarray
Coordinates of AtomGroup. NumPy array of shape=(n_atoms, 3).
n_x_bins : int.
Number of bins in grid in the `x` dimension.
n_y_bins : int.
Expand All @@ -70,7 +70,7 @@ def get_z_surface(coordinates, n_x_bins=10, n_y_bins=10, x_range=(0, 100), y_ran
-------
z_surface: np.ndarray
Surface derived from set of coordinates in grid of `x_range, y_range` dimensions.
Returns Numpy array of floats of shape (`n_x_bins`, `n_y_bins`)
Returns NumPy array of floats of shape (`n_x_bins`, `n_y_bins`)

"""

Expand Down Expand Up @@ -131,3 +131,60 @@ def normalized_grid(grid_z_coordinates, grid_norm_unit):
z_normalized = grid_z_coordinates / grid_norm_unit

return z_normalized


def interpolation_by_array(array_surface):
"""
Interpolates values contained in `array_surface` over axis
Copy link
Member

Choose a reason for hiding this comment

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

Over which axis?


Parameters
----------

array_surface: np.ndarray
Array of floats of shape (`n_x_bins`, `n_y_bins`)

Returns
-------
interpolated_array: np.ndarray
Returns interpolated array.
Array of shape (`n_x_bins` x `n_y_bins`,)

"""

# create mask for nans
mask_nans = ~np.isnan(array_surface)

# index of array_surface
index_array = np.arange(array_surface.shape[0])

# interpolate values in array
interpolated_array = np.interp(index_array,
np.flatnonzero(mask_nans),
array_surface[mask_nans])

return interpolated_array
Comment on lines +155 to +165
Copy link
Member

Choose a reason for hiding this comment

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

I'm not quite understanding this. mask_nans is 2D, with shape (n_x, n_y). index_array is 1D, with shape n_x. np.flatnonzero(mask_nans) and array_surface[mask_nans] will be 1D, of length between 0 to n_x * n_y. So I think that means that you are only ever operating on the first row of array_surface?

I guess my question is then, why is array_surface 2D? Why not just pass in one row at a time? In addition, as it is a 2D surface, why use a 1D interpolation function? Why not a 2D one?



def surface_interpolation(array_surface):
"""
Calculates interpolation of arrays.

Parameters
----------

array_surface: np.ndarray
Array of floats of shape (`n_x_bins`, `n_y_bins`)

Returns
-------

interpolated_surface: np.ndarray
Interpolated surface derived from set of coordinates
in grid of `x_range, y_range` dimensions.
Array of floats of shape (`n_x_bins`, `n_y_bins`)

"""

interpolated_surface = np.apply_along_axis(interpolation_by_array, axis=0, arr=array_surface)
ojeda-e marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Member

Choose a reason for hiding this comment

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

Ok, I see that you are applying interpolation_by_array by row, then. In that case I suggest amending the documentation of interpolation_by_array to specify that the input is 1-dimensional with length n_x, instead of n_x, n_y. Although I think doing 2D interpolation would be less arbitrary (choosing which axis is x and which is y is basically a random, right?)


return interpolated_surface
Loading