[Acquisition functions, Part 5] Added the ExpIntVar acquisition. (#239)

* Addresed code-review comments.

* [Acquisition functions, Part 5] Added the ExpIntVar acquisition.

* Moved the calculations of the omega and prior terms outside the evaluate function; Provided the option to control importance sampling.

* Optimised the computation of K.

* Optimised the computation of w.

* Added the option to choose between the NUTS and metropolis samplers for the importance sampling; Altered the term naming and their definitions..

* Altered the term names.

* Added the option for grid integration.

* Addressed the code review comments on the grid integration.

* Altered the descriptions of the integration methods.

CHANGELOG.rst

@@ -9,6 +9,7 @@ dev
 - Added the MaxVar acquisition method
 - Added the RandMaxVar acquisition method
 - Fix elfi.Distance to support scipy 1.0.0
+- Added the ExpIntVar acquisition method
 
 0.6.3 (2017-09-28)
 ------------------

elfi/methods/bo/acquisition.py

@@ -3,6 +3,7 @@
 import logging
 
 import numpy as np
+import scipy.linalg as sl
 import scipy.stats as ss
 
 import elfi.methods.mcmc as mcmc
@@ -550,6 +551,190 @@ def _evaluate_logpdf(theta):
         return batch_theta
 
 
+class ExpIntVar(MaxVar):
+    r"""The Expected Integrated Variance (ExpIntVar) acquisition method.
+
+    Essentially, we define a loss function that measures the overall uncertainty
+    in the unnormalised ABC posterior over the parameter space.
+    The value of the loss function depends on the next simulation and thus
+    the next evaluation location \theta^* is chosen to minimise the expected loss.
+
+    \theta_{t+1} = arg min_{\theta^* \in \Theta} L_{1:t}(\theta^*), where
+
+    \Theta is the parameter space, and L is the expected loss function approximated as follows:
+
+    L_{1:t}(\theta^*) \approx 2 * \sum_{i=1}^s (\omega^i * p^2(\theta^i)
+                                * w_{1:t+1})(theta^i, \theta^*), where
+
+    \omega^i is an importance weight,
+    p^2(\theta^i) is the prior squared, and
+    w_{1:t+1})(theta^i, \theta^*) is the expected variance of the unnormalised ABC posterior
+    at \theta^i after running the simulation model with parameter \theta^*
+
+    References
+    ----------
+    [1] arXiv:1704.00520 (Järvenpää et al., 2017)
+
+    """
+
+    def __init__(self, quantile_eps=.01, integration='grid', d_grid=.2,
+                 n_samples_imp=100, iter_imp=2, sampler='nuts', n_samples=2000,
+                 sigma_proposals_metropolis=None, *args, **opts):
+        """Initialise ExpIntVar.
+
+        Parameters
+        ----------
+        quantile_eps : int, optional
+            Quantile of the observed discrepancies used in setting the discrepancy threshold.
+        integration : str, optional
+            Integration method. Options:
+            - grid (points are taken uniformly): more accurate yet
+            computationally expensive in high dimensions;
+            - importance (points are taken based on the importance weight): less accurate though
+            applicable in high dimensions.
+        d_grid : float, optional
+            Grid tightness.
+        n_samples_imp : int, optional
+            Number of importance samples.
+        iter_imp : int, optional
+            Gap between acquisition iterations in performing importance sampling.
+        sampler : string, optional
+            Sampler for generating random numbers from the proposal distribution for IS.
+            (Options: metropolis, nuts.)
+        n_samples : int, optional
+            Chain length for the sampler that generates the random numbers
+            from the proposal distribution for IS.
+        sigma_proposals_metropolis : array_like, optional
+            Standard deviation proposals for tuning the metropolis sampler.
+
+        """
+        super(ExpIntVar, self).__init__(quantile_eps, *args, **opts)
+        self.name = 'exp_int_var'
+        self.label_fn = 'Expected Loss'
+        self._integration = integration
+        self._n_samples_imp = n_samples_imp
+        self._iter_imp = iter_imp
+
+        if self._integration == 'importance':
+            self.density_is = RandMaxVar(model=self.model,
+                                         prior=self.prior,
+                                         n_inits=self.n_inits,
+                                         seed=self.seed,
+                                         quantile_eps=self.quantile_eps,
+                                         sampler=sampler,
+                                         n_samples=n_samples,
+                                         sigma_proposals_metropolis=sigma_proposals_metropolis)
+        elif self._integration == 'grid':
+            grid_param = [slice(b[0], b[1], d_grid) for b in self.model.bounds]
+            self.points_int = np.mgrid[grid_param].reshape(len(self.model.bounds), -1).T
+
+    def acquire(self, n, t):
+        """Acquire a batch of acquisition points.
+
+        Parameters
+        ----------
+        n : int
+            Number of acquisitions.
+        t : int
+            Current iteration.
+
+        Returns
+        -------
+        array_like
+            Coordinates of the yielded acquisition points.
+
+        """
+        logger.debug('Acquiring the next batch of %d values', n)
+        gp = self.model
+        self.sigma2_n = gp.noise
+
+        # Updating the discrepancy threshold.
+        self.eps = np.percentile(gp.Y, self.quantile_eps * 100)
+
+        # Performing the importance sampling step every self._iter_imp iterations.
+        if self._integration == 'importance' and t % self._iter_imp == 0:
+            self.points_int = self.density_is.acquire(self._n_samples_imp)
+
+        # Obtaining the omegas_int and priors_int terms to be used in the evaluate function.
+        self.mean_int, self.var_int = gp.predict(self.points_int, noiseless=True)
+        self.priors_int = (self.prior.pdf(self.points_int)**2)[np.newaxis, :]
+        if self._integration == 'importance' and t % self._iter_imp == 0:
+            omegas_int_unnormalised = (1 / MaxVar.evaluate(self, self.points_int)).T
+            self.omegas_int = omegas_int_unnormalised / \
+                np.sum(omegas_int_unnormalised, axis=1)[:, np.newaxis]
+        elif self._integration == 'grid':
+            self.omegas_int = np.empty(len(self.points_int))
+            self.omegas_int.fill(1 / len(self.points_int))
+
+        # Initialising the attributes used in the evaluate function.
+        self.thetas_old = np.array(gp.X)
+        self._K = gp._gp.kern.K
+        self.K = self._K(self.thetas_old, self.thetas_old) + \
+            self.sigma2_n * np.identity(self.thetas_old.shape[0])
+        self.k_int_old = self._K(self.points_int, self.thetas_old).T
+        self.phi_int = ss.norm.cdf(self.eps, loc=self.mean_int.T,
+                                   scale=np.sqrt(self.sigma2_n + self.var_int.T))
+
+        # Obtaining the location where the expected loss is minimised.
+        # Note: The gradient is computed numerically as GPy currently does not
+        # directly provide the derivative computations used in Järvenpää et al., 2017.
+        theta_min, _ = minimize(self.evaluate,
+                                gp.bounds,
+                                grad=None,
+                                prior=self.prior,
+                                n_start_points=self.n_inits,
+                                maxiter=self.max_opt_iters,
+                                random_state=self.random_state)
+
+        # Using the same location for all points in the batch.
+        batch_theta = np.tile(theta_min, (n, 1))
+        return batch_theta
+
+    def evaluate(self, theta_new, t=None):
+        """Evaluate the acquisition function at the location theta_new.
+
+        Parameters
+        ----------
+        theta_new : array_like
+            Evaluation coordinates.
+        t : int, optional
+            Current iteration, (unused).
+
+        Returns
+        -------
+        array_like
+            Expected loss's term dependent on theta_new.
+
+        """
+        gp = self.model
+        n_imp, n_dim = self.points_int.shape
+        # Alter the shape of theta_new.
+        if n_dim != 1 and theta_new.ndim == 1:
+            theta_new = theta_new[np.newaxis, :]
+        elif n_dim == 1 and theta_new.ndim == 1:
+            theta_new = theta_new[:, np.newaxis]
+
+        # Calculate the integrand term w.
+        # Note: w's second term (given in Järvenpää et al., 2017) is dismissed
+        # because it is constant with respect to theta_new.
+        _, var_new = gp.predict(theta_new, noiseless=True)
+        k_old_new = self._K(self.thetas_old, theta_new)
+        k_int_new = self._K(self.points_int, theta_new).T
+        # Using the Cholesky factorisation to avoid computing matrix inverse.
+        term_chol = sl.cho_solve(sl.cho_factor(self.K), k_old_new)
+        cov_int = k_int_new - np.dot(self.k_int_old.T, term_chol).T
+        delta_var_int = cov_int**2 / (self.sigma2_n + var_new)
+        a = np.sqrt((self.sigma2_n + self.var_int.T - delta_var_int) /
+                    (self.sigma2_n + self.var_int.T + delta_var_int))
+        # Using the skewnorm's cdf to substitute the Owen's T function.
+        phi_skew_imp = ss.skewnorm.cdf(self.eps, a, loc=self.mean_int.T,
+                                       scale=np.sqrt(self.sigma2_n + self.var_int.T))
+        w = ((self.phi_int - phi_skew_imp) / 2)
+
+        loss_theta_new = 2 * np.sum(self.omegas_int * self.priors_int * w, axis=1)
+        return loss_theta_new
+
+
 class UniformAcquisition(AcquisitionBase):
     """Acquisition from uniform distribution."""
 

tests/conftest.py

@@ -12,7 +12,7 @@
 import elfi.examples.gauss
 import elfi.examples.ma2
 from elfi.methods.bo.gpy_regression import GPyRegression
-from elfi.methods.bo.acquisition import MaxVar, RandMaxVar
+from elfi.methods.bo.acquisition import ExpIntVar, MaxVar, RandMaxVar
 from elfi.methods.utils import ModelPrior
 
 elfi.clients.native.set_as_default()
@@ -129,6 +129,23 @@ def acq_randmaxvar():
     return method_acq
 
 
+@pytest.fixture()
+def acq_expintvar():
+    """Initialise an ExpIntVar fixture.
+
+    Returns
+    -------
+    ExpIntVar
+        Acquisition method.
+
+    """
+    gp, prior = _get_dependencies_acq_fn()
+
+    # Initialising the acquisition method.
+    method_acq = ExpIntVar(model=gp, prior=prior)
+    return method_acq
+
+
 def _get_dependencies_acq_fn():
     """Provide the requirements for the MaxVar-based acquisition function initialisation.
 

tests/unit/test_bo.py

@@ -197,3 +197,30 @@ def test_acq_bounds(self, acq_randmaxvar):
         for dim in range(n_dim_fixture):
             assert np.all((batch_theta[:, dim] >= bounds[dim][0]) &
                           (batch_theta[:, dim] <= bounds[dim][1]))
+
+
+class Test_ExpIntVar:
+    """Run a collection of tests for the ExpIntVar acquisition."""
+
+    def test_acq_bounds(self, acq_expintvar):
+        """Check if the acquisition is performed within the bounds.
+
+        Parameters
+        ----------
+        acq_expintvar : ExpIntVar
+            Acquisition method.
+
+        """
+        bounds = acq_expintvar.model.bounds
+        n_dim_fixture = len(acq_expintvar.model.bounds)
+        batch_size = 2
+        n_it = 2
+
+        # Acquiring points.
+        for it in range(n_it):
+            batch_theta = acq_expintvar.acquire(n=batch_size, t=it)
+
+        # Checking if the acquired points are within the bounds.
+        for dim in range(n_dim_fixture):
+            assert np.all((batch_theta[:, dim] >= bounds[dim][0]) &
+                          (batch_theta[:, dim] <= bounds[dim][1]))