Port spheroid experiment to python.

* python/spheroid.py: New file.

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()