# rerun Fernando's script to load in css
%run talktools

import random
from math import exp, sqrt, pi

def f(u):
    return(exp(-u**2/2)/sqrt(2*pi))

def loopy_monte(N):
    hits = 0
    for count in range(N):
        x = random.uniform(-5,5)
        y = random.uniform(0,0.5)
        hits = hits + (y < f(x))
        
    estimate = hits / float(N) * (0.5*10)
    return(estimate)

loopy_monte(10**6)

%timeit loopy_monte(10**6)

import numpy as np

def numpy_f(u):
    return(np.exp(-u**2/2)/np.sqrt(2*np.pi))

def numpy_monte(N):
    x = np.random.uniform(-5,5,N)
    y = np.random.uniform(0,0.5,N)
    hits = np.sum(y<numpy_f(x))

    estimate = hits / np.float(N) * (0.5*10)
    return(estimate)

numpy_monte(10**6)

%timeit numpy_monte(10**6)

from IPython import parallel

clients = parallel.Client()
clients.block = True
print(clients.ids)

dview = clients[:] # use all engines

import numpy as np

def numpy_f(u):
    return(np.exp(-u**2/2)/np.sqrt(2*np.pi))

# Redefine numpy_monte to take a seed argument
def numpy_monte(N,seed):
    np.random.seed(seed)
    x = np.random.uniform(-5,5,N)
    y = np.random.uniform(0,0.5,N)
    hits = np.sum(y<numpy_f(x))

    estimate = hits / float(N) * (0.5*10)
    return(estimate)

# First, we do this in serial
# Set up
cores = 2

# This contains the number of iterations to be performed by each engine.
# 10 times the number we did earlier. When working in parallel, you need to ensure that each engine has sufficent work to do to
# ensure that communications overheads don't swamp the benefits
N = [10**7]*cores

seeds = range(0,cores) # A list of seeds, one for each engine

# We use the map function to map numpy_monte over the lists created above
# The total number of iterations will be cores * N
%timeit estimate = sum(list(map(numpy_monte,N,seeds)))/cores
estimate = sum(list(map(numpy_monte,N,seeds)))/cores

#Display the esimate
estimate

%%px
#Define all functions on the engines

import numpy as np

def numpy_f(u):
    return(np.exp(-u**2/2)/np.sqrt(2*np.pi))

# Redefine numpy_monte to take a seed argument
def numpy_monte(N,seed):
    np.random.seed(seed)
    x = np.random.uniform(-5,5,N)
    y = np.random.uniform(0,0.5,N)
    hits = np.sum(y<numpy_f(x))

    estimate = hits / float(N) * (0.5*10)
    return(estimate)

#Do the work in parallel using the dview we created earlier
%timeit estimate = sum(dview.map(numpy_monte,N,seeds))/cores

#Print the result
estimate

%load_ext pymatbridge

%%matlab 
a = linspace(-4,4,100);
b = sin(a);
plot(a,b)

%%matlab
format long

N = 10^6;
hits = 0;
tic
for count = 1:N
    x = rand() * 10 - 5;
    y = rand() * 0.5;
    hits = hits + ( y < (exp(-x^2/2) / sqrt(2*pi)) ); 
end
estimate = hits / N * (0.5*10)
toc



%%writefile loopy_monte.m
function estimate = loopy_monte(N)

hits = 0;
tic
for count = 1:N
    x = rand() * 10 - 5;
    y = rand() * 0.5;
    hits = hits + ( y < (exp(-x^2/2) / sqrt(2*pi)) ); 
end
estimate = hits / N * (0.5*10);
toc

from pymatbridge import Matlab

# Start a MATLAB session
mlab = Matlab()
mlab.start()

# Change MATLAB's working directory to match Ipython's current working directory
current_dir = %pwd
mlab.run_code('cd ' + current_dir)

# Run the code in MATLAB and print the results
results = mlab.run_code('loopy_monte(10^6)')
print(results['content']['stdout'])

# Stop the MATLAB session
mlab.stop()


