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

Compute gradients with respect to grid spacing #116

Merged
merged 6 commits into from
Aug 20, 2024
Merged

Conversation

ctlee
Copy link
Contributor

@ctlee ctlee commented Nov 20, 2023

Description

The calculation of the gradients of the height field should respect the grid spacing. This PR adds support for variable grid spacing. It also adds some tests against parametric surfaces.

fixes #115

Todos

Notable points that this PR has either accomplished or will accomplish.

  • Make np.gradient respect grid spacing
  • Add test to compare computed mean and gaussian curvature against analytical formula for hemisphere
  • Add test to compare computed mean and gaussian curvature against analytical formula for monkey saddle
  • Fix unit tests to check against curvature calculations which respect grid spacing.

Questions

N/A

Status

  • Ready to go

Signed-off-by: Christopher T. Lee <ctlee@ucsd.edu>
@pep8speaks
Copy link

pep8speaks commented Nov 20, 2023

Hello @ctlee! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

Line 138:65: W291 trailing whitespace

Line 135:77: W291 trailing whitespace
Line 136:76: W291 trailing whitespace
Line 145:16: W291 trailing whitespace
Line 168:80: E501 line too long (99 > 79 characters)
Line 173:80: E501 line too long (85 > 79 characters)
Line 189:16: W291 trailing whitespace
Line 193:31: W291 trailing whitespace
Line 225:80: E501 line too long (81 > 79 characters)
Line 502:80: E501 line too long (119 > 79 characters)
Line 507:80: E501 line too long (119 > 79 characters)
Line 513:62: W291 trailing whitespace
Line 629:80: E501 line too long (81 > 79 characters)
Line 630:80: E501 line too long (91 > 79 characters)
Line 631:80: E501 line too long (81 > 79 characters)

Comment last updated at 2024-08-20 10:36:04 UTC

Copy link

codecov bot commented Nov 21, 2023

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 100.00%. Comparing base (44d136b) to head (35ceed9).
Report is 9 commits behind head on main.

Additional details and impacted files

@IAlibay
Copy link
Member

IAlibay commented Dec 19, 2023

@ctlee Apologies for the very late response here - unfortunately membrane curvature doesn't have an active maintainer currently so there haven't been anyone to shephred pull requests.

We will try to get someone to review your code as soon as possible. Pinging @MDAnalysis/coredevs here in case someone might want to take on an initial review to help this process along.

Copy link
Member

@hmacdope hmacdope left a comment

Choose a reason for hiding this comment

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

Initial review, :)

You may need to investigate some changes to docs, but otherwise great start,

membrane_curvature/tests/test_membrane_curvature.py Outdated Show resolved Hide resolved
)
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.

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

@IAlibay
Copy link
Member

IAlibay commented Dec 19, 2023

Updating against main should help get CI working properly here (since we've now merged #112 )

Copy link
Contributor Author

@ctlee ctlee left a comment

Choose a reason for hiding this comment

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

Thanks @hmacdope for the initial review and @IAlibay for moving this along and fixing CI. I've added a commit to change docs to numpy style.

Please advise on what to do wrt the other requested changes.

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

@hmacdope
Copy link
Member

Sounds reasonble @ctlee, I am happy to give my approval on a purely code basis, but from a scientific basis I will have to take your assertions about the correctness of the new implementation relative to the old one on trust. I will tag @ojeda-e as the principal developer as she may have a view here. I will leave my view blocking for now pending waiting for feedback. :)

@hmacdope hmacdope requested a review from ojeda-e January 5, 2024 03:07
@ctlee
Copy link
Contributor Author

ctlee commented Jan 29, 2024

Following up on this as new versions appear to be released! Folks will get incorrect curvature values if their grid spacing isn't 1. Can verify this issue in the current implementation by comparing with the parametric surface tests added in this PR.

An alternative approach is to add a check that the dx, dy is 1 and to throw an error otherwise. @ojeda-e

Copy link
Member

@ojeda-e ojeda-e left a comment

Choose a reason for hiding this comment

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

Hi @ctlee
Thanks for working on this changes! I left two question below to clarify my understanding.
Would it be possible to offer an example where*varargs is illustrated? It would be a nice addition to the docs.

membrane_curvature/curvature.py Show resolved Hide resolved
Signed-off-by: Christopher T. Lee <ctlee@ucsd.edu>
@ctlee ctlee requested a review from ojeda-e February 21, 2024 04:53
@hmacdope
Copy link
Member

@ojeda-e are we happy with this PR as is? I would be happy to go ahead, but thought would check with you?

Copy link
Member

@ojeda-e ojeda-e left a comment

Choose a reason for hiding this comment

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

Hi @ctlee @hmacdope,
thanks for following up and pinging me here. Yes, I think the changes are fine. I am approving while codecov is fixed (or removed) in the CI.

@ctlee ctlee requested a review from hmacdope August 20, 2024 10:37
@hmacdope hmacdope merged commit 2019357 into MDAnalysis:main Aug 20, 2024
16 of 18 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Gradient calculation should respect arbitrary grid spacing
5 participants