Merge pull request #4 from IBM/poi

Poi

Author: ordavidov

configs/dir_configs.yaml

@@ -0,0 +1,10 @@
+local:
+  buckets: 
+    inputs: ''
+    inputs_dest: ''
+    objectives: ''
+    objectives_dest: ''
+    data: ''
+    data_dest: ''
+    solutions: ''
+    solutions_dest: ''

doframework/api.py

@@ -322,9 +322,9 @@ def run(generate_user_solution, configs_file, **kwargs):
     args = Args(objectives, datasets, feasibility_regions, run_mode, distribute, mcmc, logger, mock, after_idle_for, rayvens_logs, alg_num_cpus, data_num_cpus)
 
     if args.run_mode == 'operator':
-        ray.init(address='auto')
+        ray.init(address='auto',ignore_reinit_error=True)
     else:
-        ray.init()
+        ray.init(ignore_reinit_error=True)
     rayvens.init(mode=args.run_mode ,release=(not args.rayvens_logs))
 
     if args.logger: print('({}) INFO ... Running simulation with args objectives={o} datasets={s} feasibility_regions={r} distribute={d} run_mode={m} logger={l}'.format('root', 

doframework/core/gp.py

@@ -14,13 +14,73 @@
 # limitations under the License.
 #
 
+from typing import List
 import numpy as np
 from scipy.stats import gaussian_kde
 from scipy.integrate import quad
-from GPy.kern import RBF
-from GPy.models import GPRegression
-from GPy.core.parameterization.priors import Gamma
-from GPy.inference.mcmc import HMC
+import GPy
+
+def plot_posteriors(samples: np.array, labels: List[str], **kwargs):    
+
+    '''
+    Plot hyper-parameter posteriors following HMC. By convention, first column of samples is RBF kernel variance chain. 
+    The last column of samples is the RBF kernel noise. The middle columns of samples are ARD length scale parameters.
+    
+            Parameters:
+                    samples (np.array): hyper-parameter HMC chains.
+                    labels (List[str]): sample column labels for graph legends.
+    '''        
+    
+    import matplotlib.pyplot as plt
+
+    num_of_plots = 3
+    num_of_posteriors = samples.shape[-1]
+    dim = num_of_posteriors - num_of_plots + 1
+
+    assert len(labels) == num_of_posteriors, \
+    f'Expected {num_of_posteriors} plot labels, but received {len(labels)}.'
+
+    fig_width = kwargs['fig_width'] if 'fig_width' in kwargs else 10
+    fig_length = kwargs['fig_length'] if 'fig_length' in kwargs else 3
+
+    _, axs = plt.subplots(num_of_plots,1,figsize=(fig_width,num_of_plots*fig_length))
+    cmap = plt.cm.get_cmap(name='Accent',lut=num_of_posteriors)
+
+    for i, J in zip(range(num_of_plots),[[0],[1+j for j in range(dim)],[dim+1]]):
+
+        modals = []
+
+        for j in J:
+            
+            s = samples[:,j]
+            xmin = s.min()*0.9
+            xmax = s.max()*1.1
+            xs = np.linspace(xmin,xmax,1000)
+
+            if ('kdes' in kwargs) and (j in kwargs['kdes']):
+                kde = kwargs['kdes'][j]
+            else:
+                kde = gaussian_kde(s)
+
+            if ('modals' in kwargs) and (j in kwargs['modals']):
+                modal = kwargs['modals'][j]
+            else:
+                density = kde.pdf(xs)
+                argmax = np.argmax(density)
+                modal = xs[argmax]
+
+            modals.append(modal)            
+            modal_density = kde.pdf(modal)[0]
+
+            axs[i].plot([modal,modal],[0,modal_density],ls='--',color=cmap(j)) 
+            axs[i].plot(xs,kde(xs),label=labels[j],lw=3,color=cmap(j))
+
+        ticks = np.sort(np.hstack([np.around(np.array(modals),2), axs[i].get_xticks()]))
+        axs[i].set_xticks(ticks[1:-1]) # ticks[1:-1]
+        axs[i].tick_params(axis="x", rotation=90, labelsize=12)
+        axs[i].legend()
+
+    plt.tight_layout()
 
 def find_modal(samples, linspace_num: int=1000):
 
@@ -35,69 +95,78 @@ def find_modal(samples, linspace_num: int=1000):
         argmax = np.argmax(density)
         modal = xs[argmax]
 
-        return modal
+        return modal, kde
 
     except:
 
-        return None
+        return None, None
 
 def gp_model(X: np.array, 
              y: np.array, 
              is_mcmc: bool=False, 
              num_samples: int=1000,
              hmc_iters: int=2,
-             linspace_num: int=1000) -> GPRegression:
+             plot_kernel_posteriors: bool=False,
+             linspace_num: int=1000) -> GPy.models.GPRegression:
 
     dim = X.shape[-1]
 
     if is_mcmc:        
 
-        factor = 10.0 # TODO: clever factor for numerical issues in HMC train
+        factor = 10.0 # factor for numerical issues in HMC train
 
-        kern = RBF(input_dim=dim, ARD=True)
-        model = GPRegression(factor*X,y,kernel=kern.copy())    
+        kern = GPy.kern.RBF(input_dim=dim, ARD=True)
+        model = GPy.models.GPRegression(factor*X,y,kernel=kern.copy())    
 
-        # TODO: automate prior for RBF variance
-        model.kern.variance.set_prior(Gamma.from_EV(0.1,0.1),warning=False)
+        # automate prior for RBF variance
+        model.kern.variance.set_prior(GPy.priors.Gamma.from_EV(0.1,0.1),warning=False)
 
         lengthscales = {}
+        kdes = {}
         for i in range(dim):
             kde = gaussian_kde(X[:,i])
+            kdes[i+1] = kde
             mean = quad(lambda x: x * kde.pdf(x), a=-np.inf, b=np.inf)[0]
             var = quad(lambda x: x**2 * kde.pdf(x), a=-np.inf, b=np.inf)[0] - mean**2
             lengthscales[i] = np.sqrt(var)
-            model.kern.lengthscale[[i]].set_prior(Gamma.from_EV(lengthscales[i],lengthscales[i]/2),warning=False) # data variance as length scale
+            model.kern.lengthscale[[i]].set_prior(GPy.priors.Gamma.from_EV(lengthscales[i],lengthscales[i]/2),warning=False)
 
-        hmc = HMC(model)
+        hmc = GPy.inference.mcmc.HMC(model)
         samples = hmc.sample(num_samples=num_samples,hmc_iters=hmc_iters)
 
-        modals = {}          
+        modals = {}
+        kdes = {}
         for i in range(samples.shape[-1]):
-            modal = find_modal(samples[:,i],linspace_num)
-            if modal is not None:
+            modal, kde = find_modal(samples[:,i],linspace_num)
+            if (modal is not None) and (kde is not None):
                 modals[i] = modal
+                kdes[i] = kde
 
-        kern = RBF(input_dim=dim, ARD=True)
-        model = GPRegression(X,y,kernel=kern.copy())    
+        kern = GPy.kern.RBF(input_dim=dim, ARD=True)
+        model = GPy.models.GPRegression(X,y,kernel=kern.copy())    
 
         if (0 in modals) and (dim-1 in modals):            
             model.rbf.variance = modals[0]/factor**2
             model.Gaussian_noise.variance = modals[dim-1]/factor**2
         else:
-            model = None
-            
+            raise ValueError('HMC failed. Possible unsuitable priors on kernels parameters leading to repetative samples.')
+
         for i in range(dim):
             if i in modals:
                 model.rbf.lengthscale[i] = modals[1+i]/factor
             else:
                 model.rbf.lengthscale[i] = lengthscales[i]
 
+        if plot_kernel_posteriors:
+            labels = ['RBF kernel variance']+ [f'RBF kernel lengthscale[x{i}]' for i in range(dim)] + ['RBF kernel noise']
+            plot_posteriors(samples, labels, modals=modals, kdes=kdes)
+
     else:
 
-        kern = RBF(input_dim=dim, ARD=True)
-        model = GPRegression(X,y,kernel=kern.copy())    
-        
+        kern = GPy.kern.RBF(input_dim=dim, ARD=True)
+        model = GPy.models.GPRegression(X,y,kernel=kern.copy())    
+
         # model.optimize_restarts(num_restarts=10,optimizer='lbfgs',verbose=False)
         model.optimize(optimizer='lbfgs',messages=False)
-        
-    return model
+
+    return model

doframework/core/hit_and_run.py

@@ -174,11 +174,12 @@ def in_domain(xs, A: np.array, R: Optional[float]=None, **kwargs) -> np.array:
     '''
     Check whether xs are inside the intersection of a convex polytope and a ball of radius R.
     The convex polytope is given by the matrix A such that Ax<=0 defines it.
+    When the radius is not specified, this restriction is dropped.
     
             Parameters:
                     xs (np.array): points.
                     A (np.array): matrix defining a convex polytope Ax <= 0.
-                    R (float): radius.
+                    R (float): radius (default: None).
                     
             Returns:
                     a boolean numpy array that indicates whether the points are inside the polytope.

doframework/core/inputs.py

@@ -76,7 +76,7 @@ def parse_vertex_num(sim_input: dict) -> int:
 
     return num
 
-def get_configs(configs_file, is_logger: bool=True):
+def get_configs(configs_file, is_logger: bool=False):
 
     with open(configs_file,'r') as file:
         try:
@@ -86,7 +86,7 @@ def get_configs(configs_file, is_logger: bool=True):
             if is_logger:
                 print('({}) ERROR ... Could not load configs yaml. Check your path.'.format('root'))
                 print(e)
-            raise e        
+            raise e   
 
 def legit_configs(configs):
 

doframework/core/poi.py

@@ -0,0 +1,125 @@
+#
+# Copyright IBM Corporation 2022
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+from typing import List
+
+import numpy as np
+import pandas as pd
+from scipy.stats import multivariate_normal
+from GPy.models import GPRegression
+
+from doframework.core.utils import order_stats
+
+def plot_joint_distribution(samples: np.array, **kwargs):
+
+    import seaborn as sns
+    
+    assert samples.shape[-1]==2, 'Array must be of dimension Nx2. Array is a draw from a bivariate distribution.'
+    
+    cols = kwargs['cols'] if 'cols' in kwargs else [f'x{i}' for i in range(2)]
+
+    if 'x_min' not in kwargs or 'x_max' not in kwargs:
+        x_min = samples.min(axis=0).min()*0.5
+        x_max = samples.max(axis=0).max()*1.5
+    else:
+        x_min = kwargs['x_min']
+        x_max = kwargs['x_max']
+
+    lims = np.array([x_min,x_max])[:,None]
+
+    df = pd.DataFrame(samples,columns=cols)
+    dl = pd.DataFrame(np.hstack([lims,lims]),columns=cols)
+
+    sns.set(style="white", color_codes=True)
+    sns.jointplot(data=df, x=cols[0], y=cols[1], kind="hex", xlim=lims, ylim=lims)
+    sns.lineplot(data=dl, x=cols[0], y=cols[1])
+
+class POI(object):
+    '''
+    Class for probability of improvement outcomes.
+    '''
+    
+    def __init__(self, point: np.array, probability: float, **kwargs):
+        
+        self.point = point
+        assert all([probability>=0.0,probability<=1.0]), f'Probability value should be in [0,1]. Received {probability:.2f}.'
+        self.probability = probability
+                                    
+        self.upper_bound = kwargs['upper_bound'] if 'upper_bound' in kwargs else True
+        self.reference = kwargs['reference'] if 'reference' in kwargs else np.array([])
+        self.threshold = kwargs['threshold'] if 'threshold' in kwargs else None
+            
+    def __repr__(self):
+        return 'POI('+''.join([f'point={self.point},',
+                                f' probability={self.probability},',
+                                f' upper_bound={self.upper_bound}',
+                                ','*any([self.reference.size > 0]),
+                                f' reference={self.reference}'*(self.reference.size > 0),
+                               ','*any([self.threshold is not None]),
+                                f' threshold={self.threshold}'*(self.threshold is not None)])+')'
+
+def probability_of_improvement(solutions: np.array, references: np.array, model: GPRegression,
+                               sample_num: int=100000, is_constraint: bool=False, upper_bound: bool=True, plot_joint_gaussians: bool=False,
+                               **kwargs) -> List[POI]:
+    
+    sols = np.atleast_2d(solutions)
+    d = sols.shape[-1]
+    is_minimum = not upper_bound
+
+    if is_constraint:
+
+        refs = np.atleast_2d(references.flatten()).T
+
+    else:
+
+        refs = np.atleast_2d(references)
+
+    ref_num = refs.shape[0]
+    ref_dim = refs.shape[-1]
+
+    assert ref_dim == d or ref_dim == 1, \
+    'Input reference row dimension must either be:\n(1) identical to solution row dimension (POI for objective target, is_constraint=False)\n(2) or equal to 1 (POI for constraint satisfaction, is_constraint=True).\nYour input has inferred dimension {} for solution vectors and inferred dimension {} for reference vectors (is_constraint={}).'.format(d,ref_dim,is_constraint)
+
+    sols_rep = np.tile(sols, (1,ref_num)).reshape(ref_num*sols.shape[0],sols.shape[-1])
+    refs_rep = np.tile(refs, (sols.shape[0],1))
+
+    N = sols_rep.shape[0]*(sols_rep.shape[0] == refs_rep.shape[0])
+
+    pois = []
+
+    for i in range(N):
+
+        if is_constraint:
+            mu, cov = model.predict(np.vstack([sols_rep[i]]),full_cov=True)
+            samples = multivariate_normal(mean=mu.flatten(),cov=cov).rvs(size=sample_num)
+            samples = np.hstack([samples[:,None],np.tile(refs_rep[i:i+1],(samples.size,1))])
+        else:
+            mu, cov = model.predict(np.vstack([sols_rep[i],refs_rep[i]]),full_cov=True)
+            samples = multivariate_normal(mean=mu.flatten(),cov=cov).rvs(size=sample_num)
+
+        if is_constraint:            
+            pois.append(POI(sols_rep[i],order_stats(samples,is_minimum),upper_bound=upper_bound,threshold=refs_rep[i]))
+        else:
+            pois.append(POI(sols_rep[i],order_stats(samples,is_minimum),upper_bound=upper_bound,reference=refs_rep[i]))
+
+        if plot_joint_gaussians and not is_constraint:
+
+            kwargs = {'cols': ['f({})'.format(np.around(sols_rep[i],2)),
+                               'f({})'.format(np.around(refs_rep[i]),2)]}
+            
+            plot_joint_distribution(samples=samples, **kwargs)
+
+    return pois

doframework/core/pwl.py

@@ -1,3 +1,19 @@
+#
+# Copyright IBM Corporation 2022
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
 import itertools
 import numpy as np
 from numpy import linalg