Feature/refactor for reduced complexity

.gitignore

@@ -1,2 +1,3 @@
 *.pyc
 __pycache__
+.idea/

README.md

@@ -1,6 +1,14 @@
-# `excursion` — Efficient Excursion Set Estimation 
-
+[![Build Status](https://travis-ci.org/irinaespejo/excursion.svg?branch=master)](https://travis-ci.org/irinaespejo/excursion)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.1634427.svg)](https://zenodo.org/badge/latestdoi/146087019)
+[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/irinaespejo/excursion/master?urlpath=https%3A%2F%2Fgithub.com%2Firinaespejo%2Fexcursion%2Fblob%2Fmaster%2Fexamples%2Ftutorial_1D.ipynb)
+[![Documentation Status](https://readthedocs.org/projects/excursion/badge/?version=latest)](https://excursion.readthedocs.io/en/latest/?badge=latest)
+[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
+[![License](https://img.shields.io/badge/License-Apache%202.0-blue.svg)](https://opensource.org/licenses/Apache-2.0)
+
+
+
+# `excursion` — Efficient Excursion Set Estimation 
+[ReadTheDocs](https://excursion.readthedocs.io/en/latest/)
 
 This package implements a Bayesian Optimization procedure based on Gaussian Processes to efficiently determine excursion sets (or equivalently iso-surfaces) of one or many expensive black-box functions.
 

conf.py

@@ -0,0 +1,52 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# This file only contains a selection of the most common options. For a full
+# list see the documentation:
+# https://www.sphinx-doc.org/en/master/usage/configuration.html
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+import os
+import sys
+
+# sys.path.insert(0, os.path.abspath('.'))
+
+
+# -- Project information -----------------------------------------------------
+
+project = "excursion"
+copyright = "2020, Lukas Heinrich, Irina Espejo, Giles Louppe, Kyle Cranmer"
+author = "Lukas Heinrich, Irina Espejo, Giles Louppe, Kyle Cranmer"
+
+
+# -- General configuration ---------------------------------------------------
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = []
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ["_templates"]
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ["_build", "Thumbs.db", ".DS_Store"]
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = "alabaster"
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ["_static"]

excursion/__init__.py

@@ -1,29 +1,5 @@
-import numpy as np
-from scipy.linalg import cho_solve
-from scipy.stats import norm
-from sklearn.gaussian_process import GaussianProcessRegressor
-from sklearn.gaussian_process.kernels import RBF
-from sklearn.gaussian_process.kernels import ConstantKernel
-from sklearn.gaussian_process.kernels import WhiteKernel
-from sklearn.gaussian_process.kernels import Matern
-
-
-def get_gp(X, y, alpha=10**-7, kernel_name='const_rbf'):
-    if kernel_name == 'const_rbf':
-        length_scale = [1.]*X.shape[-1]
-        kernel = ConstantKernel() * RBF(length_scale_bounds=[0.1, 100.0], length_scale = length_scale)
-    elif kernel_name == 'tworbf_white':
-        kernel = ConstantKernel() * RBF(length_scale_bounds=[1e-2,100]) + \
-                 ConstantKernel() * RBF(length_scale_bounds=[100., 1000.0]) + \
-                 WhiteKernel(noise_level_bounds=[1e-7,1e-4])
-    elif kernel_name == 'onerbf_white':
-        kernel = ConstantKernel() * RBF(length_scale_bounds=[1e-2,100]) + WhiteKernel(noise_level_bounds=[1e-7,1e-1])
-    else:
-        raise RuntimeError('unknown kernel')
-    gp = GaussianProcessRegressor(kernel=kernel,
-                                  n_restarts_optimizer=10,
-                                  alpha=alpha,
-                                  random_state=1234)
-    gp.fit(X, y.ravel())
-    return gp
-
+from .plotting import plot
+from .excursion import ExcursionProblem
+from .optimizer import Optimizer
+from .learner import Learner
+__all__ = ["plot", "ExcursionProblem", "Optimizer", "Learner"]

excursion/acquisition/MES.py

@@ -0,0 +1,62 @@
+import torch
+from .utils import cdf
+from .base import AcquisitionFunction
+
+
+class MES(AcquisitionFunction):
+    def __init__(self, device=None, dtype=None, batch=False, ):
+        self._prev_acq_point_index = []
+        self.device = device
+        self.dtype = dtype
+        self.batch = batch
+        self.acq_vals = None
+
+    def acquire(self, gp, thresholds, X_pointsgrid):
+    # compute predictive posterior of Y(x) | trin data
+        likelihood = gp.likelihood
+        gp.eval()
+        likelihood.eval()
+
+    # ok
+        prediction = likelihood(gp(X_pointsgrid))
+        pred_mean = prediction.mean
+
+        pred_stdev = torch.sqrt(prediction.variance)
+
+        num_points = X_pointsgrid.size()[0]
+        entropy_grid = torch.zeros(num_points,).to(device=self.device, dtype=self.dtype)
+
+        for j in range(len(thresholds) - 1):
+            p_j = cdf(pred_mean, pred_stdev, thresholds[j + 1]) - cdf(pred_mean, pred_stdev, thresholds[j])
+
+            entropy_grid[p_j > 0] -= torch.log(torch.exp(p_j[p_j > 0])) * torch.log(torch.exp(torch.log(p_j[p_j > 0])))
+
+        acq_cand_vals = entropy_grid
+        self.acq_vals = torch.clone(entropy_grid)
+
+        if self.batch:
+            print("Will implement this later for batched results\n")
+            pass
+
+        return X_pointsgrid[self.get_first_max_index(gp, X_pointsgrid, acq_cand_vals)]
+
+    def get_first_max_index(self, gp, X_pointsgrid, acq_cand_vals):
+        X_train = gp.train_inputs[0].to(device=self.device, dtype=self.dtype)
+        X_train = X_train.tolist()
+        new_index = torch.argmax(acq_cand_vals)
+
+        # if the index is not already picked nor in the training set
+        # accept it ans remove from future picks
+        return self._check_prev_acq(new_index, X_train, X_pointsgrid, acq_cand_vals)
+
+    # recursion helper function
+    def _check_prev_acq(self, new_index, X_train, X_pointsgrid, acq_cand_vals):
+        new_X = X_pointsgrid[new_index]
+        if (new_index not in self._prev_acq_point_index) and (
+                new_X.tolist() not in X_train):
+            self._prev_acq_point_index.append(new_index.item())
+            return new_index.item()
+        else:
+            acq_cand_vals[new_index] = torch.Tensor([(-1.0) * float("Inf")])
+            new_index = torch.argmax(acq_cand_vals)
+            return self._check_prev_acq(new_index, X_train, X_pointsgrid, acq_cand_vals)

excursion/acquisition/PES.py

@@ -0,0 +1,103 @@
+import torch
+from .utils import truncated_std_conditional
+from .utils import h_normal
+from .base import AcquisitionFunction
+
+
+class PES(AcquisitionFunction):
+    def __init__(self, device=None, dtype=None, batch=False, ):
+        self._prev_acq_point_index = []
+        self.device = device
+        self.dtype = dtype
+        self.batch = batch
+        self.acq_vals = None
+
+    def _acquire(self, gp, thresholds, X_pointsgrid, x_candidate):
+
+        # compute predictive posterior of Y(x) | train data
+        kernel = gp.covar_module
+        likelihood = gp.likelihood
+        gp.eval()
+        likelihood.eval()
+        X_all = torch.cat((x_candidate, X_pointsgrid))
+        Y_pred_all = likelihood(gp(X_all))
+        Y_pred_grid = torch.distributions.Normal(
+            loc=Y_pred_all.mean[1:], scale=(Y_pred_all.variance[1:]) ** 0.5)
+
+        # vector of expected value H1 under S(x) for each x in X_grid
+        # E_S_H1 = torch.zeros(len(X_grid)).to(device=self.device, dtype=self.dtype)
+        E_S_H1 = torch.zeros(len(X_pointsgrid)).to(device=self.device, dtype=self.dtype)
+
+        for j in range(len(thresholds) - 1):
+            # vector of sigma(Y(x_candidate)|S(x)=j) truncated
+            trunc_std_j = truncated_std_conditional(
+                Y_pred_all, thresholds[j], thresholds[j + 1]
+            )
+            H1_j = h_normal(trunc_std_j)
+
+            # vector of p(S(x)=j)
+            p_j = Y_pred_grid.cdf(thresholds[j + 1]) - Y_pred_grid.cdf(thresholds[j])
+            mask = torch.where(p_j == 0.0)
+            H1_j[mask] = 0.0
+            # print(p_j)
+            # print(H1_j)
+            E_S_H1 += p_j * H1_j  # expected value of H1 under S(x)
+
+        # entropy of Y(x_candidate)
+        H0 = h_normal(Y_pred_all.variance[0] ** 0.5)
+
+        # info gain on the grid, vector
+        info_gain = H0 - E_S_H1
+
+        info_gain[~torch.isfinite(info_gain)] = 0.0  # just to avoid NaN
+
+        # cumulative info gain over grid
+        cumulative_info_gain = info_gain.sum()
+        return cumulative_info_gain.item()
+
+    def acquire(self, gp, thresholds, X_pointsgrid):
+        """
+        Calculates information gain of choosing x_candidate as next point to evaluate.
+        Performs this calculation with the Predictive Entropy Search approximation. Roughly,
+        PES(x_candidate) = int dx { H[Y(x_candidate)] - E_{S(x=j)} H[Y(x_candidate)|S(x=j)] }
+        Notation: PES(x_candidate) = int dx H0 - E_Sj H1
+
+        """
+
+        # compute predictive posterior of Y(x) | train data
+        # likelihood = gp.likelihood
+        # gp.eval()
+        # likelihood.eval()
+        acquisition_values = []
+
+        # for x_candidate in X_grid:
+        for x_candidate in X_pointsgrid:
+            x_candidate = x_candidate.view(1, -1).to(device=self.device, dtype=self.dtype)
+            acquisition_values.append(self._acquire(gp, thresholds, X_pointsgrid, x_candidate))
+
+        acq_cand_vals = torch.Tensor(acquisition_values).to(device=self.device, dtype=self.dtype)
+        self.acq_vals = torch.clone(acq_cand_vals)
+
+        return X_pointsgrid[self.get_first_max_index(gp, X_pointsgrid, acq_cand_vals)]
+
+
+    def get_first_max_index(self, gp, X_pointsgrid, acq_cand_vals):
+        X_train = gp.train_inputs[0].to(device=self.device, dtype=self.dtype)
+        X_train = X_train.tolist()
+        new_index = torch.argmax(acq_cand_vals)
+
+        # if the index is not already picked nor in the training set
+        # accept it ans remove from future picks
+        return self._check_prev_acq(new_index, X_train, X_pointsgrid, acq_cand_vals)
+
+    # recursion helper function
+    def _check_prev_acq(self, new_index, X_train, X_pointsgrid, acq_cand_vals):
+        new_X = X_pointsgrid[new_index]
+        if (new_index not in self._prev_acq_point_index) and (
+                new_X.tolist() not in X_train):
+            self._prev_acq_point_index.append(new_index.item())
+            return new_index.item()
+        else:
+            acq_cand_vals[new_index] = torch.Tensor([(-1.0) * float("Inf")])
+            new_index = torch.argmax(acq_cand_vals)
+            return self._check_prev_acq(new_index, X_train, X_pointsgrid, acq_cand_vals)

excursion/acquisition/__init__.py

@@ -0,0 +1,9 @@
+from .MES import MES
+from .PES import PES
+from .base import AcquisitionFunction
+
+__all__ = [
+    "PES",
+    "MES",
+    "AcquisitionFunction",
+]

excursion/acquisition/approximations.py

@@ -0,0 +1,75 @@
+# from scipy.stats import norm
+# import torch
+# import numpy as np
+#
+#
+# def h_normal_gpytorch(s):
+#     """ Entropy of a normal distribution """
+#     return torch.log(s * (2 * np.e * np.pi) ** 0.5)
+#
+#
+# def approx_mi_vec_gpytorch(mu, cov, thresholds):
+#     # Expectation Propagation
+#     mu1 = mu[:, 0]
+#     std1 = cov[:, 0, 0] ** 0.5
+#     mu2 = mu[:, 1]
+#     std2 = cov[:, 1, 1] ** 0.5
+#     rho = cov[:, 0, 1] / (std1 * std2)
+#
+#     std_sx = []
+#
+#     for j in range(len(thresholds) - 1):
+#         alpha_j = (thresholds[j] - mu2) / std2
+#         beta_j = (thresholds[j + 1] - mu2) / std2
+#         alpha_j = alpha_j.detach().numpy()
+#         beta_j = beta_j.detach().numpy()
+#
+#         c_j = norm.cdf(beta_j) - norm.cdf(alpha_j)
+#
+#         # \sigma(Y(X)|S(x')=j)
+#         b_phi_b = beta_j * norm.pdf(beta_j)
+#         b_phi_b[~np.isfinite(beta_j)] = 0.0
+#         a_phi_a = alpha_j * norm.pdf(alpha_j)
+#         a_phi_a[~np.isfinite(alpha_j)] = 0.0
+#
+#         alpha_j = torch.tensor(alpha_j)
+#         beta_j = torch.tensor(beta_j)
+#
+#         mu_cond = mu1 - std1 * rho / torch.tensor(c_j) * (
+#             torch.tensor(norm.pdf(beta_j)) - torch.tensor(norm.pdf(alpha_j))
+#         )
+#         var_cond = (
+#             mu1 ** 2
+#             - 2
+#             * mu1
+#             * std1
+#             * (
+#                 rho
+#                 / torch.tensor(c_j)
+#                 * (torch.tensor(norm.pdf(beta_j)) - torch.tensor(norm.pdf(alpha_j)))
+#             )
+#             + std1 ** 2
+#             * (
+#                 1.0
+#                 - (rho ** 2 / torch.tensor(c_j))
+#                 * (torch.tensor(b_phi_b) - torch.tensor(a_phi_a))
+#             )
+#             - mu_cond ** 2
+#         )
+#
+#         std_sx_j = var_cond ** 0.5
+#
+#         std_sx.append(std_sx_j)
+#
+#     # Entropy
+#     h = h_normal_gpytorch(std1)
+#
+#     for j in range(len(thresholds) - 1):
+#         p_j = norm(mu2.detach().numpy(), std2.detach().numpy()).cdf(
+#             thresholds[j + 1]
+#         ) - norm(mu2.detach().numpy(), std2.detach().numpy()).cdf(thresholds[j])
+#         print("pj ", p_j)
+#         dec = torch.tensor(p_j) * h_normal_gpytorch(std_sx[j])
+#         h[p_j > 0.0] -= dec[p_j > 0.0]
+#
+#     return h

excursion/acquisition/base.py

@@ -0,0 +1,28 @@
+from collections import defaultdict
+
+
+class AcquisitionFunction(object):
+    """All acquisition functions used by the library should implement this interface
+    to be used by the optimizer object."""
+    def acquire(self, gp, thresholds, meshgrid):
+        raise NotImplemented
+
+    def set_params(self, **params):
+        """
+        Set the parameters of this acquisition function.
+        Parameters
+        ----------
+        **params : dict
+            Function parameters.
+        Returns
+        -------
+        self : object
+            Function instance.
+        """
+        if not params:
+            # Simple optimization to gain speed (inspect is slow)
+            return self
+        for key, value in params.items():
+            setattr(self, key, value)
+
+        return self