Parallelize examples with joblib (closes #5)

examples/README.md

@@ -2,6 +2,25 @@
 
 These are some example scripts to demonstrate the various simulations that can be done, and to verify the simulator by reproducing results of already-published works.
 
+Scripts double as tests: each one that has expected output defines a
+`reference_table` (the values its computed `table` is checked against, in the
+same order as the script's columns) and a `tolerance` (the absolute
+comparison tolerance).  `tests/test_examples.py` runs each script and checks
+`table` against `reference_table`.  Where a script reproduces the published
+results, the reference values come from the paper; where it does not (see
+issues #88 and #91), the reference is the script's "Typical result" as a
+regression guard until the discrepancy is fixed.
+
+## Parallel Monte Carlo
+
+Most scripts that run many elections use [joblib](https://joblib.readthedocs.io/) to parallelize embarrassingly parallel simulation loops. The usual pattern is:
+
+- **`batch_size = 100`** — each worker runs this many elections before returning, which cuts down on scheduling overhead.
+- **`Parallel(n_jobs=-3, verbose=5, backend='loky')`** — use all but two CPU cores, print progress, and pin the loky backend so each worker draws elections with its own RNG state (fork-based backends would duplicate the shared RNG across workers).
+- **Aggregation** — workers return partial `Counter` or `defaultdict` results that the main process merges (order does not matter).
+
+These are Monte Carlo verification scripts: they do not fix a global random seed, so exact counts differ from run to run; each worker draws elections independently.
+
 ## Wikipedia
 
 ### Likelihood of a Condorcet cycle

examples/merrill_1984_fig_2c_2d.py

@@ -41,11 +41,11 @@
 running as many simulations, however the Coombs results are consistently
 high.
 """
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -59,6 +59,11 @@
 corr = 0.5
 D = 2
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Runoff': runoff, 'Hare': irv,
                   'Borda': borda, 'Coombs': coombs, 'Black': black}
 
@@ -87,14 +92,13 @@
     'Plurality': {2: 100.0, 3: 51.3,  4: 36.2,  5: 21.0,  7: 7.8},
 }
 
-for fig, disp, ymin, orig in (('2.c', 1.0, 50, merrill_fig_2c),
-                              ('2.d', 0.5, 0, merrill_fig_2d)):
+table = {}
 
+def simulate_batch(disp):
+    """Run one batch of elections and return the partial tallies."""
     condorcet_winner_count = {key: Counter() for key in (
         ranked_methods.keys() | rated_methods.keys() | {'CW'})}
-    start_time = time.monotonic()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         for n_cands in n_cands_list:
             v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                      disp=disp)
@@ -115,9 +119,22 @@
                     if method(utilities, tiebreaker='random') == CW:
                         condorcet_winner_count[name][n_cands] += 1
 
-    elapsed_time = time.monotonic() - start_time
-    print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)),
-          '\n')
+    return condorcet_winner_count
+
+
+for fig, disp, ymin, orig in (('2.c', 1.0, 50, merrill_fig_2c),
+                              ('2.d', 0.5, 0, merrill_fig_2d)):
+
+    jobs = [delayed(simulate_batch)(disp)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results
+    condorcet_winner_count = {key: Counter() for key in (
+        ranked_methods.keys() | rated_methods.keys() | {'CW'})}
+    for result in results:
+        for method, counter in result.items():
+            condorcet_winner_count[method].update(counter)
 
     plt.figure(f'Figure {fig}. {n_voters} voters, {n_elections} elections')
     plt.title(f'Figure {fig}: Condorcet Efficiency under Spatial-Model '
@@ -128,24 +145,25 @@
 
     for method in ('Black', 'Coombs', 'Borda', 'Approval', 'Hare', 'Runoff',
                    'Plurality'):
         x, y = zip(*sorted(orig[method].items()))
         plt.plot(x, y, ':', lw=0.8)
 
     # Restart color cycle, so result colors match
     plt.gca().set_prop_cycle(None)
 
-    table = []
+    table[fig] = []
 
     # Of those elections with CW, likelihood that method chooses CW
     x_cw, y_cw = zip(*sorted(condorcet_winner_count['CW'].items()))
     for method in ('Black', 'Coombs', 'Borda', 'Approval', 'Hare', 'Runoff',
                    'Plurality'):
         x, y = zip(*sorted(condorcet_winner_count[method].items()))
         CE = np.array(y)/y_cw
         plt.plot(x, CE*100, '-', label=method)
-        table.append([method, *CE*100])
+        table[fig].append([method, *CE*100])
 
-    print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
+    print(tabulate(table[fig], ["Method", *x], tablefmt="pipe",
+                   floatfmt='.1f'))
     print()
 
     plt.plot([], [], 'k:', lw=0.8, label='Merrill')  # Dummy plot for label
@@ -155,3 +173,32 @@
     plt.ylim(ymin, 102)
     plt.xlim(1.8, 7.2)
     plt.show()
+
+# Regression reference from the "Results with 500_000 elections" tables in the
+# docstring, ordered by n_cands_list.  These sims do not reproduce Merrill's
+# published figures exactly (discrepancies up to ~7%, see issue #88), so this
+# is only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    '2.c': {
+        'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 99.4, 98.6, 97.8, 96.9, 96.0),
+        'Borda': (100.0, 91.4, 89.2, 87.1, 85.7, 84.6),
+        'Approval': (100.0, 85.9, 79.8, 73.9, 70.1, 66.8),
+        'Hare': (100.0, 94.1, 86.6, 78.9, 71.7, 65.2),
+        'Runoff': (100.0, 94.1, 87.1, 79.7, 72.8, 66.1),
+        'Plurality': (100.0, 80.6, 67.6, 57.4, 49.3, 42.6),
+    },
+    '2.d': {
+        'Black': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 98.2, 95.9, 93.4, 90.9, 88.4),
+        'Borda': (100.0, 89.2, 86.3, 83.8, 82.1, 80.8),
+        'Approval': (100.0, 84.0, 76.9, 71.5, 67.8, 64.7),
+        'Hare': (100.0, 72.2, 50.3, 35.8, 26.0, 19.7),
+        'Runoff': (100.0, 72.2, 50.6, 35.3, 24.4, 16.9),
+        'Plurality': (100.0, 55.9, 34.7, 21.5, 13.5, 8.5),
+    },
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0

examples/merrill_1984_fig_2c_2d_updated.py

@@ -45,11 +45,11 @@
 running as many simulations, however the Coombs results are consistently
 high.
 """
-import time
 from collections import Counter
 
 import matplotlib.pyplot as plt
 import numpy as np
+from joblib import Parallel, delayed
 from tabulate import tabulate
 
 from elsim.elections import normal_electorate, normed_dist_utilities
@@ -64,6 +64,11 @@
 corr = 0.5
 D = 2
 
+# Simulate more than just one election per worker to improve efficiency
+batch_size = 100
+n_batches = n_elections // batch_size
+assert n_batches * batch_size == n_elections
+
 ranked_methods = {'Plurality': fptp, 'Top-2 Runoff': runoff, 'Hare RCV': irv,
                   'Borda': borda, 'Coombs': coombs, 'Condorcet RCV': black}
 
@@ -76,14 +81,13 @@
                      star(honest_normed_scores(utilities, 5), tiebreaker),
                  }
 
-for fig, disp, ymin in (('2.c', 1.0, 50),
-                        ('2.d', 0.5, 0)):
+table = {}
 
+def simulate_batch(disp):
+    """Run one batch of elections and return the partial tallies."""
     condorcet_winner_count = {key: Counter() for key in (
         ranked_methods.keys() | rated_methods.keys() | {'CW'})}
-    start_time = time.monotonic()
-
-    for _iteration in range(n_elections):
+    for _iteration in range(batch_size):
         for n_cands in n_cands_list:
             v, c = normal_electorate(n_voters, n_cands, dims=D, corr=corr,
                                      disp=disp)
@@ -104,27 +108,41 @@
                     if method(utilities, tiebreaker='random') == CW:
                         condorcet_winner_count[name][n_cands] += 1
 
-    elapsed_time = time.monotonic() - start_time
-    print('Elapsed:', time.strftime("%H:%M:%S", time.gmtime(elapsed_time)),
-          '\n')
+    return condorcet_winner_count
+
+
+for fig, disp, ymin in (('2.c', 1.0, 50),
+                        ('2.d', 0.5, 0)):
+
+    jobs = [delayed(simulate_batch)(disp)] * n_batches
+    print(f'{len(jobs)} tasks total:')
+    results = Parallel(n_jobs=-3, verbose=5, backend='loky')(jobs)
+
+    # Aggregate results
+    condorcet_winner_count = {key: Counter() for key in (
+        ranked_methods.keys() | rated_methods.keys() | {'CW'})}
+    for result in results:
+        for method, counter in result.items():
+            condorcet_winner_count[method].update(counter)
 
     plt.figure(f'Figure {fig}. {n_voters} voters, {n_elections} elections',
                figsize=(8, 6.5))
     plt.title(f'Figure {fig}: Condorcet Efficiency under Spatial-Model '
               f'Assumptions [Disp: {disp}]')
 
-    table = []
+    table[fig] = []
 
     # Of those elections with CW, likelihood that method chooses CW
     x_cw, y_cw = zip(*sorted(condorcet_winner_count['CW'].items()))
     for method in ('Condorcet RCV', 'Coombs', 'STAR', 'Borda', 'Score',
                    'Approval (opt.)', 'Hare RCV', 'Top-2 Runoff', 'Plurality'):
         x, y = zip(*sorted(condorcet_winner_count[method].items()))
         CE = np.array(y)/y_cw
         plt.plot(x, CE*100, '-', label=method)
-        table.append([method, *CE*100])
+        table[fig].append([method, *CE*100])
 
-    print(tabulate(table, ["Method", *x], tablefmt="pipe", floatfmt='.1f'))
+    print(tabulate(table[fig], ["Method", *x], tablefmt="pipe",
+                   floatfmt='.1f'))
     print()
 
     plt.legend()
@@ -133,3 +151,36 @@
     plt.ylim(ymin, 102)
     plt.xlim(1.8, 7.2)
     plt.show()
+
+# Regression reference from the "Results with 100_000 elections" tables in the
+# docstring, ordered by n_cands_list.  These sims do not reproduce Merrill's
+# published figures exactly (discrepancies up to ~7%, see issue #88), so this
+# is only a regression guard against breaking the current close-enough output
+# until that discrepancy is fixed.
+reference_table = {
+    '2.c': {
+        'Condorcet RCV': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 99.4, 98.6, 97.8, 96.9, 96.0),
+        'STAR': (100.0, 97.8, 94.7, 92.1, 89.9, 88.0),
+        'Borda': (100.0, 91.4, 89.2, 87.1, 85.8, 84.8),
+        'Score': (100.0, 88.7, 84.6, 82.7, 81.3, 79.9),
+        'Approval (opt.)': (100.0, 86.0, 79.7, 73.9, 70.5, 67.0),
+        'Hare RCV': (100.0, 94.1, 86.6, 79.0, 71.3, 65.1),
+        'Top-2 Runoff': (100.0, 94.1, 87.1, 79.9, 72.8, 65.9),
+        'Plurality': (100.0, 80.6, 67.8, 57.3, 49.1, 42.7),
+    },
+    '2.d': {
+        'Condorcet RCV': (100.0, 100.0, 100.0, 100.0, 100.0, 100.0),
+        'Coombs': (100.0, 98.3, 95.7, 93.4, 90.8, 88.5),
+        'STAR': (100.0, 96.5, 91.5, 87.2, 83.1, 79.8),
+        'Borda': (100.0, 89.2, 86.4, 83.9, 82.2, 80.7),
+        'Score': (100.0, 86.2, 80.5, 77.4, 74.5, 72.1),
+        'Approval (opt.)': (100.0, 83.7, 76.9, 71.6, 67.6, 64.6),
+        'Hare RCV': (100.0, 72.4, 50.4, 35.6, 26.2, 19.8),
+        'Top-2 Runoff': (100.0, 72.4, 50.5, 35.3, 24.5, 17.1),
+        'Plurality': (100.0, 56.3, 34.8, 21.5, 13.5, 8.6),
+    },
+}
+
+# Absolute tolerance (percentage points) for test_examples.py
+tolerance = 5.0