Feature/scikit pr code

.gitignore

@@ -1,2 +1,2 @@
 *.pyc
-__pycache__
+__pycache__

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/__init__.py

@@ -0,0 +1,7 @@
+from .skPES import SKPES
+from .base import AcquisitionFunction
+
+__all__ = [
+    "SKPES",
+    "AcquisitionFunction",
+]

excursion/acquisition/base.py

@@ -0,0 +1,30 @@
+from collections import defaultdict
+import numpy as np
+
+class AcquisitionFunction(object):
+    """
+    All acquisition functions used by this library should implement this abstract interface in order to be compatible
+    with the Optimizer class.
+    """
+    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

excursion/acquisition/skPES.py

@@ -0,0 +1,37 @@
+from .base import AcquisitionFunction
+from .utils import info_gain
+import numpy as np
+import os
+
+
+class SKPES(AcquisitionFunction):
+    def __init__(self, batch=False):
+        self.batch = batch
+        self.acq_vals = None
+
+    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
+
+        """
+
+        self.acq_vals = self._acquire(gp, thresholds, X_pointsgrid)
+
+        X_train = gp.X_train_.tolist()
+        for i, cacq in enumerate(X_pointsgrid[np.argsort(self.acq_vals)]):
+            if cacq.tolist() not in X_train:
+                newx = cacq
+                return newx
+
+    def _acquire(self, gp, thresholds, X_pointsgrid):
+        try:
+            from joblib import Parallel, delayed
+            nparallel = int(os.environ.get('EXCURSION_NPARALLEL', os.cpu_count()))
+            result = Parallel(nparallel)(
+                delayed(info_gain)(x_candidate, gp, thresholds, X_pointsgrid) for x_candidate in X_pointsgrid)
+            return np.asarray(result)
+        except ImportError:
+            return np.array([info_gain(x_candidate, gp, thresholds, X_pointsgrid) for x_candidate in X_pointsgrid])

excursion/acquisition/utils.py

@@ -0,0 +1,76 @@
+import numpy as np
+from scipy.linalg import cho_solve
+from scipy.stats import norm
+
+
+def approx_mi_vec(mu, cov, thresholds):
+    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
+        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
+
+        mu_cond = mu1 - std1 * rho / c_j * (norm.pdf(beta_j) - norm.pdf(alpha_j))
+        var_cond = (mu1 ** 2 - 2 * mu1 * std1 * (rho / c_j * (norm.pdf(beta_j) - norm.pdf(alpha_j))) +
+                    std1 ** 2 * (1. - (rho ** 2 / c_j) * (b_phi_b - a_phi_a)) -
+                    mu_cond ** 2)
+        std_sx_j = var_cond ** 0.5
+
+        std_sx.append(std_sx_j)
+
+    # Entropy
+    CONSTANT = (2 * np.e * np.pi) ** 0.5
+
+    h = np.log(std1 * CONSTANT)
+    for j in range(len(thresholds) - 1):
+        p_j = norm(mu2, std2).cdf(thresholds[j+1]) - norm(mu2, std2).cdf(thresholds[j])
+        dec = p_j * np.log(std_sx[j] * CONSTANT)
+        h[p_j > 0.0] -= dec[p_j > 0.0]
+
+    return h
+
+
+def info_gain(x_candidate, gp, thresholds, meanX):
+    n_samples = len(meanX)
+    X_all = np.concatenate([np.array([x_candidate]), meanX]).reshape(1 + n_samples, -1)
+    K_trans_all = gp.kernel_(X_all, gp.X_train_)
+    y_mean_all = K_trans_all.dot(gp.alpha_) + gp._y_train_mean
+    v_all = cho_solve((gp.L_, True), K_trans_all.T)
+
+    mus = np.zeros((n_samples, 2))
+    mus[:, 0] = y_mean_all[0]
+    mus[:, 1] = y_mean_all[1:]
+
+    covs = np.zeros((n_samples, 2, 2))
+    c = gp.kernel_(X_all[:1], X_all)
+    covs[:, 0, 0] = c[0, 0]
+    covs[:, 1, 1] = c[0, 0]
+    covs[:, 0, 1] = c[0, 1:]
+    covs[:, 1, 0] = c[0, 1:]
+
+    x_train_len = len(gp.X_train_)
+    K_trans_all_repack = np.zeros((n_samples, 2, x_train_len))
+    K_trans_all_repack[:, 0, :] = K_trans_all[0, :]
+    K_trans_all_repack[:, 1, :] = K_trans_all[1:]
+    v_all_repack = np.zeros((n_samples, x_train_len, 2))
+    v_all_repack[:, :, 0] = v_all[:, 0]
+    v_all_repack[:, :, 1] = v_all[:, 1:].T
+    covs -= np.einsum('...ij,...jk->...ik', K_trans_all_repack, v_all_repack)
+
+    mi = approx_mi_vec(mus, covs, thresholds)
+    mi[~np.isfinite(mi)] = 0.0
+
+    return -np.mean(mi)