#!/usr/bin/env python

# To execute:
#
#   mpirun -np 4 demo.py

from scipy import *
import pypar    # The Python-MPI interface
import random

def g(x):
  return 1./(1.0 - cos(x[0])*cos(x[1])*cos(x[2]))/pi**3

exact = 1.3932039296856768591842462603255


if __name__ == '__main__':
    
    mpi_size = pypar.size()
    my_rank =  pypar.rank()
    my_name =  pypar.get_processor_name()
    print "I am proc %d of %d on node %s" %(my_rank, mpi_size, my_name)
    
    random.jumpahead(my_rank)   # different start on each processor

    N = 100000
    ave=0.0
    av2=0.0
    for i in range(N):
        f = g(array([random.random(),random.random(),random.random()])*pi)
        ave += f
        av2 += f*f
  
    ave=ave/N   # now we have average
    av2=av2/N   #     and average^2
    Vol = pi**3 # Volume of the region
    Int = Vol*ave         # Integrals is Vol*<f>
    Int2 = Vol*Vol*av2    # For error we also need Vol^2 * <f^2>
    Error = sqrt((Int2-Int*Int)/N) # This is standard deviation  sigma^2 =  ( (<f>*Vol)^2 - <f^2>*Vol^2 )/N
    print "Integral=", Int, "Error=", Error, "approximation-exact=", Int-exact
    
    res=array([Int,Int2])
    sum_res=zeros(2)
    pypar.reduce(res, pypar.SUM, 0, buffer=sum_res)

    if my_rank==0:
        Int = sum_res[0]/mpi_size
        Int2 = sum_res[1]/mpi_size
        Error = sqrt((Int2-Int*Int)/(N*mpi_size))
        print "Final Integral=", Int, "Error=", Error, "approximation-exact=", Int-exact
    
    pypar.finalize()