Port integral experiments to python.

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

1 files changed, 74 insertions, 0 deletions

diff --git a/python/distribution-integrals.py b/python/distribution-integrals.py
new file mode 100644
index 0000000..95287fa
--- /dev/null
+++ b/python/distribution-integrals.py
@@ -0,0 +1,74 @@
+# 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(integrand, oracle, true_integral, n, rtol):
+    print(n, end=' ', flush=True)
+    integral, samples = nsmc.integral(integrand, oracle, true_integral(n), n, rtol)
+    return samples
+
+def average_samples_for_rtol(integrand, oracle, true_integral, rtol, trials, dimension):
+    for n in dimension:
+        yield mean_thunk(lambda: samples_for_dimension_rtol(integrand, oracle, true_integral, n, rtol), trials)
+
+def experiment(integrand, oracle, true_integral, rtols, trials, dimension, filename):
+    print(filename)
+    for rtol in rtols:
+        print(f'rtol = {rtol}:', end=' ', flush=True)
+        samples = list(average_samples_for_rtol(integrand, oracle, true_integral, 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) domain
+max_extent = 1.0
+oracle = nsmc.make_uniform_extent_oracle(0, max_extent)
+
+# Gaussian integrand
+integrand = nsmc.gaussian_integrand()
+true_integral = lambda n: nsmc.gaussian_true_integral(n, max_extent)
+experiment(integrand, oracle, true_integral, rtols, trials, dimension, 'gaussian.png')
+
+# Polynomial integrand
+coefficients = [1, -1.5, 0.6875, -0.09375]
+integrand = nsmc.polynomial_integrand(coefficients)
+true_integral = lambda n: nsmc.polynomial_true_integral(n, max_extent, coefficients)
+experiment(integrand, oracle, true_integral, rtols, trials, dimension, 'polynomial.png')
+
+# x-coordinate integrand
+integrand = nsmc.x_coordinate_integrand()
+true_integral = lambda n: nsmc.x_coordinate_true_integral(n, max_extent)
+experiment(integrand, oracle, true_integral, rtols, trials, dimension, 'x-coordinate.png')