%%writefile vector_monte.m
function estimate = vector_monte(N)
% A vectorised version of our monte carlo code
tic
x = rand(1,N) * 10 - 5;
y = rand(1,N) * 0.5;
hits = ( sum(y < (exp(-x.^2/2) / sqrt(2*pi))) ); 
estimate = hits / N * (0.5*10);
toc

from pymatbridge import Matlab

# Start a MATLAB session
mlab = Matlab()
mlab.start()

# Change MATLAB's working directory to match Ipython's current working directory
current_dir = %pwd
mlab.run_code('cd ' + current_dir)

# Run the code in MATLAB and print the results
results = mlab.run_code('vector_monte(10^6)')
print(results['content']['stdout'])

# Stop the MATLAB session
mlab.stop()

%%writefile serial_monte.m
function estimate = par_monte(N,jobs)
tic
estimate = 0;
for i=1:jobs
    x = rand(1,N) * 10 - 5;
    y = rand(1,N) * 0.5;
    hits = ( sum(y < (exp(-x.^2/2) / sqrt(2*pi))) ); 
    estimate = estimate + (hits / N * (0.5*10));
end
estimate = estimate / jobs
toc

from pymatbridge import Matlab

# Start a MATLAB session
mlab = Matlab()
mlab.start()

# Change MATLAB's working directory to match Ipython's current working directory
current_dir = %pwd
mlab.run_code('cd ' + current_dir)

# Run the code in MATLAB and print the results
results = mlab.run_code('serial_monte(5*10^6,100)')
print(results['content']['stdout'])

# Stop the MATLAB session
mlab.stop()

%%writefile par_monte.m
function estimate = par_monte(N,jobs)
tic
estimate = 0;
parfor i=1:jobs
    x = rand(1,N) * 10 - 5;
    y = rand(1,N) * 0.5;
    hits = ( sum(y < (exp(-x.^2/2) / sqrt(2*pi))) ); 
    estimate = estimate + (hits / N * (0.5*10));
end
estimate = estimate / jobs
toc

from pymatbridge import Matlab

# Start a MATLAB session
mlab = Matlab()
mlab.start()

# Change MATLAB's working directory to match Ipython's current working directory
current_dir = %pwd
mlab.run_code('cd ' + current_dir)
# Open the parallel pool
mlab.run_code('parpool local')

# Run the code in MATLAB and print the results
results = mlab.run_code('par_monte(5*10^6,100)')
print(results['content']['stdout'])

# Stop the MATLAB session
mlab.stop()

%load_ext rmagic

%%R

f = function(u) exp(-u^2/2)/sqrt(2*pi)
x = seq(-5, 5, length.out = 100)
plot(x, f(x), type="l", ylim=c(0, 0.5))
abline(h=0, lty=4, col=4);abline(h=0.5, lty=3, col=4); abline(v=c(-5, 5), lty=4, col=4)

%%R -o estimate
# start the clock
start_time = proc.time()
f = function(u) exp(-u^2/2)/sqrt(2*pi)

simulate_pt = function(i){
    x = runif(1,-5,5); y = runif(1,0,0.5);
    y < f(x)
}

N = 10^6; hits = 0
for(i in 1:N)
 hits = hits + simulate_pt() 
(estimate = hits/N*(0.5*10))

#Stop the clock
print(proc.time() - start_time)
print(estimate)


%%writefile loop_monte.R
# start the clock
start_time = proc.time()
f = function(u) exp(-u^2/2)/sqrt(2*pi)

simulate_pt = function(i){
    x = runif(1,-5,5); y = runif(1,0,0.5);
    y < f(x)
}

N = 10^6; hits = 0
for(i in 1:N)
 hits = hits + simulate_pt() 
(estimate = hits/N*(0.5*10))

#Stop the clock
print(proc.time() - start_time)
print(estimate)

!R --vanilla < loop_monte.R

%%writefile vector_monte.R
start_time = proc.time()
N=10^6
f = function(u) exp(-u^2/2)/sqrt(2*pi)
estimate = sum(runif(N, 0, 0.5) < f(runif(N, -5, 5)))/N*(0.5*10)

print(proc.time() - start_time)
print(estimate)

!R --vanilla < vector_monte.R