# -*- coding: utf-8 -*-
"""
Created on Fri Sep  9 09:59:13 2022

@author: Chris
Lab7_Q6
"""


import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import qmc
    
def region(x,y):
    return x**2+y**2<1**2

def f(x,y):
    f = 1 #x*y
    return f

def monty2d(f,region,a,b,M):
    #performs 2d Monty-Carlo integration.
    # Here:
    # f is a 2D function f(x,y)
    # region is a 2D function that returns true if x,y is in the region
    # a,b are the external boundaries for the Monty-Carlo integration
    # M gives the total number of points via N = 2**M
    
    

    #generate numbers using uniform distribution
    #N = 1000
    #p = a+(b-a)*np.random.uniform(size=[N,2])

    #generate Sobol sequence

    sampler = qmc.Sobol(d=2) # set up the sequence
    p = a+(b-a)*sampler.random_base2(M) #generate 2^M points

    xp, yp = p[:,0],p[:,1]

    plt.plot(xp,yp,'x')
    plt.gca().set_aspect('equal', adjustable='box')

    xr = xp*region(xp,yp)
    yr = yp*region(xp,yp)

    plt.plot(xr,yr,'or')
    
    N = 2**M
    fr = f(xp,yp)*region(xp,yp)
    
    I = (b-a)*(b-a)/N*sum(fr)
    
    return I
    
a,b = -1,1
M = 10
I = monty2d(f,region,a,b,M)

print(I)