From 9db689e734e4df13d8d069d40032e3b0a4ca767c Mon Sep 17 00:00:00 2001 From: Geoffrey Woollard Date: Sat, 2 Apr 2022 10:19:19 -0700 Subject: [PATCH 1/6] test_binary_mask_3d --- tests/test_expectation_maximization.py | 30 +++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/tests/test_expectation_maximization.py b/tests/test_expectation_maximization.py index 52af1fc..9b14649 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): @@ -198,6 +198,30 @@ def test_compute_fsc(test_ir, n_pix): assert fsc_1.shape == (n_pix // 2,) +def test_binary_mask_3d(test_ir): + """Test binary_mask_3d. + + Sums shell through an axis, then converts to circle, + then checks if circle/square ratio agrees with largest + circle inscribed in square. Should be pi/4 by trigonometry, + in the limit of infinite n_pix. Use high n_pix so good approx. + """ + n_pix = 512 + + center = (n_pix // 2, n_pix // 2, n_pix // 2) + radius = n_pix // 2 + shape = (n_pix, n_pix, n_pix) + for fill in [True, False]: + mask = test_ir.binary_mask_3d( + center, radius, shape, fill=fill, shell_thickness=1 + ) + + for axis in [0, 1, 2]: + circle = mask.sum(axis=axis) > 0 + circle_to_square_ratio = circle.mean() + assert np.isclose(circle_to_square_ratio, np.pi / 4, atol=1e-3) + + def test_expand_1d_to_3d(test_ir, n_pix): """Test expansion of 1D array into spherical shell.""" arr_1d = np.ones(n_pix // 2) From 3c99c56168090881066463d2f0c6d9b8199d93d8 Mon Sep 17 00:00:00 2001 From: "deepsource-autofix[bot]" <62050782+deepsource-autofix[bot]@users.noreply.github.com> Date: Sat, 2 Apr 2022 17:25:45 +0000 Subject: [PATCH 2/6] Format code with black --- tests/test_expectation_maximization.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_expectation_maximization.py b/tests/test_expectation_maximization.py index 9b14649..ba25a20 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 e61818c29197c415a3a6944e6f4c2bc2873137d2 Mon Sep 17 00:00:00 2001 From: Geoffrey Woollard Date: Sat, 2 Apr 2022 10:28:45 -0700 Subject: [PATCH 3/6] resolve deepsource "Unused variable found" --- .../expectation_maximization.py | 20 ++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index 4571e60..6f83445 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -141,6 +141,11 @@ def iterative_refinement(self, wiener_small_number=0.01, count_norm_const=1): ) = IterativeRefinement.compute_bayesian_weights( particles_f_1[particle_idx], slices_conv_ctfs_1, sigma ) + print( + "log z_norm_const_1={}, em_loss_1={}".format( + z_norm_const_1, em_loss_1 + ) + ) ( bayes_factors_2, z_norm_const_2, @@ -148,6 +153,11 @@ def iterative_refinement(self, wiener_small_number=0.01, count_norm_const=1): ) = IterativeRefinement.compute_bayesian_weights( particles_f_2[particle_idx], slices_conv_ctfs_2, sigma ) + print( + "log z_norm_const_2={}, em_loss_2={}".format( + z_norm_const_2, em_loss_2 + ) + ) for one_slice_idx in range(bayes_factors_1.shape[0]): xyz = xyz_rotated[one_slice_idx] @@ -213,7 +223,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): @@ -361,7 +371,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 @@ -422,7 +432,7 @@ def compute_bayesian_weights(particle, slices, sigma): ) slices_norm = np.linalg.norm(slices, axis=(1, 2)) ** 2 particle_norm = np.linalg.norm(particle) ** 2 - scale = -((2 * sigma**2) ** -1) + scale = -((2 * sigma ** 2) ** -1) log_bayesian_weights = scale * (slices_norm - 2 * corr_slices_particle) offset_safe = log_bayesian_weights.max() bayesian_weights = np.exp(log_bayesian_weights - offset_safe) @@ -540,8 +550,8 @@ def binary_mask_3d(center, radius, shape, fill=True, shell_thickness=1): a, b, c = center nx0, nx1, nx2 = shape x0, x1, x2 = np.ogrid[-a : nx0 - a, -b : nx1 - b, -c : nx2 - c] - r2 = x0**2 + x1**2 + x2**2 - mask = r2 <= radius**2 + r2 = x0 ** 2 + x1 ** 2 + x2 ** 2 + mask = r2 <= radius ** 2 if not fill and radius - shell_thickness > 0: mask_outer = mask mask_inner = r2 <= (radius - shell_thickness) ** 2 From 1958a68a3c74ac6c809ecdf5412e8a1399c240b8 Mon Sep 17 00:00:00 2001 From: "deepsource-autofix[bot]" <62050782+deepsource-autofix[bot]@users.noreply.github.com> Date: Sat, 2 Apr 2022 17:29:29 +0000 Subject: [PATCH 4/6] Format code with black --- .../iterative_refinement/expectation_maximization.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/reconstructSPI/iterative_refinement/expectation_maximization.py b/reconstructSPI/iterative_refinement/expectation_maximization.py index 6f83445..ebaa780 100644 --- a/reconstructSPI/iterative_refinement/expectation_maximization.py +++ b/reconstructSPI/iterative_refinement/expectation_maximization.py @@ -223,7 +223,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): @@ -371,7 +371,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 @@ -432,7 +432,7 @@ def compute_bayesian_weights(particle, slices, sigma): ) slices_norm = np.linalg.norm(slices, axis=(1, 2)) ** 2 particle_norm = np.linalg.norm(particle) ** 2 - scale = -((2 * sigma ** 2) ** -1) + scale = -((2 * sigma**2) ** -1) log_bayesian_weights = scale * (slices_norm - 2 * corr_slices_particle) offset_safe = log_bayesian_weights.max() bayesian_weights = np.exp(log_bayesian_weights - offset_safe) @@ -550,8 +550,8 @@ def binary_mask_3d(center, radius, shape, fill=True, shell_thickness=1): a, b, c = center nx0, nx1, nx2 = shape x0, x1, x2 = np.ogrid[-a : nx0 - a, -b : nx1 - b, -c : nx2 - c] - r2 = x0 ** 2 + x1 ** 2 + x2 ** 2 - mask = r2 <= radius ** 2 + r2 = x0**2 + x1**2 + x2**2 + mask = r2 <= radius**2 if not fill and radius - shell_thickness > 0: mask_outer = mask mask_inner = r2 <= (radius - shell_thickness) ** 2 From 5451fd14f2d887d9232f9f9f71bf8a89a42fd88e Mon Sep 17 00:00:00 2001 From: Geoffrey Woollard Date: Sat, 2 Apr 2022 10:57:12 -0700 Subject: [PATCH 5/6] more tests --- tests/test_expectation_maximization.py | 52 +++++++++++++++++++++++--- 1 file changed, 46 insertions(+), 6 deletions(-) diff --git a/tests/test_expectation_maximization.py b/tests/test_expectation_maximization.py index ba25a20..cbf4a60 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): @@ -201,10 +201,15 @@ def test_compute_fsc(test_ir, n_pix): def test_binary_mask_3d(test_ir): """Test binary_mask_3d. - Sums shell through an axis, then converts to circle, + Tests the limit of infinite n_pix. Use high n_pix so good approx. + 1. Sums shell through an axis, then converts to circle, then checks if circle/square ratio agrees with largest - circle inscribed in square. Should be pi/4 by trigonometry, - in the limit of infinite n_pix. Use high n_pix so good approx. + circle inscribed in square. Should be pi/4. + + 2. Make shells at sizes r and r/2 and check ratios of perimeter + of circle (mid slice) and surface area of sphere. + + 3. Make filled sphere of sizes r and r/2 and check ratio of volume. """ n_pix = 512 @@ -221,6 +226,41 @@ def test_binary_mask_3d(test_ir): circle_to_square_ratio = circle.mean() assert np.isclose(circle_to_square_ratio, np.pi / 4, atol=1e-3) + mask = test_ir.binary_mask_3d(center, radius, shape, fill=True, shell_thickness=1) + circle = mask[n_pix // 2] + circle_to_square_ratio = circle.mean() + assert np.isclose(circle_to_square_ratio, np.pi / 4, atol=1e-3) + + r_half = radius / 2 + for shell_thickness in [1, 2]: + mask_r = test_ir.binary_mask_3d( + center, radius, shape, fill=False, shell_thickness=1 + ) + mask_r_half = test_ir.binary_mask_3d( + center, r_half, shape, fill=False, shell_thickness=1 + ) + perimeter_ratio = mask_r[n_pix // 2].sum() / mask_r_half[n_pix // 2].sum() + assert np.isclose(2, perimeter_ratio, atol=0.1) + if shell_thickness == 1: + assert np.isclose( + mask_r[n_pix // 2].sum() / (2 * np.pi * radius), 1, atol=0.1 + ) + assert np.isclose( + mask_r_half[n_pix // 2].sum() / (2 * np.pi * r_half), 1, atol=0.1 + ) + + surface_area_ratio = mask_r.sum() / mask_r_half.sum() + surface_area_ratio_analytic = (radius / r_half) ** 2 + assert np.isclose(surface_area_ratio, surface_area_ratio_analytic, atol=0.1) + + mask_r = test_ir.binary_mask_3d(center, radius, shape, fill=True, shell_thickness=1) + mask_r_half = test_ir.binary_mask_3d( + center, r_half, shape, fill=True, shell_thickness=1 + ) + volume_ratio = mask_r.sum() / mask_r_half.sum() + volume_ratio_analytic = (radius / r_half) ** 3 + assert np.isclose(volume_ratio, volume_ratio_analytic, atol=0.005) + def test_expand_1d_to_3d(test_ir, n_pix): """Test expansion of 1D array into spherical shell.""" From a868ab712e07b6cb64610acb21309e1b4a67cbf7 Mon Sep 17 00:00:00 2001 From: "deepsource-autofix[bot]" <62050782+deepsource-autofix[bot]@users.noreply.github.com> Date: Sat, 2 Apr 2022 17:57:30 +0000 Subject: [PATCH 6/6] Format code with black --- tests/test_expectation_maximization.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_expectation_maximization.py b/tests/test_expectation_maximization.py index cbf4a60..598179f 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):