Merge pull request #299 from QMCSoftware/mc_vec

Mc vec

demos/elliptic-pde.ipynb

@@ -486,7 +486,7 @@
     "        self.compute_eigenpairs()\n",
     "        self.sampler = sampler\n",
     "        self.true_measure = qp.Gaussian(self.sampler)\n",
-    "        self.dprime = 1\n",
+    "        self.rho = 1\n",
     "        super(EllipticPDE, self).__init__()\n",
     "        \n",
     "    def compute_eigenpairs(self):\n",
@@ -718,7 +718,7 @@
     "        self.leveltype = \"adaptive-multi\"\n",
     "        self.sums = np.zeros(6)\n",
     "        self.cost = 0\n",
-    "        self.dprime = 1\n",
+    "        self.rho = 1\n",
     "        super(MLEllipticPDE, self).__init__()\n",
     "    \n",
     "    def _spawn(self, level, sampler):\n",

environment.yml

@@ -8,8 +8,8 @@ dependencies:
   - sphinx=5.0.2
   - pandoc>=2.2
   - pip:
-    - numpy==1.23.0
-    - scipy==1.8.1
+    - numpy==1.23.3
+    - scipy==1.9.2
     - pandas==1.4.3
     - matplotlib==3.5.2
     - jupyter==1.0.0

qmcpy/accumulate_data/ld_transform_bayes_data.py

@@ -18,7 +18,7 @@ class LDTransformBayesData(AccumulateData):
     """
 
     def __init__(self, stopping_crit, integrand, true_measure, discrete_distrib, m_min: int, m_max: int,
-                 fbt, merge_fbt, kernel):
+                 fbt, merge_fbt, kernel, alpha):
         """
         Args:
             stopping_crit (StoppingCriterion): a StoppingCriterion instance
