'''
==========================================================================
Python for Parallelism in Introductory Computer Science Education
SC '13 HPC Educators Program
Steven Bogaerts, Wittenberg University
Joshua Stough, Washington and Lee
http://www.joshuastough.com/SC13
MIT License: see README_LICENSE.txt


file: parallelIntegration.py
author: stough
Summary: Timing sequential versus parallel numerical integration
of sin from 0 to 1 in N steps (specified by argument),
default 10000000. Riemann sum

Uses Pool/map paradigm, where we evenly divide the independent
work (e.g., [0,.25], (.25,.5], (.5,.75], (.75,1]) and have
each process of the pool compute the partial integrals and add
up the results.

$ python[3] parallelIntegration.py [10000000]

See:
http://docs.python.org/library/multiprocessing.html
http://docs.python.org/py3k/library/multiprocessing.html
==========================================================================
'''

from multiprocessing import Pool, Process, cpu_count
from math import sin
import time, sys

#Dependencies defined below main()

def main():
    """
    The main function, which times integration of sin in [0,1] in
    n steps (n provided by user).
    -time integrate(sin, 0,1,n).
    -time parallel integrate.  That is, evenly divide the domain and
     have each process in the pool independently compute the integral
     for that subdomain. Then sum the partial integrals.
    """
    n = 10000000
    if len(sys.argv) > 1:
        n = int(sys.argv[1])
        
    f = sin
    #First, time with one process:
    start = time.time()
    ans = integrate(f, 0,1,1.0/n)
    elapsed = time.time() - start
    print("sequential integral ans, time: %f, %f sec." \
          % (ans, elapsed))

    #So that cpu usage shows a lull.
    time.sleep(3)


    #Now, the parallel solver
    start = time.time()
    stpsize = 1.0/n
    #Use the number of cpu's to split the work.
    cpus = cpu_count()
    
    #make a list of tuples for the arguments to the processes.
    args = []
    endpoints = linspace(0,1,cpus+1)
    for i in range(cpus):
        args.append((f, endpoints[i], endpoints[i+1], stpsize))

    #Each element of args is now a tuple containing (f, a,b,stpsize)
    #where a and b are the partial domain for the integration.

    #Instantiate the pool of processes
    pool = Pool(processes = cpus)

    #Have each process in the pool execute wrapIntegrate(args[i])
    #as long as there are more args to consider.
    results = pool.map(wrapIntegrate, args)

    #The integral of the domain is the sum of the integrals of
    #the non-overlapping partial domains.
    ans = sum(results)
    
    elapsed = time.time() - start
    print("parallel integral ans, time: %f, %f sec." \
          % (ans, elapsed))


def integrate(f, a,b,h):
    """computes the integral of f from a to b in steps of size h"""
    s = 0
    x = a + h/2.0
    while x < b:
        s += h*f(x)
        x += h
    return s

def wrapIntegrate(fabh):
    """
    Wrapper for the integrate function.
    In that the function sent to map can take only one argument,
    and argument unpacking is not supported in the method header
    in python3, this is a workaround. Search 'calculatestar' in
    above web reference.
    """
    (f, a,b,h) = fabh
    return integrate(f, a,b,h)

def linspace(a,b,nsteps):
    """
    returns list of simple linear steps from a to b in nsteps.
    """
    ssize = float(b-a)/(nsteps-1)
    return [a + i*ssize for i in range(nsteps)]


#Call the main method if run from the command line.
if __name__ == '__main__':
    main()