Added example for parallel Monte-Carlson algorithm using MT2203 stream

Author: oleksandr-pavlyk

example/README.md

@@ -0,0 +1,26 @@
+# Parallel Monte-Carlo example
+
+Using `mkl_random` package, we use MT-2203 family of pseudo-random number generation algorithms,
+we create workers, assign them RandomState objects with different members of the family of algorithms,
+and use multiprocessing Pool to distribute chunks of MC work to them to process.
+
+Each worker gets `rs` and `n` arguments, `rs` representing RandomState object associated with the worker,
+and `n` being the size of the problem. `rs` is used to generate samples of size `n`, perform Monte-Carlo
+estimate(s) based on the sample and return.
+
+After run is complete, a generator is returns that contains results of each worker. 
+
+This data is post-processed as necessary for the application.
+
+## Stick triangle problem
+
+Code is tested to estimate the probability that 3 segments, obtained by splitting a unit stick 
+in two randomly chosen places, can be sides of a triangle. This probability is known in closed form to be $\frac{1}{4}$.
+
+## Stick tetrahedron problem
+
+Code is used to estimate the probability that 6 segments, obtained by splitting a unit stick in 
+5 random chosen places, can be sides of a tetrahedron. 
+
+The probability is not known in closed form. See
+[math.stackexchange.com/questions/351913](https://math.stackexchange.com/questions/351913/probability-that-a-stick-randomly-broken-in-five-places-can-form-a-tetrahedron) for more details.

example/arg_parsing.py

@@ -0,0 +1,36 @@
+import argparse
+
+__all__ = ['parse_arguments']
+
+def pos_int(s):
+    v = int(s)
+    if v > 0:
+        return v
+    else:
+        raise argparse.ArgumentTypeError('%r is not a positive integer' % s)
+
+
+def nonneg_int(s):
+    v = int(s)
+    if v >= 0:
+        return v
+    else:
+        raise argparse.ArgumentTypeError('%r is not a non-negative integer' % s)
+
+
+def parse_arguments():
+    argParser = argparse.ArgumentParser(
+        prog="stick_tetrahedron.py",
+        description="Monte-Carlo estimation of probability that 6 segments of a stick randomly broken in 5 places can form a tetrahedron.",
+        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+
+    argParser.add_argument('-s', '--seed',        default=7777,   type=pos_int,    help="Random seed to initialize algorithms from MT2203 family")
+    argParser.add_argument('-b', '--batch_size',  default=65536,  type=pos_int,    help="Batch size for the Monte-Carlo run")
+    argParser.add_argument('-n', '--batch_count', default=2048,   type=pos_int,    help="Number of batches executed in parallel")
+    argParser.add_argument('-p', '--processes',   default=-1,     type=int,        help="Number of processes used to execute batches")
+    argParser.add_argument('-d', '--id_offset',   default=0,      type=nonneg_int, help="Offset for the MT2203/WH algorithms id")
+    argParser.add_argument('-j', '--jump_size',   default=0,      type=nonneg_int, help="Jump size for skip-ahead")
+      
+    args = argParser.parse_args()
+
+    return args

example/fancy.py

@@ -0,0 +1,90 @@
+import numpy as np
+import mkl_random as rnd
+
+__doc__ = """
+Let's solve a classic problem of MC-estimating a probability that 3 segments of a unit stick randomly broken in 2 places can form a triangle. 
+Let $u_1$ and $u_2$ be standard uniform random variables, denoting positions where the stick has been broken.
+
+Let $w_1 = \min(u_1, u_2)$ and $w_2 = \max(u_1, u_2)$. Then, length of segments are $x_1 = w_1$, $x_2 = w_2-w_1$, $x_3 = 1-w_2$. 
+These lengths must satisfy triangle inequality.
+
+The closed form result is known to be $\frac{1}{4}$.
+
+"""
+
+def triangle_inequality(x1, x2, x3):
+    """Efficiently finds `np.less(x1,x2+x3)*np.less(x2,x1+x3)*np.less(x3,x1+x2)`"""
+    tmp_sum = x2 + x3
+    res = np.less(x1, tmp_sum)   # x1 < x2 + x3
+    np.add(x1, x3, out=tmp_sum)
+    buf = np.less(x2, tmp_sum)   # x2 < x1 + x3
+    np.logical_and(res, buf, out=res)
+    np.add(x1, x2, out=tmp_sum)
+    np.less(x3, tmp_sum, out=buf) # x3 < x1 + x2
+    np.logical_and(res, buf, out=res)
+    return res
+
+
+def mc_dist(rs, n):
+    """Monte Carlo estimate of probability on sample of size `n`, using given random state object `rs`"""
+    ws = np.sort(rs.rand(2,n), axis=0)
+    x2 = np.empty(n, dtype=np.double)
+    x3 = np.empty(n, dtype=np.double)
+
+    x1 = ws[0]
+    np.subtract(ws[1], ws[0], out=x2)
+    np.subtract(1, ws[1], out=x3)
+    mc_prob = triangle_inequality(x1, x2, x3).sum() / n
+
+    return mc_prob
+
+
+def assign_worker_rs(w_rs):
+    """Assign process local random state variable `rs` the given value"""
+    assert not 'rs' in globals(), "Here comes trouble. Process is not expected to have global variable `rs`"
+
+    global rs
+    rs = w_rs
+    # wait to ensure that the assignment takes place for each worker
+    b.wait()
+
+
+def worker_compute(w_id):
+    return mc_dist(rs, batch_size)
+
+
+if __name__ == '__main__':
+    import multiprocessing as mp
+    from itertools import repeat
+    from timeit import default_timer as timer
+
+    seed = 77777
+    n_workers = 12
+    batch_size = 1024 * 256
+    batches = 10000
+
+    t0 = timer()
+    # Create instances of RandomState for each worker process from MT2203 family of generators
+    rss = [ rnd.RandomState(seed, brng=('MT2203', idx)) for idx in range(n_workers) ]
+    # use of Barrier ensures that every worker gets one
+    b = mp.Barrier(n_workers)
+
+    with mp.Pool(processes=n_workers) as pool:
+        # map over every worker once to distribute RandomState instances
+        pool.map(assign_worker_rs, rss, chunksize=1)
+        # Perform computations on workers
+        r = pool.map(worker_compute, range(batches), chunksize=1)
+
+    # retrieve values of estimates into numpy array
+    ps = np.fromiter(r, dtype=np.double)
+    # compute sample estimator's mean and standard deviation
+    p_est = ps.mean()
+    pop_std = ps.std()
+    t1 = timer()
+
+    dig = 3 - int(np.log10(pop_std))
+    frm_str = "{0:0." + str(dig) + "f}"
+    print(("Monte-Carlo estimate of probability: " + frm_str).format(p_est))
+    print(("Population estimate of the estimator's standard deviation: " + frm_str).format(pop_std))
+    print(("Expected standard deviation of the estimator: " + frm_str).format(np.sqrt(p_est * (1-p_est)/batch_size)))
+    print("Execution time: {0:0.3f} seconds".format(t1-t0))

example/parallel_mc.py

@@ -0,0 +1,44 @@
+import multiprocessing as mp
+
+__all__ = ['parallel_mc_run']
+
+def worker_compute(w_id):
+    "Worker function executed on the spawned slave process"
+    # global _local_rs
+    return _worker_mc_compute(_local_rs)
+
+
+def assign_worker_rs(w_rs):
+    """Assign process local random state variable `rs` the given value"""
+    assert not '_local_rs' in globals(), "Here comes trouble. Process is not expected to have global variable `_local_rs`"
+
+    global _local_rs
+    _local_rs = w_rs
+    # wait to ensure that the assignment takes place for each worker
+    b.wait()
+
+def parallel_mc_run(random_states, n_workers, n_batches, mc_func):
+    """
+    Given iterable `random_states` of length `n_workers`, the number of batches `n_batches`,
+    and the function `worker_compute` to execute, return iterator with results returned by 
+    the supplied function. The function is expected to conform to signature f(worker_id), 
+    and has access to worker-local global variable `rs`, containing worker's random states.
+    """
+    # use of Barrier ensures that every worker gets one
+    global b, _worker_mc_compute
+    b = mp.Barrier(n_workers)
+    
+    _worker_mc_compute = mc_func
+    with mp.Pool(processes=n_workers) as pool:
+        # 1. map over every worker once to distribute RandomState instances
+        pool.map(assign_worker_rs, random_states, chunksize=1)
+        # 2. Perform computations on workers
+        r = pool.map(worker_compute, range(n_batches), chunksize=1)
+
+    return r
+
+
+def sequential_mc_run(random_states, n_workers, n_batches, mc_func):
+    for rs in random_states:
+        for _ in range(n_batches):
+            yield mc_func(rs)

example/parallel_random_states.py

@@ -0,0 +1,17 @@
+import mkl_random as rnd
+
+
+def build_MT2203_random_states(seed, id0, n_workers):
+    # Create instances of RandomState for each worker process from MT2203 family of generators
+    return (rnd.RandomState(seed, brng=('MT2203', id0 + idx)) for idx in range(n_workers))
+
+
+def build_SFMT19937_random_states(seed, jump_size, n_workers):
+    import copy
+    # Create instances of RandomState for each worker process from MT2203 family of generators
+    rs = rnd.RandomState(seed, brng='SFMT19937')
+    yield copy.copy(rs)
+    for _ in range(1, n_workers):
+        rs.skipahead(jump_size)
+        yield copy.copy(rs)
+

example/stick_tetrahedron.py

@@ -0,0 +1,75 @@
+import numpy as np
+from parallel_mc import parallel_mc_run, sequential_mc_run
+from parallel_random_states import build_MT2203_random_states
+from sticky_math import mc_six_piece_stick_tetrahedron_prob
+from arg_parsing import parse_arguments
+
+def mc_runner(rs):
+    return mc_six_piece_stick_tetrahedron_prob(rs, batch_size)
+
+def aggregate_mc_counts(counts, n_batches, batch_size):
+    ps = counts / batch_size
+    # compute sample estimator's mean and standard deviation
+    p_est = ps.mean()
+    p_std = ps.std()/np.sqrt(batches)
+    
+    # compute parameters for Baysean posterior of the probability
+    event_count = 0
+    nonevent_count = 0
+    for ni in counts:
+        event_count += int(ni)
+        nonevent_count += int(batch_size - ni)
+
+    assert event_count >= 0
+    assert nonevent_count >= 0
+    return (p_est, p_std, event_count, nonevent_count) 
+
+
+def print_result(p_est, p_std, mc_size):
+    dig = 3 - int(np.log10(p_std)) # only show 3 digits past width of confidence interval
+    frm_str = "{0:0." + str(dig) + "f}"
+
+    print(("Monte-Carlo estimate of probability: " + frm_str).format(p_est))
+    print(("Population estimate of the estimator's standard deviation: " + frm_str).format(p_std))
+    print(("Expected standard deviation of the estimator: " + frm_str).format(np.sqrt(p_est * (1-p_est)/mc_size)))
+    print("Total MC size: {}".format(mc_size))
+    
+
+if __name__ == '__main__':
+    import multiprocessing as mp
+    from itertools import repeat
+    from timeit import default_timer as timer
+    import sys
+
+    args = parse_arguments()
+    
+    seed = args.seed
+    n_workers = args.processes
+    if n_workers <= 0:
+        n_workers = mp.cpu_count()
+
+    batch_size = args.batch_size
+    batches = args.batch_count
+    id0 = args.id_offset
+
+    t0 = timer()
+
+    rss = build_MT2203_random_states(seed, id0, n_workers)
+    r = parallel_mc_run(rss, n_workers, batches, mc_runner)
+    # r = sequential_mc_run(rss, n_workers, batches, mc_runner)
+
+    # retrieve values of estimates into numpy array
+    counts = np.fromiter(r, dtype=np.double)
+    p_est, p_std, event_count, nonevent_count = aggregate_mc_counts(counts, batches, batch_size)
+
+    t1 = timer()
+
+
+    print("Input parameters: -s {seed} -b {batchSize} -n {numBatches} -p {processes} -d {idOffset}".format(
+        seed=args.seed, batchSize=args.batch_size, numBatches=args.batch_count, processes=n_workers, idOffset=args.id_offset))
+    print("")
+    print_result(p_est, p_std, batches * batch_size)
+    print("")
+    print("Bayesian posterior beta distribution parameters: ({0}, {1})".format(event_count, nonevent_count))
+    print("")
+    print("Execution time: {0:0.3f} seconds".format(t1-t0))