From eda1512edbec95c34b90b605cdbd75b31bdae714 Mon Sep 17 00:00:00 2001 From: Tyler Heim Date: Mon, 28 Mar 2022 14:04:12 -0700 Subject: [PATCH 01/14] Random rotations --- .../expectation_maximization.py | 30 +++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index 47c61ee..dce344c 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -1,5 +1,6 @@ """Iterative refinement with Bayesian expectation maximization.""" +import healpy as hp import numpy as np from simSPI.transfer import eval_ctf @@ -274,8 +275,33 @@ def grid_SO3_uniform(n_rotations): Array describing rotations. Shape (n_rotations, 3, 3) """ - rots = np.ones((n_rotations, 3, 3)) - return rots + try: + # If n_rotations fits the healpix rotations, + # we can do a proper uniform map. + nside = hp.npix2nside(n_rotations) + + rots = np.ones((n_rotations, 3, 3)) + return rots + except ValueError: + # If n_rotations can't be evenly distributed + # on a sphere by healpix, use random rots. + phi = np.random.random(n_rotations) * 2 * np.pi() + theta = np.arccos(np.random.random(n_rotations)) * (np.random.randint(2) * 2 - 1) + phi_matrices = np.array( + [ + [np.cos(phi), -np.sin(phi), 0], + [np.sin(phi), np.cos(phi), 0], + [0, 0, 1] + ] + ) + theta_matrices = np.array( + [ + [1, 0, 0], + [0, np.cos(theta), -np.sin(theta)], + [0, np.sin(theta), np.cos(theta)] + ] + ) + return np.dot(phi_matrices, theta_matrices) @staticmethod def generate_xy_plane(n_pix): From ca187c88aa3381e0c024dfd7ab23656c482807ef Mon Sep 17 00:00:00 2001 From: Tyler Heim Date: Mon, 28 Mar 2022 14:09:01 -0700 Subject: [PATCH 02/14] np.dot was the wrong operator. This is too, but its closer to the correct one --- reconstructSPI/iterative_refinement/expectation_maximization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index dce344c..3576fce 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -301,7 +301,7 @@ def grid_SO3_uniform(n_rotations): [0, np.sin(theta), np.cos(theta)] ] ) - return np.dot(phi_matrices, theta_matrices) + return phi_matrices * theta_matrices @staticmethod def generate_xy_plane(n_pix): From f4b7d136027196118959680454569f92a3244603 Mon Sep 17 00:00:00 2001 From: thisTyler Date: Mon, 28 Mar 2022 23:46:54 -0700 Subject: [PATCH 03/14] Implemented random version (correctly) --- .../expectation_maximization.py | 49 +++++++++---------- 1 file changed, 24 insertions(+), 25 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index 3576fce..f3da42c 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -275,33 +275,32 @@ def grid_SO3_uniform(n_rotations): Array describing rotations. Shape (n_rotations, 3, 3) """ - try: - # If n_rotations fits the healpix rotations, - # we can do a proper uniform map. - nside = hp.npix2nside(n_rotations) - - rots = np.ones((n_rotations, 3, 3)) - return rots - except ValueError: - # If n_rotations can't be evenly distributed - # on a sphere by healpix, use random rots. - phi = np.random.random(n_rotations) * 2 * np.pi() - theta = np.arccos(np.random.random(n_rotations)) * (np.random.randint(2) * 2 - 1) - phi_matrices = np.array( - [ - [np.cos(phi), -np.sin(phi), 0], - [np.sin(phi), np.cos(phi), 0], - [0, 0, 1] - ] + rots = np.zeros((n_rotations, 3, 3)) + for i in range(n_rotations): + phi, theta, psi = np.random.random(3) * 2 * np.pi + phi_matrix = np.array( + ( + (np.cos(phi), -np.sin(phi), 0), + (np.sin(phi), np.cos(phi), 0), + (0, 0, 1) + ) + ) + theta_matrix = np.array( + ( + (np.cos(theta), 0, np.sin(theta)), + (0, 1, 0), + (-np.sin(theta), 0, np.cos(theta)) + ) ) - theta_matrices = np.array( - [ - [1, 0, 0], - [0, np.cos(theta), -np.sin(theta)], - [0, np.sin(theta), np.cos(theta)] - ] + psi_matrix = np.array( + ( + (1, 0, 0), + (0, np.cos(psi), -np.sin(psi)), + (0, np.sin(psi), np.cos(psi)) + ) ) - return phi_matrices * theta_matrices + rots[i] = phi_matrix @ theta_matrix @ psi_matrix + return rots @staticmethod def generate_xy_plane(n_pix): From cbd542fb8a1a0c6fffcc769ea4e7fa61e5c58be1 Mon Sep 17 00:00:00 2001 From: "deepsource-autofix[bot]" <62050782+deepsource-autofix[bot]@users.noreply.github.com> Date: Tue, 29 Mar 2022 06:52:21 +0000 Subject: [PATCH 04/14] Format code with black --- .../iterative_refinement/expectation_maximization.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index f3da42c..469f9d0 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -282,21 +282,21 @@ def grid_SO3_uniform(n_rotations): ( (np.cos(phi), -np.sin(phi), 0), (np.sin(phi), np.cos(phi), 0), - (0, 0, 1) + (0, 0, 1), ) ) theta_matrix = np.array( ( (np.cos(theta), 0, np.sin(theta)), (0, 1, 0), - (-np.sin(theta), 0, np.cos(theta)) + (-np.sin(theta), 0, np.cos(theta)), ) ) psi_matrix = np.array( ( (1, 0, 0), (0, np.cos(psi), -np.sin(psi)), - (0, np.sin(psi), np.cos(psi)) + (0, np.sin(psi), np.cos(psi)), ) ) rots[i] = phi_matrix @ theta_matrix @ psi_matrix From 319d8604ddebf389ccd39779468b8d748dc46b7c Mon Sep 17 00:00:00 2001 From: thisTyler Date: Mon, 28 Mar 2022 23:54:00 -0700 Subject: [PATCH 05/14] Removed healpy import --- reconstructSPI/iterative_refinement/expectation_maximization.py | 1 - 1 file changed, 1 deletion(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index f3da42c..18c3078 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -1,6 +1,5 @@ """Iterative refinement with Bayesian expectation maximization.""" -import healpy as hp import numpy as np from simSPI.transfer import eval_ctf From c9d0335e9cb500bca036947e6934132ea57a30d8 Mon Sep 17 00:00:00 2001 From: Tyler Heim Date: Tue, 29 Mar 2022 14:18:04 -0700 Subject: [PATCH 06/14] Untested geomstats implementation --- .../expectation_maximization.py | 29 +++---------------- 1 file changed, 4 insertions(+), 25 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index a6459e3..0143318 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -2,6 +2,7 @@ import numpy as np from simSPI.transfer import eval_ctf +from geomstats.geometry import special_orthogonal class IterativeRefinement: @@ -274,31 +275,9 @@ def grid_SO3_uniform(n_rotations): Array describing rotations. Shape (n_rotations, 3, 3) """ - rots = np.zeros((n_rotations, 3, 3)) - for i in range(n_rotations): - phi, theta, psi = np.random.random(3) * 2 * np.pi - phi_matrix = np.array( - ( - (np.cos(phi), -np.sin(phi), 0), - (np.sin(phi), np.cos(phi), 0), - (0, 0, 1), - ) - ) - theta_matrix = np.array( - ( - (np.cos(theta), 0, np.sin(theta)), - (0, 1, 0), - (-np.sin(theta), 0, np.cos(theta)), - ) - ) - psi_matrix = np.array( - ( - (1, 0, 0), - (0, np.cos(psi), -np.sin(psi)), - (0, np.sin(psi), np.cos(psi)), - ) - ) - rots[i] = phi_matrix @ theta_matrix @ psi_matrix + geom = special_orthogonal(3, 'matrix') + rots = geom.random_uniform(n_rotations) + return rots @staticmethod From c31516445110ff4e59da4741740e41eb432d774e Mon Sep 17 00:00:00 2001 From: Tyler Heim Date: Tue, 29 Mar 2022 14:18:37 -0700 Subject: [PATCH 07/14] Geomstats import --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index e52e0d7..f48fdee 100644 --- a/environment.yml +++ b/environment.yml @@ -6,6 +6,7 @@ dependencies: - codecov - coverage - gemmi + - geomstats - numpy - numba - pillow>=8.2.0 From 69c60583ef5df24cdf484c3a3cb33d2da43d3a23 Mon Sep 17 00:00:00 2001 From: "deepsource-autofix[bot]" <62050782+deepsource-autofix[bot]@users.noreply.github.com> Date: Tue, 29 Mar 2022 21:18:52 +0000 Subject: [PATCH 08/14] Format code with black --- .../iterative_refinement/expectation_maximization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index 0143318..ffd4e3c 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -275,9 +275,9 @@ def grid_SO3_uniform(n_rotations): Array describing rotations. Shape (n_rotations, 3, 3) """ - geom = special_orthogonal(3, 'matrix') + geom = special_orthogonal(3, "matrix") rots = geom.random_uniform(n_rotations) - + return rots @staticmethod From e7a649a795158b9b6fe7a4b3f24c342e3ad96729 Mon Sep 17 00:00:00 2001 From: thisTyler Date: Wed, 30 Mar 2022 00:21:12 -0700 Subject: [PATCH 09/14] Fixed rotations to include negatives to cover entire space --- .../iterative_refinement/expectation_maximization.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index ffd4e3c..6ae119a 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -275,9 +275,10 @@ def grid_SO3_uniform(n_rotations): Array describing rotations. Shape (n_rotations, 3, 3) """ - geom = special_orthogonal(3, "matrix") + geom = special_orthogonal.SpecialOrthogonal(3, "matrix") rots = geom.random_uniform(n_rotations) - + negatives = np.tile(np.random.randint(2, size=n_rotations) * 2 - 1, (3, 3, 1)).T + rots[:] *= negatives return rots @staticmethod From e353f0fce96134bff53689bd6c18c50c115a2ec8 Mon Sep 17 00:00:00 2001 From: thisTyler Date: Sat, 2 Apr 2022 11:40:23 -0700 Subject: [PATCH 10/14] Geomstats is a pip package --- environment.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index f48fdee..9e16a27 100644 --- a/environment.yml +++ b/environment.yml @@ -6,7 +6,6 @@ dependencies: - codecov - coverage - gemmi - - geomstats - numpy - numba - pillow>=8.2.0 @@ -16,4 +15,5 @@ dependencies: - pytorch - pip : - git+https://github.com/compSPI/simSPI.git + - geomstats From 701891e916eab933262757c2058a9543e3720966 Mon Sep 17 00:00:00 2001 From: thisTyler Date: Sat, 2 Apr 2022 11:45:19 -0700 Subject: [PATCH 11/14] Import changes. verifying with black. --- .../iterative_refinement/expectation_maximization.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index 6c618a4..a9952b3 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -1,8 +1,9 @@ """Iterative refinement with Bayesian expectation maximization.""" import numpy as np -from simSPI.transfer import eval_ctf + from geomstats.geometry import special_orthogonal +from simSPI.transfer import eval_ctf class IterativeRefinement: From d51f0bc230b328583226f49b7e050c587be26c47 Mon Sep 17 00:00:00 2001 From: thisTyler Date: Sat, 2 Apr 2022 11:52:45 -0700 Subject: [PATCH 12/14] Updated docstring --- .../iterative_refinement/expectation_maximization.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index a9952b3..e415333 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -281,6 +281,10 @@ def build_ctf_array(self): def grid_SO3_uniform(n_rotations): """Generate uniformly distributed rotations in SO(3). + A note on functionality - the geomstats random_uniform library only produces + rotations onto one hemisphere. So, the rotations are randomly inverted, giving + them equal probability to fall in either hemisphere. + Parameters ---------- n_rotations : int From b251e5cd4bf98619162e70a022ba2e9925683326 Mon Sep 17 00:00:00 2001 From: thisTyler Date: Sat, 2 Apr 2022 11:55:34 -0700 Subject: [PATCH 13/14] Trying pre-commit to see if it fixes the imports --- .../iterative_refinement/expectation_maximization.py | 5 ++--- tests/test_expectation_maximization.py | 6 +++--- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index e415333..ae7d57e 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -1,7 +1,6 @@ """Iterative refinement with Bayesian expectation maximization.""" import numpy as np - from geomstats.geometry import special_orthogonal from simSPI.transfer import eval_ctf @@ -205,7 +204,7 @@ def normalize_map(map_3d, counts, norm_const): Shape (n_pix, n_pix, n_pix) map normalized by counts. """ - return map_3d * counts / (norm_const + counts**2) + return map_3d * counts / (norm_const + counts ** 2) @staticmethod def apply_noise_model(map_3d_f_norm_1, map_3d_f_norm_2): @@ -360,7 +359,7 @@ def generate_slices(map_3d_f, xy_plane, n_pix, rots): map_3d_f = np.ones_like(map_3d_f) xyz_rotated = np.ones_like(xy_plane) - size = n_rotations * n_pix**2 + size = n_rotations * n_pix ** 2 slices = np.random.normal(size=size) slices = slices.reshape((n_rotations, n_pix, n_pix)) return slices, xyz_rotated diff --git a/tests/test_expectation_maximization.py b/tests/test_expectation_maximization.py index 37ee222..4f44da6 100644 --- a/tests/test_expectation_maximization.py +++ b/tests/test_expectation_maximization.py @@ -90,7 +90,7 @@ def test_grid_SO3_uniform(test_ir, n_particles): def test_generate_xy_plane(test_ir, n_pix): """Test generation of xy plane.""" xy_plane = test_ir.generate_xy_plane(n_pix) - assert xy_plane.shape == (n_pix**2, 3) + assert xy_plane.shape == (n_pix ** 2, 3) def test_generate_slices(test_ir, n_particles, n_pix): @@ -101,9 +101,9 @@ def test_generate_slices(test_ir, n_particles, n_pix): slices, xyz_rotated = test_ir.generate_slices(map_3d, xy_plane, n_pix, rots) - assert xy_plane.shape == (n_pix**2, 3) + assert xy_plane.shape == (n_pix ** 2, 3) assert slices.shape == (n_particles, n_pix, n_pix) - assert xyz_rotated.shape == (n_pix**2, 3) + assert xyz_rotated.shape == (n_pix ** 2, 3) def test_apply_ctf_to_slice(test_ir, n_pix): From f7baaa93638814c27744301b401f63040f5f981a Mon Sep 17 00:00:00 2001 From: "deepsource-autofix[bot]" <62050782+deepsource-autofix[bot]@users.noreply.github.com> Date: Sat, 2 Apr 2022 18:55:23 +0000 Subject: [PATCH 14/14] Format code with black --- .../iterative_refinement/expectation_maximization.py | 4 ++-- tests/test_expectation_maximization.py | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index ae7d57e..fecb13d 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -204,7 +204,7 @@ def normalize_map(map_3d, counts, norm_const): Shape (n_pix, n_pix, n_pix) map normalized by counts. """ - return map_3d * counts / (norm_const + counts ** 2) + return map_3d * counts / (norm_const + counts**2) @staticmethod def apply_noise_model(map_3d_f_norm_1, map_3d_f_norm_2): @@ -359,7 +359,7 @@ def generate_slices(map_3d_f, xy_plane, n_pix, rots): map_3d_f = np.ones_like(map_3d_f) xyz_rotated = np.ones_like(xy_plane) - size = n_rotations * n_pix ** 2 + size = n_rotations * n_pix**2 slices = np.random.normal(size=size) slices = slices.reshape((n_rotations, n_pix, n_pix)) return slices, xyz_rotated diff --git a/tests/test_expectation_maximization.py b/tests/test_expectation_maximization.py index 4f44da6..37ee222 100644 --- a/tests/test_expectation_maximization.py +++ b/tests/test_expectation_maximization.py @@ -90,7 +90,7 @@ def test_grid_SO3_uniform(test_ir, n_particles): def test_generate_xy_plane(test_ir, n_pix): """Test generation of xy plane.""" xy_plane = test_ir.generate_xy_plane(n_pix) - assert xy_plane.shape == (n_pix ** 2, 3) + assert xy_plane.shape == (n_pix**2, 3) def test_generate_slices(test_ir, n_particles, n_pix): @@ -101,9 +101,9 @@ def test_generate_slices(test_ir, n_particles, n_pix): slices, xyz_rotated = test_ir.generate_slices(map_3d, xy_plane, n_pix, rots) - assert xy_plane.shape == (n_pix ** 2, 3) + assert xy_plane.shape == (n_pix**2, 3) assert slices.shape == (n_particles, n_pix, n_pix) - assert xyz_rotated.shape == (n_pix ** 2, 3) + assert xyz_rotated.shape == (n_pix**2, 3) def test_apply_ctf_to_slice(test_ir, n_pix):