diff options
| -rw-r--r-- | python/spheroid.py | 53 |
1 files changed, 53 insertions, 0 deletions
diff --git a/python/spheroid.py b/python/spheroid.py new file mode 100644 index 0000000..3d500ab --- /dev/null +++ b/python/spheroid.py @@ -0,0 +1,53 @@ +# nsmc --- n-sphere Monte Carlo method +# Copyright © 2022 Arun I <arunisaac@systemreboot.net> +# Copyright © 2022 Murugesan Venkatapathi <murugesh@iisc.ac.in> +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +# General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see +# <https://www.gnu.org/licenses/>. + +import nsmc +from math import isclose +from numpy import linspace, zeros +import matplotlib.pyplot as plt + +rtol = 0.1 +eccentricity = 2 +dimension = linspace(10, 100, 10, dtype=int) +oracle = lambda n: nsmc.make_spheroid_extent_oracle(n, eccentricity) +true_volume = lambda n: nsmc.spheroid_true_volume(n, eccentricity) + +samples_per_cone = 1000 +solid_angle_factor = 4 +solid_angle_threshold_exponent_factor = 0.25 + +def samples_for_dimension(oracle, samples_per_cone, solid_angle_factor, solid_angle_threshold_exponent_factor, dimension): + for n in dimension: + print(n, end=' ', flush=True) + mean = zeros(n) + mean[0] = 1 + vol, samples = nsmc.volume_importance(oracle(n), mean, samples_per_cone, + solid_angle_factor, solid_angle_threshold_exponent_factor) + vol = 2 * vol + print(f'{vol} {true_volume(n)} {samples}') + if not isclose(vol, true_volume(n), rel_tol=rtol): + print(f'Volume computation failed to reach tolerance {rtol} for dimension {n}') + yield samples + +samples = list(samples_for_dimension(oracle, samples_per_cone, solid_angle_factor, solid_angle_threshold_exponent_factor, dimension)) +plt.plot(dimension, samples, '-x') +plt.yscale('log') +plt.xlabel('Dimension') +plt.ylabel('Number of samples') +plt.savefig('spheroid.png') +plt.close() |