From fce9ad5f1a3a68c9d5eafb53a02369ac0eae4fdd Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Mon, 18 Dec 2023 15:10:44 -0700 Subject: [PATCH 1/3] add func to find best fit backbone vetctor --- cmeutils/geometry.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/cmeutils/geometry.py b/cmeutils/geometry.py index 089c561..6d36b1a 100644 --- a/cmeutils/geometry.py +++ b/cmeutils/geometry.py @@ -2,6 +2,33 @@ from numpy.linalg import svd +def get_backbone_vector(coordinates): + """Calculates the line of bets fit through a set of 3D coordinates. + + Best fit is calculated using numpy's Singular Value Decomposition. + + Parameters + ---------- + coordinates : numpy.ndarray, shape (N,3) + Coordinates (x,y,z) through which to fit a plane. Must be at least 3 + points. + + Returns + ------- + direction_vector : numpy.ndarray, shape (3,) + The vector of the best fit line. + """ + if coordinates.shape[0] < 2: + raise ValueError("Coordinates must contain at least 2 points.") + # Center the coordinates (subtract the mean) + centered_coordinates = coordinates - np.mean(coordinates, axis=0) + # Use PCA to find the principal components + _, _, V = np.linalg.svd(centered_coordinates) + # The first principal component (V[0]) is the vec of the best-fit line + direction_vector = V[0] + return direction_vector + + def get_plane_normal(points): """Calculate the plane which best fits a cloud of points. From bf87674ebfd2615740c5bcd7b1bd41e5f3f2b710 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 19 Dec 2023 11:54:40 -0700 Subject: [PATCH 2/3] return abs value of vector, add unit tests --- cmeutils/geometry.py | 2 +- cmeutils/tests/test_geometry.py | 19 +++++++++++++++++++ 2 files changed, 20 insertions(+), 1 deletion(-) diff --git a/cmeutils/geometry.py b/cmeutils/geometry.py index 6d36b1a..a1a2111 100644 --- a/cmeutils/geometry.py +++ b/cmeutils/geometry.py @@ -26,7 +26,7 @@ def get_backbone_vector(coordinates): _, _, V = np.linalg.svd(centered_coordinates) # The first principal component (V[0]) is the vec of the best-fit line direction_vector = V[0] - return direction_vector + return np.abs(direction_vector) def get_plane_normal(points): diff --git a/cmeutils/tests/test_geometry.py b/cmeutils/tests/test_geometry.py index 8426764..b5c59e3 100644 --- a/cmeutils/tests/test_geometry.py +++ b/cmeutils/tests/test_geometry.py @@ -3,9 +3,11 @@ import numpy as np import pytest from base_test import BaseTest +from mbuild.lib.recipes import Alkane from cmeutils.geometry import ( angle_between_vectors, + get_backbone_vector, get_plane_normal, moit, radial_grid_positions, @@ -14,6 +16,23 @@ class TestGeometry(BaseTest): + def test_backbone_vector(self): + with pytest.raises(ValueError): + coordinates = np.array([1, 1, 1]) + get_backbone_vector(coordinates) + + z_coords = np.array([[0, 0, 1], [0, 0, 2], [0, 0, 3]]) + backbone = get_backbone_vector(z_coords) + assert np.allclose(backbone, np.array([0, 0, 1]), atol=1e-5) + + x_coords = np.array([[1, 0, 0], [2, 0, 0], [3, 0, 0]]) + backbone = get_backbone_vector(x_coords) + assert np.allclose(backbone, np.array([1, 0, 0]), atol=1e-5) + + mb_chain = Alkane(n=20) + chain_backbone = get_backbone_vector(mb_chain.xyz) + assert np.allclose(chain_backbone, np.array([0, 1, 0]), atol=1e-2) + def test_moit(self): _moit = moit(points=[(-1, 0, 0), (1, 0, 0)], masses=[1, 1]) assert np.array_equal(_moit, np.array([0, 2.0, 2.0])) From 9cff1f179dd0f3c6fb10de60fe5b4990e61f3af4 Mon Sep 17 00:00:00 2001 From: chrisjonesBSU Date: Tue, 19 Dec 2023 12:12:27 -0700 Subject: [PATCH 3/3] move error check to new test --- cmeutils/tests/test_geometry.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/cmeutils/tests/test_geometry.py b/cmeutils/tests/test_geometry.py index b5c59e3..c38a577 100644 --- a/cmeutils/tests/test_geometry.py +++ b/cmeutils/tests/test_geometry.py @@ -17,10 +17,6 @@ class TestGeometry(BaseTest): def test_backbone_vector(self): - with pytest.raises(ValueError): - coordinates = np.array([1, 1, 1]) - get_backbone_vector(coordinates) - z_coords = np.array([[0, 0, 1], [0, 0, 2], [0, 0, 3]]) backbone = get_backbone_vector(z_coords) assert np.allclose(backbone, np.array([0, 0, 1]), atol=1e-5) @@ -33,6 +29,11 @@ def test_backbone_vector(self): chain_backbone = get_backbone_vector(mb_chain.xyz) assert np.allclose(chain_backbone, np.array([0, 1, 0]), atol=1e-2) + def test_backbone_vector_bad_input(self): + with pytest.raises(ValueError): + coordinates = np.array([1, 1, 1]) + get_backbone_vector(coordinates) + def test_moit(self): _moit = moit(points=[(-1, 0, 0), (1, 0, 0)], masses=[1, 1]) assert np.array_equal(_moit, np.array([0, 2.0, 2.0]))