From 5acf0ac347382705674bb596440d27cba3730bac Mon Sep 17 00:00:00 2001 From: Deepak Mallubhotla Date: Sat, 17 Sep 2022 16:08:11 -0500 Subject: [PATCH] fix: Correctly generates monte carlo version of xy model dipoles --- .../log_spaced_random_choice_xy_model.py | 9 ++- ...ultiple_dipole_fixed_magnitude_xy_model.py | 60 +++++++++++++++++++ 2 files changed, 64 insertions(+), 5 deletions(-) diff --git a/pdme/model/log_spaced_random_choice_xy_model.py b/pdme/model/log_spaced_random_choice_xy_model.py index 6576bb8..4010235 100644 --- a/pdme/model/log_spaced_random_choice_xy_model.py +++ b/pdme/model/log_spaced_random_choice_xy_model.py @@ -107,15 +107,14 @@ class LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel(DipoleModel): rng = rng_to_use shape = (monte_carlo_n, self.n_max) - theta = 2 * numpy.pi * rng.random(shape) - phi = numpy.arccos(2 * rng.random(shape) - 1) + phi = 2 * numpy.pi * rng.random(shape) p_mask = rng.binomial(1, self.prob_occupancy, shape) p_magnitude = self.pfixed * p_mask - px = p_magnitude * numpy.cos(theta) * numpy.sin(phi) - py = p_magnitude * numpy.sin(theta) * numpy.sin(phi) - pz = p_magnitude * numpy.cos(phi) + px = p_magnitude * numpy.cos(phi) + py = p_magnitude * numpy.sin(phi) + pz = p_magnitude * 0 sx = rng.uniform(self.xmin, self.xmax, shape) sy = rng.uniform(self.ymin, self.ymax, shape) diff --git a/tests/model/test_log_spaced_random_count_multiple_dipole_fixed_magnitude_xy_model.py b/tests/model/test_log_spaced_random_count_multiple_dipole_fixed_magnitude_xy_model.py index 5ee402f..a885877 100644 --- a/tests/model/test_log_spaced_random_count_multiple_dipole_fixed_magnitude_xy_model.py +++ b/tests/model/test_log_spaced_random_count_multiple_dipole_fixed_magnitude_xy_model.py @@ -97,6 +97,66 @@ def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles_invariant_xy() ) +def test_random_count_multiple_dipole_fixed_mag_model_get_dipoles_invariant_monte_carlo_xy(): + + x_min = -10 + x_max = 10 + y_min = -5 + y_max = 5 + z_min = 2 + z_max = 3 + p_fixed = 10 + max_frequency = 5 + monte_carlo_n = 20 + + model = LogSpacedRandomCountMultipleDipoleFixedMagnitudeXYModel( + x_min, x_max, y_min, y_max, z_min, z_max, 0, max_frequency, p_fixed, 1, 0.5 + ) + model.rng = numpy.random.default_rng(1234) + + dipole_arrangement = model.get_dipoles(5) + dipoles = dipole_arrangement.dipoles + + assert len(dipoles) == 1, "Should have only had one dipole generated." + for i in range(10): + dipoles = model.get_monte_carlo_dipole_inputs(monte_carlo_n, max_frequency) + + min_s = numpy.array([x_min, y_min, z_min]) + max_s = numpy.array([x_max, y_max, z_max]) + + _logger.warning(dipoles) + _logger.warning(dipoles[:, 0, 0:3]) + _logger.warning(dipoles[:, 0, 3:6]) + _logger.warning(dipoles[:, 0, 6]) + + ps = dipoles[:, 0, 0:3] + ss = dipoles[:, 0, 3:6] + ws = dipoles[:, 0, 6] + + numpy.testing.assert_equal( + numpy.logical_and(min_s < ss, max_s > ss).all(), + True, + f"Dipole location [{ss}] should have been between min [{min_s}] and max [{max_s}] bounds.", + ) + assert (ws < 10**max_frequency).all() and ( + ws > 10**0 + ).all(), "Dipole frequency should have been between 0 and max." + + norms = numpy.linalg.norm(ps, axis=1) + filtered_norms = norms[norms > 0] + numpy.testing.assert_allclose( + filtered_norms, + p_fixed, + err_msg="Should have had the expected dipole moment magnitude.", + ) + + numpy.testing.assert_allclose( + ps[:, 2], + 0, + err_msg="Should have had zero z magnitude.", + ) + + def test_random_count_multiple_dipole_shape(): x_min = -10