#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Sep 17 16:11:22 2024

@author: sam
"""

import matplotlib.pyplot as plt
from scipy.stats import qmc

def f(x,y):
    return 1

def region(x,y):
    return (x**2 + y**2) <= 1

def monty2d(f,region,a,b,M):
    sampler = qmc.Sobol(d=2)
    U = a+(b-a)*sampler.random_base2(M)
    xp, yp = U[:,0],U[:,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)**2/N)*sum(fr)
    return I



print(monty2d(f,region,-1,1,10))