init and main methods of the MISER class

miser_1.py

from numpy.random import uniform

class MISER:
    def __init__(self, miser_params):
        """
        Parameters
        ----------
        MNBS : int
            If less than MNBS evaluations are left we do
            vanilla MC.
        MNPTS : int
            Minimum number of points used for variance
            exploration.
        pfac : float in [0,1]
            Fraction of the remaining points used in each
            recursive call. Typical value is 0.1
        alpha : float
            Factor used to estimate the optimal points to be
            used on each side of a bisection. Usually set to 2
        """
        self.MNBS = miser_params["MNBS"]
        self.MNPTS = miser_params["MNPTS"]
        self.pfac = miser_params["pfac"]
        self.alpha = miser_params["alpha"]
        
        # We store the points of all the evaluations used
        # in the vanilla MC step so that we can plot them
        # later.
        self.POINTS = []
        
    def integrate(self, f, N_points, xl, xu):
        self.POINTS = []
        integral_estimate, variance_estimate = self.miser_kernel(f, N_points, xl, xu)
        
        return (integral_estimate, variance_estimate)
        
    def miser_kernel(self, f, N_points, xl, xu):
        """
        MISER algorithm.
        
        Parameters
        ----------
        f : function
            The integrand, it needs to be vectorized.
        N_points : int
            Total number of function evaluations available.
        xl : [float]
            List of lower bounds of the integral.
        xu : [float]
            List of upper bounds of the integral.
        
        Returns
        -------
        ave : float
            Estimate of the integral.
        var : float
            Estimated variance of the integral.
        """
        ndims = len(xu)

        if N_points < self.MNBS: # Straight Monte Carlo
            ave, var = self._mc_kernel(f, xu, xl, N_points)
            
            return ave, var
            
        else: # Recursion path
            # 1) Calculate number of points for this
            #    variance exploration step
            npre = int(max(self.pfac * N_points, self.MNPTS))
            
            # 2) Find best split         
            newxu, newxl, var_l, var_u = self._split_domain(f, xl, xu, npre)

            # 3) Compute new allocation of points for split regions
            N_points_remaining = N_points - npre
            N_points_l, N_points_u = self._allocate_points(var_l, var_u, N_points_remaining)

            # 4) Recurse
            ave_l, var_l = self.miser_kernel(f, N_points_l, xl, newxu)
            ave_u, var_u = self.miser_kernel(f, N_points_u, newxl, xu)
            
            # 5) Combine results
            ave = 0.5 * (ave_l + ave_u)
            var = 0.25 * (var_l + var_u)

            return (ave, var)