diff --git a/README.md b/README.md index 2ad941c0..a64c3a77 100644 --- a/README.md +++ b/README.md @@ -4,7 +4,31 @@ Surrogate models and active learning for scientific applications. ## Building -This project uses [pdm](https://pdm-project.org/en/stable/) as its package manager. With pdm installed, run `pdm install` at the root of this repository to install the dependencies. The file [pyproject.toml](pyproject.toml) has the list of dependencies and configurations for the project. Use `pdm build` to build the project. The build artifacts will be in the `dist` directory. Please, find more information about pdm in its website. +### Binaries + +The binaries for the latest version are available at https://github.com/NREL/bbopt/releases/latest. They can be installed through standard installation, e.g., + +using pip (https://pip.pypa.io/en/stable/cli/pip_install/): + +```sh +python -m pip install +https://github.com/NREL/bbopt/archive/refs/tags/v0.4.2.tar.gz +``` + +using conda (https://docs.anaconda.com/working-with-conda/packages/install-packages/): + +```sh +conda install +https://github.com/NREL/bbopt/archive/refs/tags/v0.4.2.tar.gz +``` + +### From source + +This package contains a [pyproject.toml](pyproject.toml) with the list of requirements and dependencies (More about `pyproject.toml` at https://packaging.python.org/en/latest/specifications/pyproject-toml/). This project is configured to use the package manager [pdm](https://pdm-project.org/en/stable/). With pdm installed, run `pdm install` at the root of this repository to install the dependencies. The file [pyproject.toml](pyproject.toml) has the list of dependencies and configurations for the project. Use `pdm build` to build the packages binaries in the `dist` directory. Then follow the steps in the [Binaries section](#binaries) to install. + +### For developers + +Install the package manager [pdm](https://pdm-project.org/en/stable/). To install all the dependencies listted in the [pyproject.toml](pyproject.toml), run `pdm install`. ## Documentation diff --git a/blackboxopt/__init__.py b/blackboxopt/__init__.py index 7c2e8fbc..3ca576d8 100644 --- a/blackboxopt/__init__.py +++ b/blackboxopt/__init__.py @@ -16,7 +16,7 @@ # along with this program. If not, see . __all__ = ["acquisition", "optimize", "sampling", "rbf"] -__version__ = "0.4.1" +__version__ = "0.4.2" from . import acquisition from . import optimize diff --git a/blackboxopt/acquisition.py b/blackboxopt/acquisition.py index 1c67100e..342f200e 100644 --- a/blackboxopt/acquisition.py +++ b/blackboxopt/acquisition.py @@ -31,7 +31,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np @@ -53,7 +53,7 @@ # Local imports from .sampling import NormalSampler, Sampler -from .rbf import RbfModel, RbfType +from .rbf import RbfModel, RbfKernel from .problem import ( ProblemWithConstraint, ProblemNoConstraint, @@ -376,13 +376,13 @@ def acquire( # Evaluate candidates if not listOfSurrogates: samples = surrogateModel.samples() - fx, _ = surrogateModel.eval(x) + fx, _ = surrogateModel(x) else: samples = surrogateModel[0].samples() objdim = len(surrogateModel) fx = np.empty((nCand, objdim)) for i in range(objdim): - fx[:, i], _ = surrogateModel[i].eval(x) + fx[:, i], _ = surrogateModel[i](x) # Create scaled x and scaled distx xlow = np.array([bounds[i][0] for i in range(dim)]) @@ -494,13 +494,13 @@ def acquire( # Evaluate candidates if not listOfSurrogates: samples = surrogateModel.samples() - fx, _ = surrogateModel.eval(x) + fx, _ = surrogateModel(x) else: samples = surrogateModel[0].samples() objdim = len(surrogateModel) fx = np.empty((nCand, objdim)) for i in range(objdim): - fx[:, i], _ = surrogateModel[i].eval(x) + fx[:, i], _ = surrogateModel[i](x) # Create scaled x and scaled distx xlow = np.array([bounds[i][0] for i in range(dim)]) @@ -601,7 +601,7 @@ def acquire( else np.random.choice(self.cycleLength + 2) ) if sample_stage == 0: # InfStep - minimize Mu_n - LDLt = ldl(surrogateModel.get_RBFmatrix()) + LDLt = ldl(surrogateModel.get_RBFmatrix(smoothing=0)) problem = ProblemWithConstraint( lambda x: surrogateModel.mu_measure(x, LDLt=LDLt), lambda x: self.tol @@ -636,7 +636,7 @@ def acquire( assert len(fbounds) == 2 # find min of surrogate model problem = ProblemNoConstraint( - lambda x: surrogateModel.eval(x)[0], + lambda x: surrogateModel(x)[0], bounds, surrogateModel.iindex, ) @@ -657,7 +657,7 @@ def acquire( ) # target for objective function value # use GA method to minimize bumpiness measure - LDLt = ldl(surrogateModel.get_RBFmatrix()) + LDLt = ldl(surrogateModel.get_RBFmatrix(smoothing=0)) problem = ProblemWithConstraint( lambda x: surrogateModel.bumpiness_measure( x, f_target, LDLt @@ -691,7 +691,7 @@ def acquire( assert len(fbounds) == 2 # find the minimum of RBF surface problem = ProblemNoConstraint( - lambda x: surrogateModel.eval(x)[0], + lambda x: surrogateModel(x)[0], bounds, surrogateModel.iindex, ) @@ -723,7 +723,7 @@ def acquire( fbounds[0] ) # target value # use GA method to minimize bumpiness measure - LDLt = ldl(surrogateModel.get_RBFmatrix()) + LDLt = ldl(surrogateModel.get_RBFmatrix(smoothing=0)) problem = ProblemWithConstraint( lambda x: surrogateModel.bumpiness_measure( x, f_target, LDLt @@ -871,9 +871,7 @@ def acquire( ].T # Evaluate the surrogate model on the candidate points and sort them - fcand[iStart:iEnd], _ = surrogateModel.eval( - candidates[iStart:iEnd, :] - ) + fcand[iStart:iEnd], _ = surrogateModel(candidates[iStart:iEnd, :]) ids = np.argsort(fcand[0:iEnd]) remevals -= iEnd - iStart @@ -899,26 +897,18 @@ def acquire( def func_continuous_search(x): x_ = xi.copy() x_[cindex] = x - return surrogateModel.eval(x_)[0] + return surrogateModel(x_)[0] def dfunc_continuous_search(x): x_ = xi.copy() x_[cindex] = x return surrogateModel.jac(x_)[cindex] - # def hessp_continuous_search(x, p): - # x_ = xi.copy() - # x_[cindex] = x - # p_ = np.zeros(dim) - # p_[cindex] = p - # return surrogateModel.hessp(x_, p_)[cindex] - res = minimize( func_continuous_search, xi[cindex], method="L-BFGS-B", jac=dfunc_continuous_search, - # hessp=hessp_continuous_search, bounds=cbounds, options={ "maxfun": remevals, @@ -1025,7 +1015,7 @@ def pareto_front_target(self, paretoFront: np.ndarray) -> np.ndarray: assert objdim > 1 # Create a surrogate model for the Pareto front in the objective space - paretoModel = RbfModel(RbfType.LINEAR) + paretoModel = RbfModel(kernel=RbfKernel.LINEAR) k = np.random.choice(objdim) paretoModel.update_samples( np.array([paretoFront[:, i] for i in range(objdim) if i != k]).T @@ -1046,13 +1036,13 @@ def pareto_front_target(self, paretoFront: np.ndarray) -> np.ndarray: ) def delta_f(tau): - tauk, _ = paretoModel.eval(tau) + tauk, _ = paretoModel(tau) _tau = np.concatenate((tau[0:k], tauk, tau[k:])) return -tree.query(_tau)[0] # Minimize delta_f res = differential_evolution(delta_f, boundsPareto) - tauk, _ = paretoModel.eval(res.x) + tauk, _ = paretoModel(res.x) tau = np.concatenate((res.x[0:k], tauk, res.x[k:])) return tau @@ -1175,7 +1165,7 @@ def acquire( endpoints = np.empty((objdim, dim)) for i in range(objdim): minimumPointProblem = ProblemNoConstraint( - lambda x: surrogateModels[i].eval(x)[0], bounds, iindex + lambda x: surrogateModels[i](x)[0], bounds, iindex ) res = pymoo_minimize( minimumPointProblem, @@ -1421,9 +1411,7 @@ def acquire( fnondominatedAndBestCandidates = np.concatenate( ( paretoFront, - np.array( - [s.eval(bestCandidates)[0] for s in surrogateModels] - ).T, + np.array([s(bestCandidates)[0] for s in surrogateModels]).T, ), axis=0, ) @@ -1508,7 +1496,7 @@ def acquire( cheapProblem = ProblemWithConstraint( self.fun, lambda x: np.transpose( - [surrogateModels[i].eval(x)[0] for i in range(gdim)] + [surrogateModels[i](x)[0] for i in range(gdim)] ), bounds, iindex, diff --git a/blackboxopt/optimize.py b/blackboxopt/optimize.py index 10872511..b71aff0a 100644 --- a/blackboxopt/optimize.py +++ b/blackboxopt/optimize.py @@ -31,7 +31,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from typing import Callable, Optional, Union @@ -781,9 +781,9 @@ def multistart_stochastic_response_surface( if out_local.fx < out.fx: out.x[:] = out_local.x out.fx = out_local.fx - out.samples[ - out.nfev : out.nfev + out_local.nfev, : - ] = out_local.samples + out.samples[out.nfev : out.nfev + out_local.nfev, :] = ( + out_local.samples + ) out.fsamples[out.nfev : out.nfev + out_local.nfev] = out_local.fsamples out.nfev = out.nfev + out_local.nfev @@ -898,7 +898,11 @@ def target_value_optimization( # max value of f if surrogateModel.nsamples() - out.nfev > 0: - maxf = np.max(surrogateModel.get_fsamples()[: -out.nfev]).item() + maxf = np.max( + surrogateModel.get_fsamples()[ + 0 : surrogateModel.nsamples() - out.nfev + ] + ).item() else: maxf = -np.Inf if out.nfev > 0: @@ -1681,7 +1685,7 @@ def gosac( # Evaluate the surrogate at the best candidates sCandidates = np.empty((len(bestCandidates), gdim)) for i in range(gdim): - sCandidates[:, i], _ = surrogateModels[i].eval(bestCandidates) + sCandidates[:, i], _ = surrogateModels[i](bestCandidates) # Find the minimum number of constraint violations constraintViolation = [ diff --git a/blackboxopt/problem.py b/blackboxopt/problem.py index 1568073a..ad2d4f5d 100644 --- a/blackboxopt/problem.py +++ b/blackboxopt/problem.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np @@ -161,7 +161,7 @@ def _evaluate(self, X, out): out["F"] = np.empty((x.shape[0], self.n_obj)) for i in range(self.n_obj): out["F"][:, i] = np.absolute( - self.surrogateModels[i].eval(x)[0] - self.tau[i] + self.surrogateModels[i](x)[0] - self.tau[i] ) @@ -185,4 +185,4 @@ def _evaluate(self, X, out): x = _dict_to_array(X) out["F"] = np.empty((x.shape[0], self.n_obj)) for i in range(self.n_obj): - out["F"][:, i] = self.surrogateModels[i].eval(x)[0] + out["F"][:, i] = self.surrogateModels[i](x)[0] diff --git a/blackboxopt/rbf.py b/blackboxopt/rbf.py index 3588a933..9057145a 100644 --- a/blackboxopt/rbf.py +++ b/blackboxopt/rbf.py @@ -20,22 +20,46 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from typing import Optional import warnings import numpy as np -from enum import Enum # Scipy imports from scipy.spatial.distance import cdist from scipy.linalg import solve, solve_triangular +from scipy.special import comb # Local imports from .sampling import Sampler +from .rbf_kernel import RbfKernel, KERNEL_DERIVATIVE_OVER_R_FUNC, KERNEL_FUNC -RbfType = Enum("RbfType", ["LINEAR", "CUBIC", "THINPLATE"]) + +def _order2_monomials(x: np.ndarray) -> np.ndarray: + m = x.shape[0] + dim = x.shape[1] + out = np.zeros((m, (dim * (dim + 1)) // 2)) + count = 0 + for i in range(dim): + for j in range(i, dim): + out[:, count] = x[:, i] * x[:, j] + count += 1 + return out + + +def _d_order2_monomials(x: np.ndarray) -> np.ndarray: + dim = len(x) + assert x.ndim == 1 + out = np.zeros((dim, (dim * (dim + 1)) // 2)) + count = 0 + for i in range(dim): + for j in range(i, dim): + out[i, count] += x[j] + out[j, count] += x[i] + count += 1 + return out class RbfFilter: @@ -95,12 +119,10 @@ class RbfModel: Attributes ---------- - type : RbfType, optional - Defines the function phi used in the RBF model. The options are: - - - RbfType.LINEAR: phi(r) = r. - - RbfType.CUBIC: phi(r) = r^3. - - RbfType.THINPLATE: phi(r) = r^2 * log(r). + smoothing : float, optional + Smoothing parameter. The interpolant perfectly fits the data when this + is set to 0. For large values, the interpolant approaches a least + squares fit of a polynomial with the specified degree. Default is 0. iindex : tuple, optional Indices of the input space that are integer. The default is (). filter : RbfFilter, optional @@ -110,14 +132,54 @@ class RbfModel: def __init__( self, - rbf_type: RbfType = RbfType.CUBIC, + *, + smoothing: float = 0.0, + kernel: RbfKernel = RbfKernel.CUBIC, + epsilon: float = 1.0, iindex: tuple[int, ...] = (), filter: Optional[RbfFilter] = None, ): - self.type = rbf_type + """Initialize the RBF model + + By default, the model uses a cubic kernel with no smoothing. + + Parameters + ---------- + smoothing : float, optional + Smoothing parameter. The interpolant perfectly fits the data when this + is set to 0. For large values, the interpolant approaches a least + squares fit of a polynomial with the specified degree. Default is 0. + kernel : RbfKernel + Defines the function phi used in the RBF model. The options are listed + in the RbfKernel enum. + epsilon : float, optional + Shape parameter that scales the input to the RBF. If `kernel` is + 'linear', 'thin_plate_spline', 'cubic', or 'quintic', this defaults to + 1 and can be ignored because it has the same effect as scaling the + smoothing parameter. Defaults to 1. + iindex : tuple, optional + Indices of the input space that are integer. The default is (). + filter : RbfFilter, optional + Filter used with the function values. The default is RbfFilter() which + is the identity function. + """ + + self.smoothing = smoothing self.iindex = iindex self.filter = RbfFilter() if filter is None else filter + # Set kernel and the degree of the polynomial tail + self._kernel = kernel + if kernel in (RbfKernel.LINEAR, RbfKernel.MULTIQUADRIC): + self._degree = 0 + elif kernel in (RbfKernel.CUBIC, RbfKernel.THINPLATE): + self._degree = 1 + elif kernel == RbfKernel.QUINTIC: + self._degree = 2 + else: + self._degree = None + self._eps = epsilon + self._valid_coefficients = True self._m = 0 self._x = np.array([]) @@ -126,6 +188,12 @@ def __init__( self._PHI = np.array([]) self._P = np.array([]) + # Scaling for better condition number + self._sx = np.array([]) + self._scale = np.array([]) + self._avg = np.array([]) + self._change_scale_factor = 2.0 # Change the scale factor when new scale is 2 times larger than the current one + def reserve(self, maxeval: int, dim: int) -> None: """Reserve space for the RBF model. @@ -144,11 +212,15 @@ def reserve(self, maxeval: int, dim: int) -> None: if self._x.size == 0: self._x = np.empty((maxeval, dim)) + self._sx = np.empty((maxeval, dim)) else: additional_rows = max(0, maxeval - self._x.shape[0]) self._x = np.concatenate( (self._x, np.empty((additional_rows, dim))), axis=0 ) + self._sx = np.concatenate( + (self._sx, np.empty((additional_rows, dim))), axis=0 + ) if self._fx.size == 0: self._fx = np.empty(maxeval) @@ -221,137 +293,11 @@ def pdim(self) -> int: Dimension of the polynomial tail. """ dim = self.dim() - if self.type == RbfType.LINEAR: - return 1 - elif self.type in (RbfType.CUBIC, RbfType.THINPLATE): - return 1 + dim - else: - raise ValueError("Unknown RBF type") - - def phi(self, r): - """Applies the function phi to the distance(s) r. - Parameters - ---------- - r : array_like - Distance(s) between points. - - Returns - ------- - out: array_like - Phi-value of the distances provided on input. - """ - if self.type == RbfType.LINEAR: - return r - elif self.type == RbfType.CUBIC: - return np.power(r, 3) - elif self.type == RbfType.THINPLATE: - if not hasattr(r, "__len__"): - if r > 0: - return r**2 * np.log(r) - else: - return 0 - else: - ret = np.zeros_like(r) - ret[r > 0] = np.multiply( - np.power(r[r > 0], 2), np.log(r[r > 0]) - ) - return ret + if self._degree is not None: + return int(comb(dim + self._degree, dim, exact=True)) else: - raise ValueError("Unknown RBF type") - - def dphi(self, r): - """Derivative of the function phi at the distance(s) r. - - Parameters - ---------- - r : array_like - Distance(s) between points. - - Returns - ------- - out: array_like - Derivative of the phi-value of the distances provided on input. - """ - if self.type == RbfType.LINEAR: - return np.ones(r.shape) - elif self.type == RbfType.CUBIC: - return 3 * np.power(r, 2) - elif self.type == RbfType.THINPLATE: - if not hasattr(r, "__len__"): - if r > 0: - return 2 * r * np.log(r) + r - else: - return 0 - else: - ret = np.zeros_like(r) - ret[r > 0] = ( - 2 * np.multiply(r[r > 0], np.log(r[r > 0])) + r[r > 0] - ) - return ret - else: - raise ValueError("Unknown RBF type") - - def dphiOverR(self, r): - """Derivative of the function phi divided by r at the distance(s) r. - - Parameters - ---------- - r : array_like - Distance(s) between points. - - Returns - ------- - out: array_like - Derivative of the phi-value of the distances provided on input - divided by the distance. - """ - if self.type == RbfType.LINEAR: - return np.ones(r.shape) / r - elif self.type == RbfType.CUBIC: - return 3 * r - elif self.type == RbfType.THINPLATE: - if not hasattr(r, "__len__"): - if r > 0: - return 2 * np.log(r) + 1 - else: - return 0 - else: - ret = np.zeros_like(r) - ret[r > 0] = 2 * np.log(r[r > 0]) + 1 - return ret - else: - raise ValueError("Unknown RBF type") - - def ddphi(self, r): - """Second derivative of the function phi at the distance(s) r. - - Parameters - ---------- - r : array_like - Distance(s) between points. - - Returns - ------- - out: array_like - Second derivative of the phi-value of the distances provided on input. - """ - if self.type == RbfType.LINEAR: - return np.zeros(r.shape) - elif self.type == RbfType.CUBIC: - return 6 * r - elif self.type == RbfType.THINPLATE: - if not hasattr(r, "__len__"): - if r > 0: - return 2 * np.log(r) + 3 - else: - return 0 - else: - ret = np.zeros_like(r) - ret[r > 0] = 2 * np.log(r[r > 0]) + 3 - return ret - else: - raise ValueError("Unknown RBF type") + return 0 def pbasis(self, x: np.ndarray) -> np.ndarray: """Computes the polynomial tail matrix for a given set of points. @@ -368,14 +314,20 @@ def pbasis(self, x: np.ndarray) -> np.ndarray: """ dim = self.dim() m = x.size // dim + assert self._degree is not None # Set up the polynomial tail matrix P - if self.type == RbfType.LINEAR: - return np.ones((m, 1)) - elif self.type in (RbfType.CUBIC, RbfType.THINPLATE): - return np.concatenate((np.ones((m, 1)), x.reshape(m, -1)), axis=1) - else: - raise ValueError("Invalid polynomial tail") + out = np.ones((m, 1)) + if self._degree >= 1: + out = np.concatenate((out, x.reshape(m, -1)), axis=1) + if self._degree >= 2: + out = np.concatenate( + (out, _order2_monomials(x.reshape(m, -1))), axis=1 + ) + if self._degree >= 3: + raise ValueError("Higher order polynomials are not supported") + + return out def dpbasis(self, x: np.ndarray) -> np.ndarray: """Computes the derivative of the polynomial tail matrix for a given x. @@ -391,41 +343,19 @@ def dpbasis(self, x: np.ndarray) -> np.ndarray: Derivative of the polynomial tail matrix for the input point. """ dim = self.dim() + assert self._degree is not None - if self.type == RbfType.LINEAR: - return np.zeros((1, 1)) - elif self.type in (RbfType.CUBIC, RbfType.THINPLATE): - return np.concatenate((np.zeros((1, dim)), np.eye(dim)), axis=0) - else: - raise ValueError("Invalid polynomial tail") - - def ddpbasis(self, x: np.ndarray, p: np.ndarray) -> np.ndarray: - """Computes the second derivative of the polynomial tail matrix for a - given x and direction p. - - Parameters - ---------- - x : numpy.ndarray - Point in a d-dimensional space. - p : numpy.ndarray - Direction in which the second derivative is evaluated. - - Returns - ------- - out: numpy.ndarray - Second derivative of the polynomial tail matrix for the input point - and direction. - """ - dim = self.dim() + out = np.zeros((dim, 1)) + if self._degree >= 1: + out = np.concatenate((out, np.eye(dim)), axis=1) + if self._degree >= 2: + out = np.concatenate((out, _d_order2_monomials(x).T), axis=1) + if self._degree >= 3: + raise ValueError("Higher order polynomials are not supported") - if self.type == RbfType.LINEAR: - return np.zeros((1, 1)) - elif self.type in (RbfType.CUBIC, RbfType.THINPLATE): - return np.zeros((dim + 1, dim)) - else: - raise ValueError("Invalid polynomial tail") + return out - def eval(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + def __call__(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray]: """Evaluates the model at one or multiple points. Parameters @@ -447,12 +377,17 @@ def eval(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray]: ) dim = self.dim() + phi = KERNEL_FUNC[self._kernel] + + # Scale x and samples + xscaled = (x.reshape(-1, dim) - self._avg) / self._scale + sscaled = self._sx[0 : self._m, :] # compute pairwise distances between candidates and sampled points - D = cdist(x.reshape(-1, dim), self.samples()) + D = cdist(xscaled, sscaled) * self._eps - Px = self.pbasis(x) - y = np.matmul(self.phi(D), self._coef[0 : self._m]) + np.dot( + Px = self.pbasis(xscaled) + y = np.matmul(phi(D), self._coef[0 : self._m]) + np.dot( Px, self._coef[self._m : self._m + Px.shape[1]] ) @@ -482,65 +417,25 @@ def jac(self, x: np.ndarray) -> np.ndarray: ) dim = self.dim() + dphiOverR = KERNEL_DERIVATIVE_OVER_R_FUNC[self._kernel] - # compute pairwise distances between candidates and sampled points - d = cdist(x.reshape(-1, dim), self.samples()).flatten() - - A = np.array([self.dphiOverR(d[i]) * x for i in range(d.size)]) - B = self.dpbasis(x) - - y = np.matmul(A.T, self._coef[0 : self._m]) + np.matmul( - B.T, self._coef[self._m : self._m + B.shape[0]] - ) - - return y.flatten() - - def hessp(self, x: np.ndarray, p: np.ndarray) -> np.ndarray: - r"""Evaluates the Hessian of the model at x in the direction of p. - - .. math:: - - H(f)(x) v = \sum_{i=1}^{m} \beta_i \left( - \phi''(\|x - x_i\|)\frac{(x^Tv)x}{\|x - x_i\|^2} + - \frac{\phi'(\|x - x_i\|)}{\|x - x_i\|} - \left(v - \frac{(x^Tv)x}{\|x - x_i\|^2}\right) - \right) - + \sum_{i=1}^{n} \beta_{m+i} H(p_i)(x) v. - - Parameters - ---------- - x : np.ndarray - Point in a d-dimensional space. - p : np.ndarray - Direction in which the Hessian is evaluated. - - Returns - ------- - numpy.ndarray - Value for the Hessian of the RBF model at x in the direction of p. - """ - if self._valid_coefficients is False: - raise RuntimeError( - "Invalid coefficients. Run update_coefficients() before evaluating the model." - ) - - dim = self.dim() + # Scale x and samples + xscaled = (x.reshape(-1, dim) - self._avg) / self._scale + sscaled = self._sx[0 : self._m, :] # compute pairwise distances between candidates and sampled points - d = cdist(x.reshape(-1, dim), self.samples()).flatten() + d = cdist(xscaled, sscaled).flatten() - xxTp = np.dot(p, x) * x - A = np.array( - [ - self.ddphi(d[i]) * (xxTp / (d[i] * d[i])) - + self.dphiOverR(d[i]) * (p - (xxTp / (d[i] * d[i]))) - for i in range(d.size) - ] + A = np.matmul( + np.array( + [dphiOverR(d[i] * self._eps) * xscaled for i in range(d.size)] + ), + np.diag(self._eps / self._scale), ) - B = self.ddpbasis(x, p) + B = np.matmul(np.diag(1 / self._scale), self.dpbasis(xscaled)) y = np.matmul(A.T, self._coef[0 : self._m]) + np.matmul( - B.T, self._coef[self._m : self._m + B.shape[0]] + B, self._coef[self._m : self._m + B.shape[1]] ) return y.flatten() @@ -606,30 +501,73 @@ def update_samples(self, xNew: np.ndarray, distNew=None) -> None: newm = xNew.shape[0] dim = xNew.shape[1] m = oldm + newm + phi = KERNEL_FUNC[self._kernel] if oldm > 0: assert dim == self.dim() if newm == 0: return - # Compute distances between new points and sampled points - if distNew is None: - if oldm == 0: - distNew = cdist(xNew, xNew) - else: - distNew = cdist( - xNew, np.concatenate((self.samples(), xNew), axis=0) - ) + # Update x + self._x[oldm:m, :] = xNew - self.reserve(m, dim) + # Compute new scaling factor + if len(self._scale) > 0: + xmax = np.max( + np.concatenate((xNew, self._x[0:oldm, :]), axis=0), axis=0 + ) + xmin = np.min( + np.concatenate((xNew, self._x[0:oldm, :]), axis=0), axis=0 + ) + new_scale = (xmax - xmin) / 2 + new_scale = np.where(new_scale == 0, 1, new_scale) + else: + xmax = np.max(xNew, axis=0) + xmin = np.min(xNew, axis=0) + new_scale = (xmax - xmin) / 2 + new_scale = np.where(new_scale == 0, 1, new_scale) + self._scale = new_scale + self._avg = (xmax + xmin) / 2 + + if len(self._scale) == 0 or np.all( + new_scale < self._scale * self._change_scale_factor + ): + # Scale points + xscaled = self._sx[oldm:m, :] + xscaled[:] = (xNew - self._avg) / self._scale + + # Compute distances between new points and sampled points + if distNew is None: + if oldm == 0: + distNew = cdist(xscaled, xscaled) + else: + sscaled = (self._x[0:oldm, :] - self._avg) / self._scale + distNew = cdist( + xscaled, + np.concatenate((sscaled, xscaled), axis=0), + ) - # Update matrices _PHI and _P - self._PHI[oldm:m, 0:m] = self.phi(distNew) - self._PHI[0:oldm, oldm:m] = self._PHI[oldm:m, 0:oldm].T - self._P[oldm:m, :] = self.pbasis(xNew) + self.reserve(m, dim) - # Update x - self._x[oldm:m, :] = xNew + # Update matrices _PHI and _P + self._PHI[oldm:m, 0:m] = phi(distNew * self._eps) + self._PHI[0:oldm, oldm:m] = self._PHI[oldm:m, 0:oldm].T + self._P[oldm:m, :] = self.pbasis(xscaled) + else: + # Update scaling factor + self._scale = new_scale + self._avg = (xmax + xmin) / 2 + + # Scale points + xscaled = self._sx[0:m, :] + xscaled[:] = (self._x[0:m, :] - self._avg) / self._scale + + # Recompute distances between sampled points + distNew = cdist(xscaled, xscaled) + + # Update matrices _PHI and _P + self._PHI[0:m, 0:m] = phi(distNew * self._eps) + self._P[0:m, :] = self.pbasis(xscaled) # Update m self._m = m @@ -680,12 +618,21 @@ def create_initial_design( if count > 100: raise RuntimeError("Cannot create valid initial design") + # Compute scaling factor + xmax = np.max(self.samples(), axis=0) + xmin = np.min(self.samples(), axis=0) + self._scale = (xmax - xmin) / 2 + self._scale = np.where(self._scale == 0, 1, self._scale) + self._avg = (xmax + xmin) / 2 + # Compute distances between new points and sampled points - distNew = cdist(self.samples(), self.samples()) + self._sx[0:m, :] = (self.samples() - self._avg) / self._scale + distNew = cdist(self._sx[0:m, :], self._sx[0:m, :]) - # Set matrix _PHI - self._PHI[0:m, 0:m] = self.phi(distNew) - self._PHI[0:0, 0:m] = self._PHI[0:m, 0:0].T + # Set matrices _PHI and _P + phi = KERNEL_FUNC[self._kernel] + self._PHI[0:m, 0:m] = phi(distNew * self._eps) + self._P[0:m, :] = self.pbasis(self._sx[0:m, :]) # Coefficients are not valid self._valid_coefficients = False @@ -734,18 +681,32 @@ def get_matrixP(self) -> np.ndarray: """ return self._P[0 : self._m, :] - def get_RBFmatrix(self) -> np.ndarray: + def get_RBFmatrix( + self, *, smoothing: Optional[float] = None + ) -> np.ndarray: """Get the matrix used to compute the RBF weights. + Parameters + ---------- + smoothing : float, optional + Smoothing parameter, if different from the one used in the model. + Returns ------- out: np.ndarray (m+pdim)-by-(m+pdim) matrix used to compute the RBF weights. """ + if smoothing is None: + smoothing = self.smoothing + pdim = self.pdim() return np.block( [ - [self._PHI[0 : self._m, 0 : self._m], self.get_matrixP()], + [ + self._PHI[0 : self._m, 0 : self._m] + + smoothing * np.eye(self._m), + self.get_matrixP(), + ], [self.get_matrixP().T, np.zeros((pdim, pdim))], ] ) @@ -791,11 +752,18 @@ def mu_measure(self, x: np.ndarray, xdist=None, LDLt=()) -> float: Optimization. Journal of Global Optimization 19, 201–227 (2001). https://doi.org/10.1023/A:1011255519438 """ + phi = KERNEL_FUNC[self._kernel] + # compute rbf value of the new point x + xscaled = (x - self._avg) / self._scale + sscaled = self._sx[0 : self._m, :] if xdist is None: - xdist = cdist(x.reshape(1, -1), self.samples()) + xdist = cdist(xscaled.reshape(1, -1), sscaled) newRow = np.concatenate( - (np.asarray(self.phi(xdist)).flatten(), self.pbasis(x).flatten()) + ( + np.asarray(phi(xdist * self._eps)).flatten(), + self.pbasis(xscaled).flatten(), + ) ) if LDLt: @@ -833,15 +801,15 @@ def mu_measure(self, x: np.ndarray, xdist=None, LDLt=()) -> float: l01[i] /= d0[i, i] # 3. d = \phi(0) - l_{01}^T D_0 l_{01} and \mu = 1/d - d = self.phi(0) - np.dot(l01, D0l01) + d = phi(0) - np.dot(l01, D0l01) mu = 1 / d if d != 0 else np.inf if not LDLt or mu == np.inf: # set up matrices for solving the linear system A_aug = np.block( [ - [self.get_RBFmatrix(), newRow.reshape(-1, 1)], - [newRow, self.phi(0)], + [self.get_RBFmatrix(smoothing=0), newRow.reshape(-1, 1)], + [newRow, phi(0)], ] ) @@ -857,16 +825,7 @@ def mu_measure(self, x: np.ndarray, xdist=None, LDLt=()) -> float: # Return huge value, only occurs if the matrix is ill-conditioned mu = np.inf - # Order of the polynomial tail - if self.type == RbfType.LINEAR: - m0 = 0 - elif self.type in (RbfType.CUBIC, RbfType.THINPLATE): - m0 = 1 - else: - raise ValueError("Unknown RBF type") - # Get the absolute value of mu - mu *= (-1) ** (m0 + 1) if mu < 0: # Return huge value, only occurs if the matrix is ill-conditioned return np.inf @@ -898,16 +857,16 @@ def bumpiness_measure(self, x: np.ndarray, target, LDLt=()) -> float: Optimization. Journal of Global Optimization 19, 201–227 (2001). https://doi.org/10.1023/A:1011255519438 """ - absmu = self.mu_measure(x, LDLt=LDLt) + mu = self.mu_measure(x, LDLt=LDLt) assert ( - absmu > 0 + mu > 0 ) # if absmu == 0, the linear system in the surrogate model singular - if absmu == np.inf: + if mu == np.inf: # Return huge value, only occurs if the matrix is ill-conditioned return np.inf # predict RBF value of x - yhat, _ = self.eval(x) + yhat, _ = self(x) assert yhat.size == 1 # sanity check # Compute the distance between the predicted value and the target @@ -916,5 +875,5 @@ def bumpiness_measure(self, x: np.ndarray, target, LDLt=()) -> float: # dist = tol # use sqrt(gn) as the bumpiness measure to avoid underflow - sqrtgn = np.sqrt(absmu) * dist + sqrtgn = np.sqrt(mu) * dist return sqrtgn diff --git a/blackboxopt/rbf_kernel.py b/blackboxopt/rbf_kernel.py new file mode 100644 index 00000000..b0ac5bab --- /dev/null +++ b/blackboxopt/rbf_kernel.py @@ -0,0 +1,211 @@ +from enum import Enum +from types import MappingProxyType +import numpy as np + +RbfKernel = Enum( + "RbfKernel", + [ + "LINEAR", + "CUBIC", + "THINPLATE", + "QUINTIC", + "MULTIQUADRIC", + "INVERSE_MULTIQUADRIC", + "INVERSE_QUADRATIC", + "GAUSSIAN", + ], +) + +# Kernel functions + + +def linear(r): + return np.negative(r) + + +def thin_plate_spline(r): + if not hasattr(r, "__len__"): + ret = 1 if r == 0 else r + else: + ret = np.where(np.asarray(r) == 0, 1, r) + return ret**2 * np.log(ret) + + +def cubic(r): + return np.power(r, 3) + + +def quintic(r): + return -np.power(r, 5) + + +def multiquadric(r): + return -np.sqrt(np.power(r, 2) + 1) + + +def inverse_multiquadric(r): + return 1 / np.sqrt(np.power(r, 2) + 1) + + +def inverse_quadratic(r): + return 1 / (np.power(r, 2) + 1) + + +def gaussian(r): + return np.exp(-np.power(r, 2)) + + +# Kernel derivatives + + +def d_linear(r: np.ndarray): + return -np.ones(r.shape) + + +def d_thin_plate_spline(r: np.ndarray): + ret = np.where(r == 0, 1, r) + return 2 * ret * np.log(ret) + ret + + +def d_cubic(r: np.ndarray): + return 3 * r**2 + + +def d_quintic(r: np.ndarray): + return -5 * r**4 + + +def d_multiquadric(r: np.ndarray): + return -r / np.sqrt(r**2 + 1) + + +def d_inverse_multiquadric(r: np.ndarray): + return -r / np.sqrt(r**2 + 1) ** 3 + + +def d_inverse_quadratic(r: np.ndarray): + return -2 * r / (r**2 + 1) ** 2 + + +def d_gaussian(r: np.ndarray): + return -2 * r * np.exp(-(r**2)) + + +# Kernel derivatives over r + + +def d_linear_over_r(r: np.ndarray): + return -1 / r + + +def d_thin_plate_spline_over_r(r: np.ndarray): + return 2 * np.log(r) + 1 + + +def d_cubic_over_r(r: np.ndarray): + return 3 * r + + +def d_quintic_over_r(r: np.ndarray): + return -5 * r**3 + + +def d_multiquadric_over_r(r: np.ndarray): + return -1 / np.sqrt(r**2 + 1) + + +def d_inverse_multiquadric_over_r(r: np.ndarray): + return -1 / np.sqrt(r**2 + 1) ** 3 + + +def d_inverse_quadratic_over_r(r: np.ndarray): + return -2 / (r**2 + 1) ** 2 + + +# Kernel second derivatives + + +def d2_linear(r: np.ndarray): + return np.zeros(r.shape) + + +def d2_thin_plate_spline(r: np.ndarray): + ret = np.where(r == 0, 1, r) + return 2 * np.log(ret) + 3 + + +def d2_cubic(r: np.ndarray): + return 6 * r + + +def d2_quintic(r: np.ndarray): + return -20 * r**3 + + +def d2_multiquadric(r: np.ndarray): + return -1 / np.sqrt(r**2 + 1) ** 3 + + +def d2_inverse_multiquadric(r: np.ndarray): + return (2 * r**2 - 1) / np.sqrt(r**2 + 1) ** 5 + + +def d2_inverse_quadratic(r: np.ndarray): + return (6 * r**2 - 2) / (r**2 + 1) ** 3 + + +def d2_gaussian(r: np.ndarray): + return (4 * r**2 - 2) * np.exp(-(r**2)) + + +KERNEL_FUNC = MappingProxyType( + { + RbfKernel.LINEAR: linear, + RbfKernel.CUBIC: cubic, + RbfKernel.THINPLATE: thin_plate_spline, + RbfKernel.QUINTIC: quintic, + RbfKernel.MULTIQUADRIC: multiquadric, + RbfKernel.INVERSE_MULTIQUADRIC: inverse_multiquadric, + RbfKernel.INVERSE_QUADRATIC: inverse_quadratic, + RbfKernel.GAUSSIAN: gaussian, + } +) + +KERNEL_DERIVATIVE_FUNC = MappingProxyType( + { + RbfKernel.LINEAR: d_linear, + RbfKernel.CUBIC: d_cubic, + RbfKernel.THINPLATE: d_thin_plate_spline, + RbfKernel.QUINTIC: d_quintic, + RbfKernel.MULTIQUADRIC: d_multiquadric, + RbfKernel.INVERSE_MULTIQUADRIC: d_inverse_multiquadric, + RbfKernel.INVERSE_QUADRATIC: d_inverse_quadratic, + RbfKernel.GAUSSIAN: d_gaussian, + } +) + +KERNEL_DERIVATIVE_OVER_R_FUNC = MappingProxyType( + { + RbfKernel.LINEAR: d_linear_over_r, + RbfKernel.CUBIC: d_cubic_over_r, + RbfKernel.THINPLATE: d_thin_plate_spline_over_r, + RbfKernel.QUINTIC: d_quintic_over_r, + RbfKernel.MULTIQUADRIC: d_multiquadric_over_r, + RbfKernel.INVERSE_MULTIQUADRIC: d_inverse_multiquadric_over_r, + RbfKernel.INVERSE_QUADRATIC: d_inverse_quadratic_over_r, + RbfKernel.GAUSSIAN: d_gaussian, + } +) + +KERNEL_SECOND_DERIVATIVE_FUNC = MappingProxyType( + { + RbfKernel.LINEAR: d2_linear, + RbfKernel.CUBIC: d2_cubic, + RbfKernel.THINPLATE: d2_thin_plate_spline, + RbfKernel.QUINTIC: d2_quintic, + RbfKernel.MULTIQUADRIC: d2_multiquadric, + RbfKernel.INVERSE_MULTIQUADRIC: d2_inverse_multiquadric, + RbfKernel.INVERSE_QUADRATIC: d2_inverse_quadratic, + RbfKernel.GAUSSIAN: d2_gaussian, + } +) diff --git a/blackboxopt/sampling.py b/blackboxopt/sampling.py index a99620b9..5ce3b102 100644 --- a/blackboxopt/sampling.py +++ b/blackboxopt/sampling.py @@ -31,7 +31,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/docs/conf.py b/docs/conf.py index 01b8b128..14f77fa2 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -17,7 +17,7 @@ project = "Black-box Opt" copyright = "2024, Alliance for Sustainable Energy, LLC" author = "Weslley S. Pereira" -release = "0.4.1" +release = "0.4.2" # -- General configuration --------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration diff --git a/examples/multiobj/pareto_front.ipynb b/examples/multiobj/pareto_front.ipynb index 56159921..42b0e873 100644 --- a/examples/multiobj/pareto_front.ipynb +++ b/examples/multiobj/pareto_front.ipynb @@ -29,7 +29,7 @@ "__maintainer__ = \"Weslley S. Pereira\"\n", "__email__ = \"weslley.dasilvapereira@nrel.gov\"\n", "__credits__ = [\"Weslley S. Pereira\"]\n", - "__version__ = \"0.4.1\"\n", + "__version__ = \"0.4.2\"\n", "__deprecated__ = False" ] }, diff --git a/examples/multiobj/socemo.ipynb b/examples/multiobj/socemo.ipynb index 46ebd721..be5ffbc0 100644 --- a/examples/multiobj/socemo.ipynb +++ b/examples/multiobj/socemo.ipynb @@ -28,7 +28,7 @@ "__maintainer__ = \"Weslley S. Pereira\"\n", "__email__ = \"weslley.dasilvapereira@nrel.gov\"\n", "__credits__ = [\"Weslley S. Pereira\"]\n", - "__version__ = \"0.4.1\"\n", + "__version__ = \"0.4.2\"\n", "__deprecated__ = False" ] }, @@ -580,8 +580,8 @@ " s[j].update_coefficients(res.fsamples[:i1,j])\n", "\n", " x = np.arange(-100, 100, 0.1).reshape(-1, 1)\n", - " y1 = s[0].eval(x)[0]\n", - " y2 = s[1].eval(x)[0]\n", + " y1 = s[0](x)[0]\n", + " y2 = s[1](x)[0]\n", " fx = objf(x)\n", " f1x = [f[0] for f in fx]\n", " f2x = [f[1] for f in fx]\n", diff --git a/examples/opt_with_constr/example_gosac_1.py b/examples/opt_with_constr/example_gosac_1.py index 557f3875..b5d6e5bf 100644 --- a/examples/opt_with_constr/example_gosac_1.py +++ b/examples/opt_with_constr/example_gosac_1.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Juliane Mueller", "Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from blackboxopt.optimize import gosac diff --git a/examples/opt_with_constr/example_gosac_2.py b/examples/opt_with_constr/example_gosac_2.py index 34f32e8f..e9ff34f7 100644 --- a/examples/opt_with_constr/example_gosac_2.py +++ b/examples/opt_with_constr/example_gosac_2.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Juliane Mueller", "Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from blackboxopt.optimize import gosac diff --git a/examples/sampling.ipynb b/examples/sampling.ipynb index 413766af..1a857fd2 100644 --- a/examples/sampling.ipynb +++ b/examples/sampling.ipynb @@ -28,7 +28,7 @@ "__maintainer__ = \"Weslley S. Pereira\"\n", "__email__ = \"weslley.dasilvapereira@nrel.gov\"\n", "__credits__ = [\"Weslley S. Pereira\"]\n", - "__version__ = \"0.4.1\"\n", + "__version__ = \"0.4.2\"\n", "__deprecated__ = False" ] }, diff --git a/examples/single_obj_rbf/LocalStochRBFstop.py b/examples/single_obj_rbf/LocalStochRBFstop.py index aa8edfaf..2628ded4 100644 --- a/examples/single_obj_rbf/LocalStochRBFstop.py +++ b/examples/single_obj_rbf/LocalStochRBFstop.py @@ -26,11 +26,11 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from optprogram1 import read_check_data_file -from blackboxopt.rbf import RbfType, RbfModel, MedianLpfFilter +from blackboxopt.rbf import RbfKernel, RbfModel, MedianLpfFilter from blackboxopt.optimize import stochastic_response_surface from blackboxopt.sampling import NormalSampler, Sampler, SamplingStrategy from blackboxopt.acquisition import CoordinatePerturbation @@ -48,7 +48,7 @@ NumberNewSamples = 2 data = read_check_data_file(data_file) nCand = 500 * data.dim - phifunction = RbfType.CUBIC + phifunction = RbfKernel.CUBIC m = 2 * (data.dim + 1) numstart = ( 0 # collect all objective function values of the current trial here @@ -94,7 +94,7 @@ print(samples) print("LocalStochRBFstop Start") - rbfModel = RbfModel(phifunction, filter=MedianLpfFilter()) + rbfModel = RbfModel(kernel=phifunction, filter=MedianLpfFilter()) optres = stochastic_response_surface( data.objfunction, diff --git a/examples/single_obj_rbf/compareLearningStrategies.py b/examples/single_obj_rbf/compareLearningStrategies.py index 161e1f62..c526828f 100644 --- a/examples/single_obj_rbf/compareLearningStrategies.py +++ b/examples/single_obj_rbf/compareLearningStrategies.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False @@ -45,7 +45,7 @@ def read_and_run( maxeval: int = 0, Ntrials: int = 0, NumberNewSamples: int = 0, - rbf_type: rbf.RbfType = rbf.RbfType.CUBIC, + rbf_type: rbf.RbfKernel = rbf.RbfKernel.CUBIC, filter: rbf.RbfFilter = rbf.RbfFilter(), optim_func=optimize.multistart_stochastic_response_surface, ) -> list[optimize.OptimizeResult]: @@ -69,7 +69,7 @@ def read_and_run( Number of trials. NumberNewSamples : int, optional Number of new samples per step of the optimization algorithm. - rbf_type : rbf.RbfType, optional + rbf_type : rbf.RbfKernel, optional Type of RBF to be used. Returns @@ -88,7 +88,9 @@ def read_and_run( optres = [] for j in range(Ntrials): # Create empty RBF model - rbfModel = rbf.RbfModel(rbf_type, data.iindex, filter=filter) + rbfModel = rbf.RbfModel( + kernel=rbf_type, iindex=data.iindex, filter=filter + ) acquisitionFuncIter = deepcopy(acquisitionFunc) # # Uncomment to compare with Surrogates.jl @@ -334,7 +336,7 @@ def check_set_parameters( maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, ) if 4 in comparisonList: optresList[4] = read_and_run( @@ -354,7 +356,7 @@ def check_set_parameters( maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, optim_func=optimize.stochastic_response_surface, ) if 5 in comparisonList: @@ -366,7 +368,7 @@ def check_set_parameters( maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, optim_func=optimize.target_value_optimization, ) if 6 in comparisonList: @@ -388,7 +390,7 @@ def check_set_parameters( maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, optim_func=optimize.cptv, ) if 7 in comparisonList: @@ -410,7 +412,7 @@ def check_set_parameters( maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, optim_func=optimize.cptvl, ) if 8 in comparisonList: @@ -422,7 +424,7 @@ def check_set_parameters( maxeval=100, Ntrials=3, NumberNewSamples=10, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, optim_func=optimize.target_value_optimization, ) diff --git a/examples/single_obj_rbf/data.py b/examples/single_obj_rbf/data.py index 5183a485..4c0b0c58 100644 --- a/examples/single_obj_rbf/data.py +++ b/examples/single_obj_rbf/data.py @@ -26,7 +26,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from dataclasses import dataclass diff --git a/examples/single_obj_rbf/datainput_Branin.py b/examples/single_obj_rbf/datainput_Branin.py index 56f46ba9..2bf5cd56 100644 --- a/examples/single_obj_rbf/datainput_Branin.py +++ b/examples/single_obj_rbf/datainput_Branin.py @@ -26,7 +26,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/examples/single_obj_rbf/datainput_BraninWithInteger.py b/examples/single_obj_rbf/datainput_BraninWithInteger.py index 056f2cf1..7a2d526c 100644 --- a/examples/single_obj_rbf/datainput_BraninWithInteger.py +++ b/examples/single_obj_rbf/datainput_BraninWithInteger.py @@ -26,7 +26,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/examples/single_obj_rbf/datainput_hartman3.py b/examples/single_obj_rbf/datainput_hartman3.py index da39411a..9067c1d7 100644 --- a/examples/single_obj_rbf/datainput_hartman3.py +++ b/examples/single_obj_rbf/datainput_hartman3.py @@ -26,7 +26,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/examples/single_obj_rbf/datainput_rastrigin.py b/examples/single_obj_rbf/datainput_rastrigin.py index 36ceeae7..2264c336 100644 --- a/examples/single_obj_rbf/datainput_rastrigin.py +++ b/examples/single_obj_rbf/datainput_rastrigin.py @@ -26,7 +26,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/examples/single_obj_rbf/optprogram1.py b/examples/single_obj_rbf/optprogram1.py index a12c351d..6f406538 100644 --- a/examples/single_obj_rbf/optprogram1.py +++ b/examples/single_obj_rbf/optprogram1.py @@ -31,7 +31,7 @@ "Haoyu Jia", "Weslley S. Pereira", ] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False @@ -57,7 +57,7 @@ def read_and_run( maxeval: int = 0, Ntrials: int = 0, NumberNewSamples: int = 0, - rbf_type: rbf.RbfType = rbf.RbfType.CUBIC, + rbf_type: rbf.RbfKernel = rbf.RbfKernel.CUBIC, filter: rbf.RbfFilter = rbf.RbfFilter(), PlotResult: bool = True, optim_func=optimize.multistart_stochastic_response_surface, @@ -82,7 +82,7 @@ def read_and_run( Number of trials. NumberNewSamples : int, optional Number of new samples per step of the optimization algorithm. - rbf_type : rbf.RbfType, optional + rbf_type : rbf.RbfKernel, optional Type of RBF to be used. PlotResult : bool, optional Plot the results. @@ -103,7 +103,9 @@ def read_and_run( optres = [] for j in range(Ntrials): # Create empty RBF model - rbfModel = rbf.RbfModel(rbf_type, data.iindex, filter=filter) + rbfModel = rbf.RbfModel( + kernel=rbf_type, iindex=data.iindex, filter=filter + ) acquisitionFuncIter = deepcopy(acquisitionFunc) # # Uncomment to compare with Surrogates.jl @@ -366,7 +368,7 @@ def main(args): maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, PlotResult=True, ) elif args.config == 4: @@ -387,7 +389,7 @@ def main(args): maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, PlotResult=True, optim_func=optimize.stochastic_response_surface, ) @@ -399,7 +401,7 @@ def main(args): maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, PlotResult=True, optim_func=optimize.target_value_optimization, ) @@ -421,7 +423,7 @@ def main(args): maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, PlotResult=True, optim_func=optimize.cptv, ) @@ -443,7 +445,7 @@ def main(args): maxeval=100, Ntrials=3, NumberNewSamples=1, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, PlotResult=True, optim_func=optimize.cptvl, ) @@ -455,7 +457,7 @@ def main(args): maxeval=100, Ntrials=3, NumberNewSamples=10, - rbf_type=rbf.RbfType.THINPLATE, + rbf_type=rbf.RbfKernel.THINPLATE, PlotResult=True, optim_func=optimize.target_value_optimization, ) diff --git a/examples/vlse_benchmark/vlse_bench.ipynb b/examples/vlse_benchmark/vlse_bench.ipynb index d4e57b75..23f35fe3 100644 --- a/examples/vlse_benchmark/vlse_bench.ipynb +++ b/examples/vlse_benchmark/vlse_bench.ipynb @@ -28,7 +28,7 @@ "__maintainer__ = \"Weslley S. Pereira\"\n", "__email__ = \"weslley.dasilvapereira@nrel.gov\"\n", "__credits__ = [\"Weslley S. Pereira\"]\n", - "__version__ = \"0.4.1\"\n", + "__version__ = \"0.4.2\"\n", "__deprecated__ = False" ] }, diff --git a/examples/vlse_benchmark/vlse_bench.py b/examples/vlse_benchmark/vlse_bench.py index d9c28816..7357643b 100644 --- a/examples/vlse_benchmark/vlse_bench.py +++ b/examples/vlse_benchmark/vlse_bench.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import os @@ -29,6 +29,7 @@ import time from tests.test_vlse_bench import run_optimizer from blackboxopt import optimize, acquisition, sampling +from pathlib import Path # Functions to be tested myRfuncs = ( @@ -198,9 +199,11 @@ tf = time.time() # Save the results + folder = os.path.dirname(os.path.abspath(__file__)) + "/pickle" + Path(folder).mkdir(parents=True, exist_ok=True) filepath = ( - os.path.dirname(os.path.abspath(__file__)) - + "/pickle/vlse_bench_plot_" + folder + + "/" + args.problem + "_" + args.algorithm diff --git a/pyproject.toml b/pyproject.toml index ff407e5b..f3c37275 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "blackboxopt" -version = "0.4.1" +version = "0.4.2" description = "Surrogate models and active learning for scientific applications" authors = [ {name = "Weslley da Silva Pereira", email = "weslley.dasilvapereira@nrel.gov"}, diff --git a/tests/gosac_benchmark.py b/tests/gosac_benchmark.py index 9ddeb0bc..66b9ab03 100644 --- a/tests/gosac_benchmark.py +++ b/tests/gosac_benchmark.py @@ -419,8 +419,8 @@ class Problem: lambda x: np.reshape(fRana(x) - 5, (-1, 1)), (0,), ((-9, 9), (-3 * np.pi, 3 * np.pi)), - (-8, -9.4142), - -104.3309, + # (-8, -9.4142), + # -104.3309, ) ) @@ -487,20 +487,20 @@ class Problem: # Problem 16 gosac_p.append(copy(gosac_p[3])) gosac_p[-1].iindex = () -gosac_p[-1].xmin = ( - 0.4849, - 0.8024, - 0.4851, - 0.0690, - 0.4849, - 0.8024, - 0.4851, - 0.9690, - 0.2310, - 0.9976, - 1.9990, -) -gosac_p[-1].fmin = 3.8853 +# gosac_p[-1].xmin = ( +# 0.4849, +# 0.8024, +# 0.4851, +# 0.0690, +# 0.4849, +# 0.8024, +# 0.4851, +# 0.9690, +# 0.2310, +# 0.9976, +# 1.9990, +# ) +# gosac_p[-1].fmin = 3.8853 # Problem 17 gosac_p.append(copy(gosac_p[2])) @@ -700,8 +700,8 @@ class Problem: lambda x: np.reshape(-np.sum(np.sin(10 * np.log(x)), axis=1), (-1, 1)), (), ((0.25, np.pi),) * 10, - (3, 2, 3, 2, 3, 2, 3, 3, 1, 2), - 1.1783, + # (3, 2, 3, 2, 3, 2, 3, 3, 1, 2), + # 1.1783, ) ) diff --git a/tests/test_gosac_bench.py b/tests/test_gosac_bench.py index 3a98a0fc..69fe0b1b 100644 --- a/tests/test_gosac_bench.py +++ b/tests/test_gosac_bench.py @@ -54,21 +54,22 @@ def test_gosac(problem: gosacbmk.Problem) -> None: @pytest.mark.parametrize("problem", gosacbmk.gosac_p) def test_benchmark(problem: gosacbmk.Problem) -> None: - print(problem.xmin) - print(problem.fmin) - print(problem.objf(np.asarray([problem.xmin]))[0]) - print(problem.gfun(np.asarray([problem.xmin]))[0]) - - if problem.fmin is not None and problem.fmin != 0: - assert ( - abs(problem.objf(np.asarray([problem.xmin]))[0] - problem.fmin) - / abs(problem.fmin) - <= 1e-3 - ) - else: - assert problem.objf(np.asarray([problem.xmin]))[0] == problem.fmin - - assert np.all(problem.gfun(np.asarray([problem.xmin]))[0] <= 0) + if problem.xmin is not None and problem.fmin is not None: + print(problem.xmin) + print(problem.fmin) + print(problem.objf(np.asarray([problem.xmin]))[0]) + print(problem.gfun(np.asarray([problem.xmin]))[0]) + + if problem.fmin != 0: + assert ( + abs(problem.objf(np.asarray([problem.xmin]))[0] - problem.fmin) + / abs(problem.fmin) + <= 1e-2 + ) + else: + assert problem.objf(np.asarray([problem.xmin]))[0] == problem.fmin + + assert np.all(problem.gfun(np.asarray([problem.xmin]))[0] <= 1e-2) if __name__ == "__main__": diff --git a/tests/test_optimize.py b/tests/test_optimize.py index 8be4242d..e482c7cc 100644 --- a/tests/test_optimize.py +++ b/tests/test_optimize.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/tests/test_rbf.py b/tests/test_rbf.py index 389b328a..11312fd6 100644 --- a/tests/test_rbf.py +++ b/tests/test_rbf.py @@ -20,44 +20,51 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np -import sys import pytest -from blackboxopt.rbf import MedianLpfFilter, RbfModel, RbfType +from blackboxopt.rbf import MedianLpfFilter, RbfModel +from blackboxopt.rbf_kernel import RbfKernel, KERNEL_FUNC + + +@pytest.mark.parametrize("kernel", [k for k in RbfKernel]) +def test_kernel(kernel): + testInput = [1.0, 2.0, 3.0, 4.0] + testResults = { + RbfKernel.LINEAR: [-1.0, -2.0, -3.0, -4.0], + RbfKernel.CUBIC: [1.0, 8.0, 27.0, 64.0], + RbfKernel.THINPLATE: [0.0, 2.7725887, 9.88751, 22.18071], + RbfKernel.QUINTIC: [-1.0, -32.0, -243.0, -1024.0], + RbfKernel.MULTIQUADRIC: [ + -1.4142135, + -2.236068, + -3.1622777, + -4.1231055, + ], + RbfKernel.INVERSE_MULTIQUADRIC: [ + 0.70710677, + 0.4472136, + 0.31622776, + 0.24253562, + ], + RbfKernel.INVERSE_QUADRATIC: [0.5, 0.2, 0.1, 0.05882353], + RbfKernel.GAUSSIAN: [ + 3.67879450e-01, + 1.83156393e-02, + 1.23409802e-04, + 1.12535176e-07, + ], + } + phi = KERNEL_FUNC[kernel] + + np.testing.assert_array_almost_equal(phi(testInput), testResults[kernel]) class TestRbfModel: rbf_model = RbfModel() - def test_phi(self): - self.rbf_model.type = RbfType.LINEAR - r_linear = np.array([1.0, 2.0, 3.0]) - result_linear = self.rbf_model.phi(r_linear) - expected_linear = np.array([1.0, 2.0, 3.0]) - np.testing.assert_array_equal(np.array(result_linear), expected_linear) - assert self.rbf_model.phi(4.0) == 4.0 - - self.rbf_model.type = RbfType.CUBIC - r_cubic = np.array([1.0, 2.0, 3.0]) - result_cubic = self.rbf_model.phi(r_cubic) - expected_cubic = np.array([1.0, 8.0, 27.0]) - np.testing.assert_array_equal(np.array(result_cubic), expected_cubic) - assert self.rbf_model.phi(4.0) == 64.0 - - self.rbf_model.type = RbfType.THINPLATE - r_thinplate = np.array([1.0, 2.0, 3.0]) - result_thinplate = self.rbf_model.phi(r_thinplate) - expected_thinplate = np.array([0.0, 2.77258872, 9.8875106]) - np.testing.assert_allclose( - np.array(result_thinplate), expected_thinplate - ) - assert self.rbf_model.phi(4.0) == ( - 4 * 4 * np.log(4 + sys.float_info.min) - ) - def test_dim(self): assert self.rbf_model.dim() == 0 diff --git a/tests/test_sampling.py b/tests/test_sampling.py index 221ae184..5068c5d5 100644 --- a/tests/test_sampling.py +++ b/tests/test_sampling.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False import numpy as np diff --git a/tests/test_vlse_bench.py b/tests/test_vlse_bench.py index 747d1b0c..28064162 100644 --- a/tests/test_vlse_bench.py +++ b/tests/test_vlse_bench.py @@ -20,7 +20,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from copy import deepcopy @@ -139,7 +139,7 @@ def objf(x: np.ndarray) -> np.ndarray: # Surrogate model with median low-pass filter rbfModel = rbf.RbfModel( - rbf.RbfType.CUBIC, iindex, filter=rbf.MedianLpfFilter() + kernel=rbf.RbfKernel.CUBIC, iindex=iindex, filter=rbf.MedianLpfFilter() ) # Update acquisition strategy, using maxEval and nArgs for the problem diff --git a/tests/vlse_benchmark/__init__.py b/tests/vlse_benchmark/__init__.py index 13ad92f2..129511d6 100644 --- a/tests/vlse_benchmark/__init__.py +++ b/tests/vlse_benchmark/__init__.py @@ -27,7 +27,7 @@ __maintainer__ = "Weslley S. Pereira" __email__ = "weslley.dasilvapereira@nrel.gov" __credits__ = ["Sonja Surjanovic", "Derek Bingham", "Weslley S. Pereira"] -__version__ = "0.4.1" +__version__ = "0.4.2" __deprecated__ = False from rpy2.robjects import r