#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Oct  4 15:59:00 2024

@author: sam
"""

from scipy.stats import qmc



def monty3d(f,region,a,b,M):
# defines function monty3d to perform 3d Monte-Carlo integration
# with arguments, a 3 variable function f to integrate over
# a 3 variable function region which defines the region of random number sampling for Monte-Carlo integration
# a and b, start and end points for the range of the sampled random numbers
# M number of random points N to generate where N = 2^M
    
    sampler = qmc.Sobol(d=3)
    # sets up the sequence
    
    U = a+(b-a)*sampler.random_base2(M)
    # a+(b-a) shifts the sampled random numbers to cover the range [a,b]
    # sampler.random_base2(M) generates 2^M 3 dimensional points
    
    xp,yp,zp = U[:,0],U[:,1],U[:,2]
    # splits random numbers into 3 parts, to call on each dimension individually
    # this is to be able to use each component for each 3 dimension input for functions f and region
    
    fr = f(xp,yp,zp)*region(xp,yp,zp)
    # multiplies the value of f at each random number by region at that point to check if it is in the region
    # which returns the value of each f(xp,yp,zp) if it is within the region, and 0 otherwise
    
    I = ((b-a)**3/2**M)*sum(fr)
    # Monte-Carlo integral formula
    # range for points sampled raised to the third power because we are integrating over a 3 dimensional region ((b-a)**3)
    # divided by the total number of points sampled (2^M)
    # multipied by the sum of the points sampled that lie within the region (points fr)
    
    return I
    # outputs the value of the integral by Monte-Carlo integration