#!/usr/bin/python3
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import p_roots #for pulling zeroes of Legendre orthogonal polynomials
from scipy.stats import qmc

def monty3d(f,region,a,b,M): #Calculates integral of f on arbitrary domains by MC methods
    rng = qmc.Sobol(d=3) #Sobol random number generator rng ~ U([0,1]^3)
    samp = a+(b-a)*rng.random_base2(M) #Generate 2^M vectors distributed as U([a,b]^3)
    xp,yp,zp = samp[:,0],samp[:,1],samp[:,2] #assign the random values to each basis element
    fr= f(xp,yp,zp)*region(xp,yp,zp) #Returns the function at random points, returns 0 for points outside the region
    return sum(fr)*((b-a)**3)/(2**M) #Volume (Lebesgue measure) of [a,b]^3 is (b-a)^3