Port volume experiments to python.

* python/distribution-volumes.py: New file.

1 files changed, 77 insertions, 0 deletions

diff --git a/python/distribution-volumes.py b/python/distribution-volumes.py
new file mode 100644
index 0000000..57e6c03
--- /dev/null
+++ b/python/distribution-volumes.py
@@ -0,0 +1,77 @@
+# 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 numpy import linspace
+from statistics import mean
+import matplotlib.pyplot as plt
+
+def mean_thunk(f, trials):
+    return mean([f() for i in range(trials)])
+
+def samples_for_dimension_rtol(oracle, true_volume, n, rtol):
+    print(n, end=' ', flush=True)
+    volume, samples = nsmc.volume(oracle(n), true_volume(n), n, rtol)
+    return samples
+
+def average_samples_for_rtol(oracle, true_volume, rtol, trials, dimension):
+    for n in dimension:
+        yield mean_thunk(lambda: samples_for_dimension_rtol(oracle, true_volume, n, rtol), trials)
+
+def experiment(oracle, true_volume, rtols, trials, dimension, filename):
+    print(filename)
+    for rtol in rtols:
+        print(f'rtol = {rtol}:', end=' ', flush=True)
+        samples = list(average_samples_for_rtol(oracle, true_volume, rtol, trials, dimension))
+        plt.plot(dimension, samples, '-x', label=f'rtol = {rtol}')
+        print()
+    plt.yscale('log')
+    plt.xlabel('Dimension')
+    plt.ylabel('Average number of samples')
+    plt.legend()
+    plt.savefig(filename)
+    plt.close()
+
+# Number of trials to average each dimension, rtol point over
+trials = 2
+
+# Dimensions and relative error tolerances to plot
+dimension = linspace(10, 100, 10, dtype=int)
+rtols = [0.2, 0.1, 0.05]
+
+# uniform(0,1)
+oracle = lambda n: nsmc.make_uniform_extent_oracle(0, 1)
+true_volume = lambda n: nsmc.uniform_true_volume(0, 1, n)
+experiment(oracle, true_volume, rtols, trials, dimension, 'uniform.png')
+
+# beta(2,2)
+oracle = lambda n: nsmc.make_beta_extent_oracle(2, 2)
+true_volume = lambda n: nsmc.beta_true_volume(2, 2, n)
+experiment(oracle, true_volume, rtols, trials, dimension, 'beta.png')
+
+# arcsine
+oracle = lambda n: nsmc.make_beta_extent_oracle(0.5, 0.5)
+true_volume = lambda n: nsmc.beta_true_volume(0.5, 0.5, n)
+experiment(oracle, true_volume, rtols, trials, dimension, 'arcsine.png')
+
+# cube (pretty slow for large dimensions and tight tolerances)
+oracle = lambda n: nsmc.make_cube_extent_oracle(n, 1.0)
+true_volume = lambda n: nsmc.cube_true_volume(n, 1.0)
+dimension = linspace(10, 40, 4, dtype=int)
+rtols = [0.2, 0.1]
+experiment(oracle, true_volume, rtols, trials, dimension, 'cube.png')