diff --git a/.github/workflows/test_hiopbbpy.yml b/.github/workflows/test_hiopbbpy.yml index e4d9289fe..32bbf7f0c 100644 --- a/.github/workflows/test_hiopbbpy.yml +++ b/.github/workflows/test_hiopbbpy.yml @@ -14,19 +14,39 @@ jobs: uses: actions/checkout@v4 - name: Set up Python - uses: actions/setup-python@v4 # Set up Python environment + uses: actions/setup-python@v4 with: - python-version: '3.13' + python-version: "3.13" + - name: Set up Conda + uses: conda-incubator/setup-miniconda@v2 + with: + auto-activate-base: true + channels: conda-forge + + - name: Install IPOPT via conda + run: | + conda init bash + conda install -c conda-forge cyipopt pkg-config + export LD_LIBRARY_PATH=$CONDA_PREFIX/lib:$LD_LIBRARY_PATH - - name: Install HiOpBBPy and its dependencies + - name: Verify cyipopt installation run: | + source ~/.bashrc + conda activate base + python -c "import cyipopt; print('cyipopt OK')" + + - name: Install HiOpBBPy + run: | + source ~/.bashrc + conda activate base + export SKIP_CYIPOPT=1 python -m pip install --upgrade pip pip install . - name: Run Tests run: | - python src/Drivers/hiopbbpy/BODriver.py 10 + source ~/.bashrc + conda activate base + python src/Drivers/hiopbbpy/BODriverCI.py 10 continue-on-error: false - - diff --git a/setup.py b/setup.py index 9fbaed53d..c211af059 100644 --- a/setup.py +++ b/setup.py @@ -6,20 +6,26 @@ ''' import sys +import os import numpy as np from setuptools import setup, find_packages +install_requires = ["smt"] + +if os.getenv("SKIP_CYIPOPT", "0") != "1": + install_requires.append("cyipopt") + metadata = dict( name="hiopbbpy", - version="0.0.3", + version="0.0.4", 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"], + install_requires=install_requires, python_requires=">=3.9", zip_safe=False, url="https://github.com/LLNL/hiop", diff --git a/src/Drivers/hiopbbpy/BODriver.py b/src/Drivers/hiopbbpy/BODriverCI.py similarity index 93% rename from src/Drivers/hiopbbpy/BODriver.py rename to src/Drivers/hiopbbpy/BODriverCI.py index d05f53161..2bcc386f2 100644 --- a/src/Drivers/hiopbbpy/BODriver.py +++ b/src/Drivers/hiopbbpy/BODriverCI.py @@ -70,6 +70,15 @@ min_obj[prob_type][acq_type] = np.inf y_opt[prob_type][acq_type] = np.zeros(num_repeat) + options = { + 'acquisition_type': acq_type, + 'bo_maxiter': 20, + 'opt_solver': 'SLSQP', #"SLSQP" "IPOPT" + 'solver_options': { + 'maxiter': 200 + } + } + print("Problem name: ", problem.name) print("Acquisition type: ", acq_type) @@ -84,8 +93,8 @@ 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) + bo = BOAlgorithm(problem, gp_model, x_train, y_train, options = options) + bo.optimize() # Retrieve optimal objec y_opt[prob_type][acq_type][n_repeat] = bo.getOptimalObjective() diff --git a/src/Drivers/hiopbbpy/BODriverEX.py b/src/Drivers/hiopbbpy/BODriverEX.py new file mode 100644 index 000000000..daf3b2720 --- /dev/null +++ b/src/Drivers/hiopbbpy/BODriverEX.py @@ -0,0 +1,95 @@ +""" + Code description: + for a 2D example LpNormProblem + 1) randomly sample training points + 2) define a Kriging-based Gaussian-process (smt backend) + trained on said data + 3) determine the minimizer via BOAlgorithm + +Authors: Tucker Hartland + Nai-Yuan Chiang +""" + +import sys +import os +import numpy as np +import matplotlib.pyplot as plt +import warnings +warnings.filterwarnings("ignore") +from LpNormProblem import LpNormProblem +from hiopbbpy.surrogate_modeling import smtKRG +from hiopbbpy.opt import BOAlgorithm +from hiopbbpy.problems import BraninProblem + + +# 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]) + +### 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 + +prob_type_l = ["LpNorm"] # ["LpNorm", "Branin"] +acq_type_l = ["LCB"] # ["LCB", "EI"] + +def con_eq(x): + return x[0] + x[1] - 4 + +def con_jac_eq(x): + return np.array([1.0, 1.0]) + +def con_ineq(x): + return x[0] - x[1] + +def con_jac_ineq(x): + return np.array([1.0, -1.0]) + +# only 'trust-constr' method supports vector-valued constraints +user_constraint = [{'type': 'ineq', 'fun': con_ineq, 'jac': con_jac_ineq}, + {'type': 'eq', 'fun': con_eq, 'jac': con_jac_eq}] + +retval = 0 +for prob_type in prob_type_l: + print() + if prob_type == "LpNorm": + problem = LpNormProblem(nx, xlimits, constraints=None) + else: + problem = BraninProblem(constraints=None) + problem.set_constraints(user_constraint) + + for acq_type in acq_type_l: + print("Problem name: ", problem.name) + print("Acquisition type: ", acq_type) + + ### 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) + + options = { + 'acquisition_type': acq_type, + 'bo_maxiter': 10, + 'opt_solver': 'IPOPT', #"SLSQP" "IPOPT" + 'solver_options': { + 'max_iter': 200, + 'print_level': 1 + } + } + + # Instantiate and run Bayesian Optimization + bo = BOAlgorithm(problem, gp_model, x_train, y_train, options = options) #EI or LCB + bo.optimize() + +sys.exit(retval) + + + + diff --git a/src/Drivers/hiopbbpy/LpNormProblem.py b/src/Drivers/hiopbbpy/LpNormProblem.py index 2ce3b2fc1..3084bd802 100644 --- a/src/Drivers/hiopbbpy/LpNormProblem.py +++ b/src/Drivers/hiopbbpy/LpNormProblem.py @@ -8,9 +8,9 @@ from hiopbbpy.problems.problem import Problem class LpNormProblem(Problem): - def __init__(self, ndim, xlimits, p=2.0): + def __init__(self, ndim, xlimits, p=2.0, constraints=None): name = "LpNormProblem" - super().__init__(ndim, xlimits, name=name) + super().__init__(ndim, xlimits, name=name, constraints=constraints) self.p = p def _evaluate(self, x): diff --git a/src/hiopbbpy/opt/__init__.py b/src/hiopbbpy/opt/__init__.py index 163941272..5692e3515 100644 --- a/src/hiopbbpy/opt/__init__.py +++ b/src/hiopbbpy/opt/__init__.py @@ -1,5 +1,6 @@ from .boalgorithm import (BOAlgorithmBase, BOAlgorithm) from .acquisition import (acquisition, LCBacquisition, EIacquisition) +from .optproblem import (IpoptProb) __all__ = [ "BOAlgorithmBase" @@ -7,4 +8,5 @@ "acquisition" "LCBacquisition" "EIacquisition" + "IpoptProb" ] diff --git a/src/hiopbbpy/opt/acquisition.py b/src/hiopbbpy/opt/acquisition.py index 9bc602b63..07a496e76 100644 --- a/src/hiopbbpy/opt/acquisition.py +++ b/src/hiopbbpy/opt/acquisition.py @@ -14,28 +14,41 @@ class acquisition(object): def __init__(self, gpsurrogate): assert isinstance(gpsurrogate, GaussianProcess) # add something here self.gpsurrogate = gpsurrogate + self.has_gradient = False # Abstract method to evaluate the acquisition function at x. def evaluate(self, x: np.ndarray) -> np.ndarray: raise NotImplementedError("Child class of acquisition should implement method evaluate") + # Abstract method to evaluate the gradient of acquisition function at x. + def eval_g(self, x: np.ndarray) -> np.ndarray: + raise NotImplementedError("Child class of acquisition should implement method evaluate") # A subclass of acquisition, implementing the Lower Confidence Bound (LCB) acquisition function. class LCBacquisition(acquisition): def __init__(self, gpsurrogate, beta=3.0): super().__init__(gpsurrogate) self.beta = beta + self.has_gradient = True # Method to evaluate the acquisition function at x. def evaluate(self, x : np.ndarray) -> np.ndarray: mu = self.gpsurrogate.mean(x) - sig = self.gpsurrogate.variance(x) - return mu - self.beta * np.sqrt(sig) + sig2 = self.gpsurrogate.variance(x) + return mu - self.beta * np.sqrt(sig2) + + def eval_g(self, x: np.ndarray) -> np.ndarray: + mu = self.gpsurrogate.mean(x) + sig2 = self.gpsurrogate.variance(x) + dsig2_dx = self.gpsurrogate.variance_gradient(x) + dmu_dx = self.gpsurrogate.mean_gradient(x) + return dmu_dx - 0.5 * self.beta * dsig2_dx / np.sqrt(sig2) # A subclass of acquisition, implementing the Expected improvement (EI) acquisition function. class EIacquisition(acquisition): def __init__(self, gpsurrogate): super().__init__(gpsurrogate) + self.has_gradient = True # Method to evaluate the acquisition function at x. def evaluate(self, x : np.ndarray) -> np.ndarray: @@ -47,12 +60,41 @@ def evaluate(self, x : np.ndarray) -> np.ndarray: retval = [] if sig.size == 1 and np.abs(sig) > 1e-12: - arg0 = (y_min - pred) / sig - retval = (y_min - pred) * norm.cdf(arg0) + sig * norm.pdf(arg0) + z = (y_min - pred) / sig + retval = (y_min - pred) * norm.cdf(z) + sig * norm.pdf(z) retval *= -1. elif sig.size == 1 and np.abs(sig) <= 1e-12: retval = 0.0 elif sig.size > 1: - NotImplementedError("TODO --- Not implemented yet!") + raise NotImplementedError("TODO --- Not implemented yet!") return retval + + def eval_g(self, x: np.ndarray) -> np.ndarray: + y_data = self.gpsurrogate.training_y + y_min = y_data[np.argmin(y_data[:, 0])] + + mean = self.gpsurrogate.mean(x) + sig2 = self.gpsurrogate.variance(x) + sig = np.sqrt(sig2) + + grad_EI = None + if sig.size == 1 and np.abs(sig) > 1e-12: + dmean_dx = self.gpsurrogate.mean_gradient(x) + dsig2_dx = self.gpsurrogate.variance_gradient(x) + dsig_dx = 0.5 * dsig2_dx / sig + + z = (y_min - mean) / sig + ncdf = norm.cdf(z) + npdf = norm.pdf(z) + EI = (y_min - mean) * ncdf + sig * npdf + + dz_dx = -dmean_dx / sig - (y_min - mean) * dsig_dx / sig**2 + grad_EI = -dmean_dx * ncdf + dsig_dx * npdf + grad_EI *= -1. + elif sig.size == 1 and np.abs(sig) <= 1e-12: + grad_EI = 0.0 + elif sig.size > 1: + raise NotImplementedError("TODO --- Not implemented yet!") + + return grad_EI \ No newline at end of file diff --git a/src/hiopbbpy/opt/boalgorithm.py b/src/hiopbbpy/opt/boalgorithm.py index ad882e4d7..081572f41 100644 --- a/src/hiopbbpy/opt/boalgorithm.py +++ b/src/hiopbbpy/opt/boalgorithm.py @@ -9,9 +9,11 @@ from numpy.random import uniform from scipy.optimize import minimize from scipy.stats import qmc +import warnings from ..surrogate_modeling.gp import GaussianProcess from .acquisition import LCBacquisition, EIacquisition from ..problems.problem import Problem +from .optproblem import IpoptProb # A base class defining a general framework for Bayesian Optimization class BOAlgorithmBase: @@ -20,7 +22,7 @@ def __init__(self): 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.bo_maxiter = 20 # Maximum number of Bayesian 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 # save some internal member train @@ -62,23 +64,49 @@ def getOptimalObjective(self): # A subclass of BOAlgorithmBase implementing a full Bayesian Optimization workflow class BOAlgorithm(BOAlgorithmBase): - def __init__(self, gpsurrogate, xtrain, ytrain, acquisition_type = "LCB"): + def __init__(self, prob:Problem, gpsurrogate:GaussianProcess, xtrain, ytrain, + user_grad = None, + options = None): super().__init__() + assert isinstance(gpsurrogate, GaussianProcess) - assert acquisition_type in ["LCB", "EI"] + self.setTrainingData(xtrain, ytrain) - self.setAcquisitionType(acquisition_type) + self.prob = prob self.gpsurrogate = gpsurrogate - self.n_iter = 20 - self.method = "SLSQP" self.bounds = self.gpsurrogate.get_bounds() - self.constraints = () - self.options = {"maxiter": 200} - self.acqf_minimizer_callback = None + self.fun_grad = None + + if options and 'bo_maxiter' in options: + self.bo_maxiter = options['bo_maxiter'] + assert self.bo_maxiter > 0, f"Invalid bo_maxiter: {self.bo_maxiter }" + + if options and 'solver_options' in options: + self.solver_options = options['solver_options'] + else: + self.solver_options = {"maxiter": 200} + + if options and 'acquisition_type' in options: + acquisition_type = options['acquisition_type'] + assert acquisition_type in ["LCB", "EI"], f"Invalid acquisition_type: {acquisition_type}" + else: + acquisition_type = "LCB" + self.setAcquisitionType(acquisition_type) + + if options and 'opt_solver' in options: + opt_solver = options['opt_solver'] + assert opt_solver in ["SLSQP", "IPOPT"], f"Invalid opt_solver: {opt_solver}" + else: + opt_solver = "SLSQP" + self.set_method(opt_solver) + + if user_grad: + self.fun_grad = user_grad + # Method to set up a callback function to minimize the acquisition function def _setup_acqf_minimizer_callback(self): - self.acqf_minimizer_callback = lambda fun, x0: pyminimize(fun, x0, self.method, self.bounds, self.constraints, self.options) + self.acqf_minimizer_callback = lambda fun, x0: minimizer(fun, x0, self.opt_solver, self.bounds, self.prob.constraints, self.solver_options) # Method to train the GP model def _train_surrogate(self, x_train, y_train): @@ -95,39 +123,44 @@ def _find_best_point(self, x_train, y_train, x0 = None): else: raise NotImplementedError("No implemented acquisition_type associated to"+self.acquisition_type) - acqf_callback = lambda x: float(np.array(acqf.evaluate(np.atleast_2d(x))).flat[0]) - + acqf_obj_callback = lambda x: float(np.array(acqf.evaluate(np.atleast_2d(x))).flat[0]) + acqf_callback = {'obj': acqf_obj_callback} + if acqf.has_gradient == True: + acqf_grad_callback = lambda x: np.array(acqf.eval_g(np.atleast_2d(x))) + acqf_callback['grad'] = acqf_grad_callback + x_all = [] y_all = [] for ii in range(self.n_start): success = False # Generate random starting point if x0 is not provided if x0 is None and self.prob is not None: - x0 = self.prob.sample(1) + x0 = self.prob.sample(1)[0] 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: x_all.append(xopt) y_all.append(yout) - + + if not x_all: + raise RuntimeError("Optimization failed for all initial points — no solution found.") + best_xopt = x_all[np.argmin(np.array(y_all))] return best_xopt # Set the optimization method def set_method(self, method): - self.method = method + self.opt_solver = method - # Set the user options for the Bayesian optimization - def set_options(self, options): - self.options = options + # Set the options for the internal optimization solver + def set_options(self, solver_options): + self.solver_options = solver_options # Method to perform Bayesian optimization - def optimize(self, prob:Problem): - self.problem = prob + def optimize(self): x_train = self.xtrain y_train = self.ytrain @@ -138,15 +171,15 @@ def optimize(self, prob:Problem): self.x_hist = [] self.y_hist = [] - for i in range(self.n_iter): + for i in range(self.bo_maxiter): print(f"*****************************") - print(f"Iteration {i+1}/{self.n_iter}") + print(f"Iteration {i+1}/{self.bo_maxiter}") # Get a new sample point x_new = self._find_best_point(x_train, y_train) # Evaluate the new sample point - y_new = prob.evaluate(np.atleast_2d(x_new)) + y_new = self.prob.evaluate(np.atleast_2d(x_new)) # Update training set x_train = np.vstack([x_train, x_new]) @@ -159,29 +192,48 @@ def optimize(self, prob:Problem): print(f"Sampled point X: {x_new.flatten()}, Observation Y: {y_new.flatten()}") # Save the optimal results and all the training data - self.idx_opt = np.argmin(y_train) - self.x_opt = x_train[self.idx_opt] - self.y_opt = y_train[self.idx_opt] + self.idx_opt = np.argmin(self.y_hist) + self.x_opt = self.y_hist[self.idx_opt] + self.y_opt = self.y_hist[self.idx_opt] self.setTrainingData(x_train, y_train) print() print() - if self.idx_opt < n_init_sample: - print(f"Optimal at initial sample: {self.idx_opt+1}") - else: - print(f"Optimal at BO iteration: {self.idx_opt-n_init_sample+1} ") + print(f"Optimal at BO iteration: {self.idx_opt+1} ") + #if self.idx_opt < n_init_sample: + # print(f"Optimal at initial sample: {self.idx_opt+1}") + #else: + # 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. +def minimizer(fun, x0, method, bounds, constraints, solver_options): + if method != "IPOPT": + if 'grad' in fun: + y = minimize(fun['obj'], x0, method=method, bounds=bounds, jac=fun['grad'], constraints=constraints, options=solver_options) + else: + y = minimize(fun['obj'], x0, method=method, bounds=bounds, constraints=constraints, options=solver_options) + success = y.success + if not success: + print(y.message) + xopt = y.x + yopt = y.fun + else: + ipopt_prob = IpoptProb(fun['obj'], fun['grad'], constraints, bounds, solver_options) + sol, info = ipopt_prob.solve(x0) + + status = info.get('status', -999) + msg = info.get('status_msg', b'unknown error') + if status == 0: + # ipopt returns 0 as success + success = True + else: + warnings.warn(f"Ipopt failed to solve the problem. Status msg: {msg}") + success = False + yopt = info['obj_val'] + xopt = sol -# Find the minimum of the input objective `fun`, using the minimize function from SciPy. -def pyminimize(fun, x0, method, bounds, constraints, options): - y = minimize(fun, x0, method=method, bounds=bounds, constraints=constraints, options=options) - success = y.success - if not success: - print(y.message) - xopt = y.x - yopt = y.fun return xopt, yopt, success diff --git a/src/hiopbbpy/opt/optproblem.py b/src/hiopbbpy/opt/optproblem.py new file mode 100644 index 000000000..0f4573aca --- /dev/null +++ b/src/hiopbbpy/opt/optproblem.py @@ -0,0 +1,79 @@ +import numpy as np +import cyipopt + +""" +Convert a Scipy optimization problem to an Ipopt problem. + +Parameters: + objective (callable): Objective function. + gradient (callable): Gradient of the objective. + constraints_list: Scipy-styple list of dicts with 'type', 'fun', and optional 'jac'. + xbounds (list of tuple): Variable bounds [(x0_lb, x0_ub), ...] + +Returns: + Ipopt-compatible prob and bounds + +Authors: Tucker Hartland + Nai-Yuan Chiang +""" +class IpoptProb: + def __init__(self, objective, gradient, constraints_list, xbounds, solver_options=None): + self.constraints_list = constraints_list + self.eval_f = objective + self.eval_g = gradient + self.xl = [b[0] for b in xbounds] + self.xu = [b[1] for b in xbounds] + self.cl = [] + self.cu = [] + self.nvar = len(xbounds) + self.ncon = len(self.constraints_list) + self.ipopt_options = solver_options + + for con in constraints_list: + if con['type'] == 'eq': + self.cl.append(0.0) + self.cu.append(0.0) + elif con['type'] == 'ineq': + self.cl.append(0.0) + self.cu.append(np.inf) + else: + raise ValueError(f"Unknown constraint type: {con['type']}") + + self.nlp = cyipopt.Problem( + n=self.nvar, + m=self.ncon, + problem_obj=self, + lb=self.xl, + ub=self.xu, + cl=self.cl, + cu=self.cu + ) + + def objective(self, x): + return self.eval_f(x) + + def gradient(self, x): + return self.eval_g(x) + + def constraints(self, x): + return np.array([con['fun'](x) for con in self.constraints_list]) + + def jacobian(self, x): + jacs = [] + for con in self.constraints_list: + if 'jac' in con: + jacs.append(con['jac'](x)) + else: + raise ValueError("Jacobian not provided for constraint.") + return np.vstack(jacs) + + def solve(self, x0, solver_options=None): + ipopt_options = self.ipopt_options + if solver_options is not None: + ipopt_options = solver_options + if ipopt_options is not None: + for key, value in ipopt_options.items(): + self.nlp.add_option(key, value) + + # Solve the optimization problem + return self.nlp.solve(x0) diff --git a/src/hiopbbpy/problems/BraninProblem.py b/src/hiopbbpy/problems/BraninProblem.py index a4485b22e..d45a7575f 100644 --- a/src/hiopbbpy/problems/BraninProblem.py +++ b/src/hiopbbpy/problems/BraninProblem.py @@ -19,11 +19,11 @@ # define the Branin problem class class BraninProblem(Problem): - def __init__(self): + def __init__(self, constraints=None): ndim = 2 xlimits = np.array([[-5.0, 10], [0.0, 15]]) name = 'Branin' - super().__init__(ndim, xlimits, name=name) + super().__init__(ndim, xlimits, name=name, constraints=constraints) def _evaluate(self, x: np.ndarray) -> np.ndarray: @@ -40,17 +40,10 @@ def _evaluate(self, x: np.ndarray) -> np.ndarray: arg1 = (x[:,1] - b * x[:,0]**2 + c * x[:,0] - r) y[:,0] = arg1**2 + s * (1 - t) * np.cos(x[:,0]) + s - ''' - # compute derivatives - dy_dx0 = 2*arg1*(-2*b*x[:,0]+c) - s*(1-t)*np.sin(x[:,0]) - dy_dx1 = 2*arg1 - - dy_dx = np.array([dy_dx0, dy_dx1]) - ''' - return y + diff --git a/src/hiopbbpy/problems/problem.py b/src/hiopbbpy/problems/problem.py index c0895ae33..1d0bae18f 100644 --- a/src/hiopbbpy/problems/problem.py +++ b/src/hiopbbpy/problems/problem.py @@ -10,7 +10,7 @@ # define the general optimization problem class class Problem: - def __init__(self, ndim, xlimits, name=None): + def __init__(self, ndim, xlimits, name=None, constraints=None): self.ndim = ndim self.xlimits = xlimits assert self.xlimits.shape[0] == ndim @@ -19,6 +19,7 @@ def __init__(self, ndim, xlimits, name=None): else: self.name = name self.sampler = qmc.LatinHypercube(ndim) + self.constraints = constraints def _evaluate(self, x: np.ndarray) -> np.ndarray: """ @@ -80,6 +81,9 @@ def sample(self, nsample: int) -> np.ndarray: return xsample + def set_constraints(self, constraints: dict): + self.constraints = constraints + diff --git a/src/hiopbbpy/surrogate_modeling/gp.py b/src/hiopbbpy/surrogate_modeling/gp.py index 9663a75d0..e1cf1955b 100644 --- a/src/hiopbbpy/surrogate_modeling/gp.py +++ b/src/hiopbbpy/surrogate_modeling/gp.py @@ -86,4 +86,36 @@ def train(self, x: np.ndarray, y: np.ndarray) -> np.ndarray: y : ndarray[n, 1] """ - NotImplementedError("Child class of GaussianProcess should implement method train") \ No newline at end of file + NotImplementedError("Child class of GaussianProcess should implement method train") + + # Abstract method for computing the gradient of the mean of the GP at a given input x + def mean_gradient(self, x: np.ndarray) -> np.ndarray: + """ + evaluation of the gradien of GP mean + + Parameters + --------- + x : ndarray[n, nx] + + Returns + ------- + ndarray[n, 1] + Gradient of Mean of GP at x + """ + raise NotImplementedError("Child class of GaussianProcess should implement method mean_gradient") + + # Abstract method for computing the gradient of variance of the GP at a given input x + def variance_gradient(self, x: np.ndarray) -> np.ndarray: + """ + evaluation of the gradient of GP variance + + Parameters + --------- + x: ndarray[n, nx] + + Returns + ------- + ndarray[n, n] + Gradient of variance of GP at w.r.t. x + """ + raise NotImplementedError("Child class of GaussianProcess should implement method variance_gradient") diff --git a/src/hiopbbpy/surrogate_modeling/krg.py b/src/hiopbbpy/surrogate_modeling/krg.py index 8cdb6e787..a9227b7a4 100644 --- a/src/hiopbbpy/surrogate_modeling/krg.py +++ b/src/hiopbbpy/surrogate_modeling/krg.py @@ -32,12 +32,12 @@ def __init__(self, theta, xlimits, ndim, corr="pow_exp", noise0=None, random_sta def mean(self, x): if not self.trained: - ValueError("must train kriging model before utilizing it to predict mean or variances") + raise ValueError("must train kriging model before utilizing it to predict mean or variances") return self.surrogatesmt.predict_values(x) def variance(self, x): if not self.trained: - ValueError("must train kriging model before utilizing it to predict mean or variances") + raise ValueError("must train kriging model before utilizing it to predict mean or variances") return self.surrogatesmt.predict_variances(x) def train(self, x, y): @@ -47,7 +47,16 @@ def train(self, x, y): self.surrogatesmt.train() self.trained = True - - - - + def mean_gradient(self, x: np.ndarray) -> np.ndarray: + if not self.trained: + raise ValueError("must train kriging model before utilizing it to predict gradient") + assert (np.size(x,-1) == self.ndim) + gradient = [ + self.surrogatesmt._predict_derivatives(x, kx) for kx in range(self.ndim) + ] + return np.atleast_2d(gradient).T + + def variance_gradient(self, x: np.ndarray) -> np.ndarray: + if not self.trained: + raise ValueError("must train kriging model before utilizing it to predict gradient") + return self.surrogatesmt.predict_variance_gradient(x)