@@ -49,9 +49,9 @@ def __init__(self, stopping_crit, integrand, true_measure, discrete_distrib, m_m
         # quantile value for the error bound
         if self.errbd_type == 'full_Bayes':
             # degrees of freedom = 2^mmin - 1
-            self.uncert = -tnorm.ppf(self.stopping_crit.alpha / 2, (2 ** m_min) - 1)
+            self.uncert = -tnorm.ppf(alpha / 2, (2 ** m_min) - 1)
         else:
-            self.uncert = -gaussnorm.ppf(self.stopping_crit.alpha / 2)
+            self.uncert = -gaussnorm.ppf(alpha / 2)
 
         # Set Attributes
         self.m_min = m_min

qmcpy/discrete_distribution/digital_net_b2/sobol_scipy.py

@@ -1,6 +1,8 @@
-from numpy.core.numeric import isscalar
 from .._discrete_distribution import DiscreteDistribution
 from ...util import ParameterError, ParameterWarning
+from numpy.core.numeric import isscalar
+from numpy import *
+import warnings
 
 class SobolSciPy(DiscreteDistribution):
     """
@@ -9,9 +11,9 @@ class SobolSciPy(DiscreteDistribution):
     >>> s = SobolSciPy(2,seed=7)
     >>> s.gen_samples(4)
     array([[5.79259991e-01, 7.40284680e-01],
-           [4.15829644e-02, 6.92069530e-04],
-           [4.78844851e-01, 7.75258362e-01],
-           [8.92499685e-01, 4.83783960e-01]])
+           [4.15829662e-02, 6.92069530e-04],
+           [4.78844853e-01, 7.75258361e-01],
+           [8.92499692e-01, 4.83783960e-01]])
     >>> s.gen_samples(1)
     array([[0.57925999, 0.74028468]])
     >>> s.gen_samples(1)

qmcpy/integrand/__init__.py

@@ -6,6 +6,6 @@
 from .custom_fun import CustomFun
 from .ml_call_options import MLCallOptions
 from .box_integral import BoxIntegral
-from .sobol_indices import SobolIndices,SensitivityIndices
+from .sensitivity_indices import SobolIndices,SensitivityIndices
 from .ishigami import Ishigami
 from .bayesian_lr_coeffs import BayesianLRCoeffs

qmcpy/integrand/_integrand.py

@@ -10,17 +10,19 @@
 class Integrand(object):
     """ Integrand abstract class. DO NOT INSTANTIATE. """
 
-    def __init__(self, dprime, parallel):
+    def __init__(self, rho, eta, parallel):
         """
         Args:
-            dprime (tuple): function output dimension shape.
+            rho (tuple): individual solution shape.
+            eta (tuple): combined solution shape. 
             parallel (int): If parallel is False, 0, or 1: function evaluation is done in serial fashion.
                 Otherwise, parallel specifies the number of CPUs used by multiprocessing.Pool.
                 Passing parallel=True sets the number of CPUs equal to os.cpu_count().
         """
         prefix = 'A concrete implementation of Integrand must have '
         self.d = self.true_measure.d
-        self.dprime = (dprime,) if isinstance(dprime,int) else tuple(dprime)
+        self.rho = (rho,) if isinstance(rho,int) else tuple(rho)
+        self.eta = (eta,) if isinstance(eta,int) else tuple(eta)
         cpus = os.cpu_count()
         self.parallel = cpus if parallel is True else int(parallel)
         self.parallel = 0 if self.parallel==1 else self.parallel
@@ -77,7 +79,7 @@ def f(self, x, periodization_transform='NONE', compute_flags=None, *args, **kwar
             ndarray: length n vector of function evaluations
         """
         periodization_transform = 'NONE' if periodization_transform is None else periodization_transform.upper()
-        compute_flags = tile(1,self.dprime) if compute_flags is None else atleast_1d(compute_flags)
+        compute_flags = tile(1,self.rho) if compute_flags is None else atleast_1d(compute_flags)
         n,d = x.shape
         # parameter checks
         if self.discrete_distrib.mimics != 'StdUniform' and periodization_transform!='NONE':
@@ -113,7 +115,7 @@ def f(self, x, periodization_transform='NONE', compute_flags=None, *args, **kwar
             xp[xp<=0] = self.EPS
             xp[xp>=1] = 1-self.EPS
         # function evaluation with chain rule
-        y = empty((n,)+self.dprime,dtype=float)
+        y = empty((n,)+self.rho,dtype=float)
         if self.true_measure == self.true_measure.transform:
             # jacobian*weight/pdf will cancel so f(x) = g(\Psi(x))
             xtf = self.true_measure._transform(xp) # get transformed samples, equivalent to self.true_measure._transform_r(x)
@@ -138,14 +140,14 @@ def _g(self,t,compute_flags,*args,**kwargs):
             y = concatenate(y,dtype=float)
         else:
             y = self._g2(t,comb_args=(args,kwargs))
-        y = y.reshape((n,)+self.dprime)
+        y = y.reshape((n,)+self.rho)
         return y
 
     def _g2(self,t,comb_args=((),{})):
         args = comb_args[0]
         kwargs = comb_args[1]
         t = atleast_2d(t)
-        if self.dprime==(1,):
+        if self.rho==(1,):
             kwargs = dict(kwargs)
             del kwargs['compute_flags']
         y = self.g(t,*args,**kwargs)
@@ -159,8 +161,8 @@ def bound_fun(self, bound_low, bound_high):
         individually.
 
         Args:
-            bound_low (ndarray): length Integrand.dprime lower error bound
-            bound_high (ndarray): length Integrand.dprime upper error bound
+            bound_low (ndarray): length Integrand.rho lower error bound
+            bound_high (ndarray): length Integrand.rho upper error bound
 
         Return:
             (tuple) containing
@@ -169,21 +171,28 @@ def bound_fun(self, bound_low, bound_high):
             - (ndarray): upper bound on function combining estimates
             - (ndarray): bool flags to override sufficient combined integrand estimation, e.g., when approximating a ratio of integrals, if the denominator's bounds straddle 0, then returning True here forces ratio to be flagged as insufficiently approximated.
         """
+        if self.rho!=self.eta:
+            raise ParameterError('''
+                Set bound_fun explicity. 
+                The default bound_fun is the identity map. 
+                Since individual solution dimension rho != combined solution dimension eta, 
+                QMCPy cannot infer a reasonable bound function.''')
         return bound_low,bound_high
 
-    def dependency(self, flags_comb):
+    def dependency(self, comb_flags):
         """
-        takes a vector of indicators of weather of not
+        Takes a vector of indicators of weather of not
         the error bound is satisfied for combined integrands and which returns flags for individual integrands.
         For example, if we are taking the ratio of 2 individual integrands, then getting flag_comb=True means the ratio
         has not been approximated to within the tolerance, so the dependency function should return [True,True]
         indicating that both the numerator and denominator integrands need to be better approximated.
         Args:
-            flags_comb (bool ndarray): flags indicating weather the combined integrals are insufficiently approximated
+            comb_flags (bool ndarray): flags indicating weather the combined integrals are insufficiently approximated
 
         Return:
-            (bool ndarray): length (Integrand.dprime) flags for individual integrands"""
-        return flags_comb
+            (bool ndarray): length (Integrand.rho) flags for individual integrands
+        """
+        return comb_flags if self.rho==self.eta else  tile((comb_flags==False).any(),self.rho)
 
     def _dimension_at_level(self, level):
         """

qmcpy/integrand/asian_option.py

@@ -83,7 +83,7 @@ def __init__(self, sampler, volatility=0.5, start_price=30., strike_price=35.,\
             self.leveltype = 'single'
             self.parent = False
             self.parameters += ['dim_frac']
-        super(AsianOption,self).__init__(dprime=1,parallel=False)    
+        super(AsianOption,self).__init__(rho=1,eta=1,parallel=False)    
 
     def _get_discounted_payoffs(self, stock_path, dimension):
         """

qmcpy/integrand/bayesian_lr_coeffs.py

@@ -14,10 +14,10 @@ class BayesianLRCoeffs(Integrand):
     >>> x = blrcoeffs.discrete_distrib.gen_samples(2**10)
     >>> y = blrcoeffs.f(x)
     >>> y.shape
-    (1024, 6)
+    (1024, 2, 3)
     >>> y.mean(0)
-    array([ 0.04639394, -0.01440543, -0.05498496,  0.02176581,  0.02176581,
-            0.02176581])
+    array([[ 0.04639394, -0.01440543, -0.05498496],
+           [ 0.02176581,  0.02176581,  0.02176581]])
     """
 
     def __init__(self, sampler, feature_array, response_vector, prior_mean=0, prior_covariance=10):
@@ -46,17 +46,16 @@ def __init__(self, sampler, feature_array, response_vector, prior_mean=0, prior_
         if self.response_vector.shape!=(obs,) or ((self.response_vector!=0)&(self.response_vector!=1)).any():
             ParameterError("response_vector must have the same length as feature_array and contain only 0 or 1 enteries.")
         self.feature_array = column_stack((self.feature_array,ones((obs,1))))
-        self.dprime = 2*self.num_coeffs
-        super(BayesianLRCoeffs,self).__init__(dprime=self.dprime,parallel=False)
+        super(BayesianLRCoeffs,self).__init__(rho=(2,self.num_coeffs),eta=self.num_coeffs,parallel=False)
 
     def g(self, x, compute_flags):
         z = x@self.feature_array.T
         z1 = z*self.response_vector
         with errstate(over='ignore'):
             den = exp(sum(z1-log(1+exp(z)),1))[:,None]
-        y = zeros((len(x),2*self.num_coeffs),dtype=float)
-        y[:,:self.num_coeffs] = x*den
-        y[:,self.num_coeffs:] = den
+        y = zeros((len(x),2,self.num_coeffs),dtype=float)
+        y[:,0] = x*den
+        y[:,1] = den
         return y
 
     def _spawn(self, level, sampler):
@@ -68,13 +67,13 @@ def _spawn(self, level, sampler):
             prior_covariance = self.prior_covariance)
 
     def bound_fun(self, bound_low, bound_high):
-        num_bounds_low,den_bounds_low = bound_low[:self.num_coeffs],bound_low[self.num_coeffs:]
-        num_bounds_high,den_bounds_high = bound_high[:self.num_coeffs],bound_high[self.num_coeffs:]
+        num_bounds_low,den_bounds_low = bound_low[0],bound_low[1]
+        num_bounds_high,den_bounds_high = bound_high[0],bound_high[1]
         comb_bounds_low = minimum.reduce([num_bounds_low/den_bounds_low,num_bounds_high/den_bounds_low,num_bounds_low/den_bounds_high,num_bounds_high/den_bounds_high])
         comb_bounds_high = maximum.reduce([num_bounds_low/den_bounds_low,num_bounds_high/den_bounds_low,num_bounds_low/den_bounds_high,num_bounds_high/den_bounds_high])
         violated = (den_bounds_low<=0)*(0<=den_bounds_high)
         comb_bounds_low[violated],comb_bounds_high[violated] = -inf,inf
         return comb_bounds_low,comb_bounds_high
 
-    def dependency(self, flags_comb):
-        return hstack((flags_comb,flags_comb))
+    def dependency(self, comb_flags):
+        return vstack((comb_flags,comb_flags))

qmcpy/integrand/box_integral.py

@@ -44,12 +44,12 @@ def __init__(self, sampler, s=array([1,2])):
         self.s = array([s]) if isscalar(s) else array(s)
         self.sampler = sampler
         self.true_measure = Uniform(self.sampler, lower_bound=0., upper_bound=1.)
-        super(BoxIntegral,self).__init__(dprime=len(self.s),parallel=False) # output dimensions per sample
+        super(BoxIntegral,self).__init__(rho=len(self.s),eta=len(self.s),parallel=False) # output dimensions per sample
 
     def g(self, t, **kwargs):
-        compute_flags = kwargs['compute_flags'] if 'compute_flags' in kwargs else ones(self.dprime,dtype=int)
+        compute_flags = kwargs['compute_flags'] if 'compute_flags' in kwargs else ones(self.rho,dtype=int)
         n,d = t.shape
-        y = zeros((n,)+self.dprime,dtype=float)
+        y = zeros((n,)+self.rho,dtype=float)
         for j in range(len(self.s)):
             if compute_flags[j] == 1: 
                 y[:,j] = (t**2).sum(1)**(self.s[j]/2)

qmcpy/integrand/custom_fun.py

@@ -10,7 +10,7 @@ class CustomFun(Integrand):
     >>> cf = CustomFun(
     ...     true_measure = Gaussian(DigitalNetB2(2,seed=7),mean=[1,2]),
     ...     g = lambda x: x[:,0]**2*x[:,1],
-    ...     dprime = 1)
+    ...     rho = 1)
     >>> x = cf.discrete_distrib.gen_samples(2**10)
     >>> y = cf.f(x)
     >>> y.shape
@@ -20,7 +20,7 @@ class CustomFun(Integrand):
     >>> cf = CustomFun(
     ...     true_measure = Uniform(DigitalNetB2(3,seed=7),lower_bound=[2,3,4],upper_bound=[4,5,6]),
     ...     g = lambda x,compute_flags=None: x,
-    ...     dprime = 3)
+    ...     rho = 3)
     >>> x = cf.discrete_distrib.gen_samples(2**10)
     >>> y = cf.f(x)
     >>> y.shape
@@ -29,12 +29,12 @@ class CustomFun(Integrand):
     array([3., 4., 5.])
     """
 
-    def __init__(self, true_measure, g, dprime=(1,), parallel=False):
+    def __init__(self, true_measure, g, rho=1, parallel=False):
         """
         Args:
             true_measure (TrueMeasure): a TrueMeasure instance. 
             g (function): a function handle. 
-            dprime (tuple): function output dimension shape.
+            rho (tuple): individual solution dimensions.
             parallel (int): If parallel is False, 0, or 1: function evaluation is done in serial fashion. 
                 Otherwise, parallel specifies the number of CPUs used by multiprocessing.Pool. 
                 Passing parallel=True sets the number of CPUs equal to os.cpu_count(). 
@@ -43,8 +43,8 @@ def __init__(self, true_measure, g, dprime=(1,), parallel=False):
         self.parameters = []
         self.true_measure = true_measure
         self.sampler = self.true_measure
-        self.__g = g
-        super(CustomFun,self).__init__(dprime,parallel)
+        self.__g = g            
+        super(CustomFun,self).__init__(rho=rho,eta=rho,parallel=parallel)
 
     def g(self, t, *args, **kwargs):
         return self.__g(t,*args,**kwargs)
@@ -53,5 +53,5 @@ def _spawn(self, level, sampler):
         return CustomFun(
             true_measure = sampler,
             g = self.__g,
-            dprime = self.dprime,
+            rho = self.rho,
             parallel = self.parallel)

qmcpy/integrand/european_option.py

@@ -56,7 +56,7 @@ def __init__(self, sampler, volatility=0.5, start_price=30, strike_price=35,
         self.call_put = call_put.lower()
         if self.call_put not in ['call','put']:
             raise ParameterError("call_put must be either 'call' or 'put'")
-        super(EuropeanOption,self).__init__(dprime=1,parallel=False)  
+        super(EuropeanOption,self).__init__(rho=1,eta=1,parallel=False)  
 
     def g(self, t):
         """ See abstract method. """

qmcpy/integrand/ishigami.py

@@ -41,7 +41,7 @@ def __init__(self,sampler, a=7, b=.1):
         self.a = a
         self.b = b
         self.true_measure = Uniform(self.sampler, lower_bound=-pi, upper_bound=pi)
-        super(Ishigami,self).__init__(dprime=1,parallel=False)
+        super(Ishigami,self).__init__(rho=1,eta=1,parallel=False)
 
     def g(self, t):
         y = (1+self.b*t[:,2]**4)*sin(t[:,0])+self.a*sin(t[:,1])**2

qmcpy/integrand/keister.py

@@ -46,7 +46,7 @@ def __init__(self, sampler):
         """
         self.sampler = sampler
         self.true_measure = Gaussian(self.sampler,mean=0,covariance=1/2)
-        super(Keister,self).__init__(dprime=1,parallel=False)
+        super(Keister,self).__init__(rho=1,eta=1,parallel=False)
 
     def g(self, t):
         d = t.shape[1]

qmcpy/integrand/linear0.py

@@ -24,7 +24,7 @@ def __init__(self, sampler):
         """
         self.sampler = sampler
         self.true_measure = Uniform(self.sampler, lower_bound=-.5, upper_bound=.5)
-        super(Linear0,self).__init__(dprime=1,parallel=False)
+        super(Linear0,self).__init__(rho=1,eta=1,parallel=False)
 
     def g(self, t):
         y = t.sum(1)

qmcpy/integrand/ml_call_options.py

@@ -70,7 +70,7 @@ def __init__(self, sampler, option='european', volatility=.2,
         self.g_submodule = getattr(self,'_g_'+self.option)
         self.level = _level
         self.max_level = inf
-        super(MLCallOptions,self).__init__(dprime=1,parallel=False)
+        super(MLCallOptions,self).__init__(rho=1,eta=1,parallel=False)
         #if self.discrete_distrib.low_discrepancy and self.option=='asian':
         #    raise ParameterError('MLCallOptions does not support LD sequence for Asian Option')