from scipy import *
from numpy.random import *
import math
import matplotlib.pyplot as pyplot
import numpy as np
from scipy.integrate import tplquad


def expf3d(x1, x2, x3):
    return 1.0*(np.exp(x1**2+x2**2+x3**2))

def mc3d(f, N_list):
    res = []
    for N in N_list:
        res.append(sum(f(uniform(0,1,N), uniform(0,1,N),
                                uniform(0,1,N)))/N)
    return array(res)

def avg(f, N):
    res = f()
    for i in range(1,N):
        res += f()
    return res/N

##########################################

pyplot.figure()


print mc3d(expf3d, [1000000])

N_range = logspace(0, 6, 10)
#print N_range
#print mc2d(f2d, N_range)
#val, abserr = tplquad(ff, bounds_x[0], bounds_x[1], lambda y: bounds_y[0], lambda y: bounds_y[1], \
 #                     lambda y, z: bounds_z[0], lambda y, z: bounds_z[1])
val3d, abserr = tplquad(lambda x,y,z: expf3d(x, y, z), 0, 1, lambda y: 0, lambda y: 1, lambda y, z: 0, lambda y, z: 1)

#print mc3d(expf3d, N_range)

pyplot.xscale("log")
pyplot.yscale("log")
pyplot.plot(N_range, avg(lambda: abs(mc3d(expf3d, N_range) - val3d), 20), "rx", markersize=3)
#pyplot.plot(N_range, avg(lambda: abs(mc2d(f2d, N_range) - pi), 20), "bo", markersize=3)
pyplot.show()

#figure.savefig("pi.pdf")