Add github workflow tests for HiOpBBPy (#716)

* test CI

* remove ipopt as a dependency

* update some parameters

* address Cosmin's comments

* test CI

* address comments

* fix bug, turn on CI

* save y_opt from 1000 runs

* address comments

* moving saved data from the base hiop directory in order to avoid overcrowding, obtaining the path of the python __file__ in order that we can load the saved data according to the python file and not from where the file is run from. This enables us to be able to run (from hiop base dir) python src/Drivers/hiopbbpy/BODriver.py and (from hiop/src/Drivers/hiopbbpy) python BODriver.py. Previously there was an error when I tried running the latter case.

---------

Co-authored-by: Tucker Hartland <tucker.hartland@gmail.com>

Author: nychiang

.github/workflows/test_hiopbbpy.yml

@@ -0,0 +1,32 @@
+name: test_hiopbbpy
+
+on: [push]
+
+jobs:
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest]
+
+    steps:
+      - name: Checkout code
+        uses: actions/checkout@v4
+
+      - name: Set up Python
+        uses: actions/setup-python@v4  # Set up Python environment
+        with:
+          python-version: '3.13'
+
+
+      - name: Install HiOpBBPy and its dependencies
+        run: |
+          python -m pip install --upgrade pip
+          pip install .
+
+      - name: Run Tests
+        run: |
+          python src/Drivers/hiopbbpy/BODriver.py 10
+        continue-on-error: false
+          
+

pyproject.toml

@@ -1,3 +1,3 @@
 [build-system]
-requires = ["setuptools", "wheel", "numpy", "smt", "cyipopt"]
+requires = ["setuptools", "wheel", "numpy", "smt"]
 build-backend = "setuptools.build_meta"

setup.py

@@ -12,14 +12,14 @@
 
 metadata = dict(
         name="hiopbbpy",
-        version="0.0.2",
+        version="0.0.3",
         description="HiOp black box optimization (hiopbbpy)",
         author="Tucker hartland et al.",
         author_email="hartland1@llnl.gov",
         license="BSD-3",
         packages=find_packages(where="src"),
         package_dir={"": "src"},
-        install_requires=["smt", "cyipopt"],
+        install_requires=["smt"],
         python_requires=">=3.9",
         zip_safe=False,
         url="https://github.com/LLNL/hiop",

src/Drivers/hiopbbpy/BODriver.py

@@ -5,8 +5,13 @@
         2) define a Kriging-based Gaussian-process (smt backend)
            trained on said data
         3) determine the minimizer via BOAlgorithm
+
+Authors:    Tucker Hartland <hartland1@llnl.gov>
+            Nai-Yuan Chiang <chiang7@llnl.gov>
 """
 
+import sys
+import os
 import numpy as np
 import matplotlib.pyplot as plt
 import warnings
@@ -17,34 +22,106 @@
 from hiopbbpy.problems import BraninProblem
 
 
-### parameters
-n_samples = 5 # number of the initial samples to train GP
-theta = 1.e-2 # hyperparameter for GP kernel
+# Get user input for the number of repetitions from command-line arguments
+if len(sys.argv) != 2 or int(sys.argv[1]) < 0:
+    num_repeat = 1
+else:
+    num_repeat = int(sys.argv[1])
 
-nx = 2 # dimension of the problem
+### parameters
+n_samples = 5  # number of the initial samples to train GP
+theta = 1.e-2  # hyperparameter for GP kernel
+nx = 2         # dimension of the problem
 xlimits = np.array([[-5, 5], [-5, 5]]) # bounds on optimization variable
 
-#problem = LpNormProblem(nx, xlimits)
-problem = BraninProblem()
-
-print(problem.name, " problem")
-
-### initial training set
-x_train = problem.sample(n_samples)
-y_train = problem.evaluate(x_train)
-
-# Define the GP surrogate model
-gp_model = smtKRG(theta, xlimits, nx)
-gp_model.train(x_train, y_train)
-
-acquisition_types = ["LCB", "EI"]
-for acquisition_type in acquisition_types:
-    print("acquisition type: ", acquisition_type)
-    
-    # Instantiate and run Bayesian Optimization
-    bo = BOAlgorithm(gp_model, x_train, y_train, acquisition_type = acquisition_type) #EI or LCB
-    bo.optimize(problem)
-    
-    # Retrieve optimal point
-    x_opt, y_opt = bo.getOptimalPoint()
-    print()
+### saved solutions --- from 1000 repetitions
+saved_min_obj = {"LpNorm": {"LCB": 0.0007586314501994839, "EI": 0.002094016049616341}, "Branin": {"LCB": 0.3979820338569908, "EI": 0.39789916461969455}}
+saved_mean_obj = {"LpNorm": {"LCB": 0.018774638321851504, "EI": 0.11583915178648867}, "Branin": {"LCB": 0.5079001079219421, "EI": 0.4377466109837465}}
+saved_max_obj = {"LpNorm": {"LCB": 0.0755173754382861, "EI": 0.4175676394969743}, "Branin": {"LCB": 1.107240543567082, "EI": 0.7522382699410031}}
+script_dir = os.path.dirname(__file__)
+saved_yopt = np.load(script_dir + "/yopt_20iter_1000run.npy",allow_pickle=True).item()
+
+prob_type_l = ["LpNorm", "Branin"]
+acq_type_l = ["LCB", "EI"]
+
+mean_obj = {}
+max_obj = {}
+min_obj = {}
+y_opt = {}
+
+retval = 0
+for prob_type in prob_type_l:
+   print()
+   if prob_type == "LpNorm":
+      problem = LpNormProblem(nx, xlimits)
+   else:
+      problem = BraninProblem()
+
+   if prob_type not in mean_obj:
+      mean_obj[prob_type] = {}
+      max_obj[prob_type] = {}
+      min_obj[prob_type] = {}
+      y_opt[prob_type] = {}
+
+   for acq_type in acq_type_l:
+      if acq_type not in mean_obj[prob_type]:
+         mean_obj[prob_type][acq_type] = 0
+         max_obj[prob_type][acq_type] = -np.inf
+         min_obj[prob_type][acq_type] = np.inf
+         y_opt[prob_type][acq_type] = np.zeros(num_repeat)
+
+      print("Problem name: ", problem.name)
+      print("Acquisition type: ", acq_type)
+   
+      for n_repeat in range(num_repeat):
+         print("Run: ", n_repeat, "/", num_repeat)
+         ### initial training set
+         x_train = problem.sample(n_samples)
+         y_train = problem.evaluate(x_train)
+
+         ### Define the GP surrogate model
+         gp_model = smtKRG(theta, xlimits, nx)
+         gp_model.train(x_train, y_train)
+   
+         # Instantiate and run Bayesian Optimization
+         bo = BOAlgorithm(gp_model, x_train, y_train, acquisition_type = acq_type) #EI or LCB
+         bo.optimize(problem)
+         
+         # Retrieve optimal objec
+         y_opt[prob_type][acq_type][n_repeat] = bo.getOptimalObjective()
+         
+         mean_obj[prob_type][acq_type] += y_opt[prob_type][acq_type][n_repeat]
+         max_obj[prob_type][acq_type] = max(max_obj[prob_type][acq_type], y_opt[prob_type][acq_type][n_repeat])
+         min_obj[prob_type][acq_type] = min(min_obj[prob_type][acq_type], y_opt[prob_type][acq_type][n_repeat])
+
+#if(num_repeat >= 1000 ):
+#   np.save("yopt_20iter_1000run.npy", y_opt)
+
+# Define percentiles
+left_percentile = 1  # 5 or 1
+right_percentile = 100 - left_percentile  # 95 or 99
+
+print("Summary:")
+for prob_type in prob_type_l:
+   for acq_type in acq_type_l:
+      mean_obj[prob_type][acq_type] /= num_repeat
+      print("(Min,Mean,Max) Opt.Obj for", prob_type, "-", acq_type, ":\t(", min_obj[prob_type][acq_type], ",",mean_obj[prob_type][acq_type], ",", max_obj[prob_type][acq_type], ")")
+   
+      ### verify the results with the saved results
+      left_value = np.percentile(saved_yopt[prob_type][acq_type], left_percentile)
+      right_value = np.percentile(saved_yopt[prob_type][acq_type], right_percentile)
+
+      is_failed = (y_opt[prob_type][acq_type] < left_value) | (y_opt[prob_type][acq_type] > right_value)
+      num_fail = np.sum(is_failed)
+
+      # currently ci test is only applied to num_repeat == 10
+      if num_fail >= 3 and num_repeat == 10:
+         print(num_fail, "test(s) fail(s):", y_opt[prob_type][acq_type][is_failed])
+         print("Recorded (min, mean, max): (", saved_min_obj[prob_type][acq_type], ",", saved_mean_obj[prob_type][acq_type], ",", saved_max_obj[prob_type][acq_type], ")")
+         retval = 1
+
+sys.exit(retval)
+
+
+
+

src/Drivers/hiopbbpy/yopt_20iter_1000run.npy

@@ -43,7 +43,7 @@ def evaluate(self, x : np.ndarray) -> np.ndarray:
         y_min = y_data[np.argmin(y_data[:, 0])]
 
         pred = self.gpsurrogate.mean(x)
-        sig = self.gpsurrogate.variance(x)
+        sig = np.sqrt(self.gpsurrogate.variance(x))
 
         retval = []
         if sig.size == 1 and np.abs(sig) > 1e-12:

src/hiopbbpy/opt/acquisition.py

@@ -8,6 +8,7 @@
 import numpy as np
 from numpy.random import uniform
 from scipy.optimize import minimize
+from scipy.stats import qmc
 from ..surrogate_modeling.gp import GaussianProcess
 from .acquisition import LCBacquisition, EIacquisition
 from ..problems.problem import Problem
@@ -18,6 +19,7 @@ def __init__(self):
         self.acquisition_type = "LCB" # Type of acquisition function (default = "LCB")
         self.xtrain = None            # Training data
         self.ytrain = None            # Training data
+        self.prob   = None            # Problem structure
         self.n_iter = 20              # Maximum number of optimization steps
         self.n_start = 10             # estimating acquisition global optima by determining local optima n_start times and then determining the discrete max of that set
         self.q = 1                    # batch size
@@ -48,11 +50,15 @@ def getOptimizationHistory(self):
         y_hist = np.array(self.y_hist, copy=True)
         return x_hist, y_hist
 
-    # Method to return the optimal solution and objective
+    # Method to return the optimal solution 
     def getOptimalPoint(self):
         x_opt = np.array(self.x_opt, copy=True)
+        return x_opt
+
+    # Method to return the optimal objective
+    def getOptimalObjective(self):
         y_opt = np.array(self.y_opt, copy=True)
-        return x_opt, y_opt
+        return y_opt
 
 # A subclass of BOAlgorithmBase implementing a full Bayesian Optimization workflow
 class BOAlgorithm(BOAlgorithmBase):
@@ -96,8 +102,11 @@ def _find_best_point(self, x_train, y_train, x0 = None):
         for ii in range(self.n_start):
             success = False
             # Generate random starting point if x0 is not provided
-            if x0 is None:
+            if x0 is None and self.prob is not None:
+                x0 = self.prob.sample(1)
+            else:
                 x0 = np.array([uniform(b[0], b[1]) for b in self.bounds])
+                
             xopt, yout, success = self.acqf_minimizer_callback(acqf_callback, x0)
 
             if success:
@@ -118,6 +127,7 @@ def set_options(self, options):
 
     # Method to perform Bayesian optimization
     def optimize(self, prob:Problem):
+      self.problem = prob
       x_train = self.xtrain
       y_train = self.ytrain
 
@@ -162,6 +172,8 @@ def optimize(self, prob:Problem):
           print(f"Optimal at BO iteration: {self.idx_opt-n_init_sample+1} ")
 
       print(f"Optimal point: {self.x_opt.flatten()}, Optimal value: {self.y_opt}")
+      print()
+
 
 
 # Find the minimum of the input objective `fun`, using the minimize function from SciPy.