28954437142

Committed 08 Jul 2026 03:24PM UTC coverage: 87.552% (+0.3%) from 87.28%

Build # 28954437142

Build Type

Pull #23

github

Committed by

web-flow

Commit Message

Merge 15e9c59e8 into acbef679d

Pull Request Pull Request #23: feat: add MO-SEGO acquisition strategy

Coverage Stats

665 of 773 new or added lines in 17 files covered. (86.03%)

3 existing lines in 1 file now uncovered.

4403 of 5029 relevant lines covered (87.55%)

1.75 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

64.12

/src/smt_optim/utils/multi_obj.py

import numpy as np
from scipy.spatial.distance import cdist

from moocore import hypervolume as moocore_hv
from pymoo.core.problem import Problem as PymooProblem


def get_pareto_mask(Y: np.ndarray) -> np.ndarray:

    n = Y.shape[0]
    is_pareto = np.ones(n, dtype=bool)

    for i in range(n):
        if not is_pareto[i]:
            continue

        # A point is dominated if another point is <= in all objectives
        # and strictly < in at least one
        dominates = np.all(Y <= Y[i], axis=1) & np.any(Y < Y[i], axis=1)

        # If any point dominates i -> i is not Pareto
        if np.any(dominates):
            is_pareto[i] = False

    return is_pareto


def get_pareto_front(Y: np.ndarray) -> np.ndarray:
    """
    Return the non-dominated objective vectors from ``Y``.

    Parameters
    ----------
    Y : np.ndarray
        Array of shape ``(n_samples, n_objectives)`` containing objective
        values for each sample.

    Returns
    -------
    np.ndarray
        Array of shape ``(n_pareto, n_objectives)`` containing the
        non-dominated objective vectors (the Pareto front).

    Notes
    -----
    Assumes a minimization problem for all objectives and no constraints.
    """

    pareto_mask = get_pareto_mask(Y)
    pareto = Y[pareto_mask]
    return pareto


def get_pf_from_dataset(
    dataset, ctol: float = 1e-4, fid: int = -1, return_dict: bool = False
) -> np.ndarray | dict:

    data_dict = dataset.export_as_dict()
    obj = data_dict["obj"]
    rscv = data_dict["rscv"]
    fidelity = data_dict["fidelity"]

    feas_mask = rscv <= ctol

    if fid == -1:
        fid = np.max(fidelity)

    fid_mask = fidelity == fid

    fid_feas_mask = fid_mask & feas_mask

    pareto_mask = get_pareto_mask(obj[fid_feas_mask])

    pareto_front = obj[fid_feas_mask][pareto_mask]

    if return_dict:
        data = {
            "x": data_dict["x"][fid_feas_mask][pareto_mask],
            "obj": pareto_front,
            "cstr": data_dict["cstr"][fid_feas_mask][pareto_mask],
            "rscv": rscv[fid_feas_mask][pareto_mask],
            "fidelity": fidelity[fid_feas_mask][pareto_mask],
        }

        return data

    return pareto_front


def hypervolume_2d(pf: np.ndarray, ref: np.ndarray) -> float:
    """
    Compute the 2D hypervolume indicator  he hypervolume of the Pareto front.

    Parameters
    ----------
    pf: np.ndarray of shape (num_points, 2)
        Pareto front.

    ref: np.ndarray of shape (2, )
        Reference objective values.

    Returns
    -------
    float
        Hypervolume indicator value.

    Notes:
        Assume both objective are minimized.
    """

    if pf.shape[1] != 2 or ref.shape[0] != 2:
        raise Exception(
            "Current hypervolume implementation is only for bi-objective optimization."
        )

    sorted_idx = np.argsort(pf[:, 0])
    sorted_pf = pf[sorted_idx]

    hv = 0.0
    prev_f2 = ref[1]

    for idx in range(sorted_pf.shape[0]):
        f1 = sorted_pf[idx, 0]
        f2 = sorted_pf[idx, 1]

        width = ref[0] - f1
        height = prev_f2 - f2

        if width > 0 and height > 0:
            hv += width * height

        prev_f2 = min(prev_f2, f2)

    return hv


def hypervolume(pf: np.ndarray, ref: np.ndarray) -> float:
    """
    Compute the hypervolume indicator of the Pareto front.

    Uses the `moocore` implementation:
    https://multi-objective.github.io/moocore/python/reference/generated/moocore.hypervolume.html

    Parameters
    ----------
    pf: np.ndarray of shape (num_points, 2)
        Pareto front.

    ref: np.ndarray of shape (2, )
        Reference objective values.

    Returns
    -------
    float
        Hypervolume indicator value.

    Notes:
        Assume both objective are minimized.
    """

    return moocore_hv(pf, ref=ref)


def spacing(pf: np.ndarray) -> float:
    """
    Compute the spacing indicator of the Pareto front (Schott, 1995). A lower value is better.

    Parameters
    ----------
    pf: np.ndarray of shape (num_points, num_objectives)
        Pareto front.

    Returns
    -------
    float
        Spacing indicator value.

    Notes:
        Assume both objective are minimized.
    """

    num_pf = pf.shape[0]

    if num_pf <= 1:
        return np.nan

    distances = cdist(pf, pf, "cityblock")
    np.fill_diagonal(distances, np.inf)

    d1 = np.min(distances, axis=1)
    d1_mean = np.mean(d1)

    value = np.sqrt(1 / (num_pf - 1) * np.sum((d1_mean - d1) ** 2))

    return value


class PymooStateWrapper(PymooProblem):
    def __init__(self, state, scaled: bool = False, train: bool = True):
        """
        Create a Pymoo problem instance for use with Pymoo algorithms, such as NSGA-II for multi-objective optimization.

        The quantities of interest are modeled using the surrogate models.

        Parameters
        ----------
        state : State
        scaled : bool
            Normalize the input variables and standardize the quantities of interest
            (before training the surrogate models, if applicable).
        train : bool
            Train the surrogate models.

        Notes
        -----
        When using `PymooStateWrapper` after a Bayesian optimization (BO) iteration,
        `train` should be set to `True`, as the surrogate models are not trained on
        the most recently evaluated infill point.
        """

        self.state = state
        self.scaled = scaled
        self.train = train

        if not self.state.problem.design_space.is_all_cont:
            raise ValueError(
                "PymooStateWrapper currently requires a continuous optimization problem."
            )

        if self.train:
            self.state.scale_dataset(self.scaled)
            self.state.build_models()

        prob = self.state.problem

        if not self.scaled:
            l_bounds = []
            u_bounds = []

            for idx, var in enumerate(prob.design_space.design_variables):
                l_bounds.append(var.lower)
                u_bounds.append(var.upper)

            self.l_bounds = np.array(l_bounds)
            self.u_bounds = np.array(u_bounds)

        else:
            self.l_bounds = np.zeros(prob.num_dim)
            self.u_bounds = np.ones(prob.num_dim)

        self.f_callables = []
        self.g_callables = []
        self.num_g = 0
        self.h_callables = []
        self.num_h = 0

        ineq_bounds = []
        self.num_g_lower = 0
        self.num_g_upper = 0

        for o_id, o_config in enumerate(prob.obj_configs):
            self.f_callables.append(
                lambda x, f=self.state.obj_models[o_id].predict_values: f(x).ravel()
            )

        for c_id, c_config in enumerate(prob.cstr_configs):
            if c_config.equal is not None:
                self.h_callables.append(
                    lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.equal: (
                        f(x).ravel() - self.state.cstr_equal[c_id]
                    )
                )
                self.num_h += 1

            else:
                self.g_callables.append(
                    lambda x, f=self.state.cstr_models[c_id].predict_values: f(
                        x
                    ).ravel()
                )

                ineq_bounds.append(np.full(2, np.nan))

                if c_config.lower is not None:
                    ineq_bounds[-1][0] = self.state.cstr_lower[c_id]
                    self.num_g += 1

                if c_config.upper is not None:
                    ineq_bounds[-1][1] = self.state.cstr_upper[c_id]
                    self.num_g += 1

        if len(self.g_callables):
            self.ineq_bounds = np.array(ineq_bounds)

            self.g_lower_mask = np.where(np.isnan(self.ineq_bounds[:, 0]), False, True)
            self.g_upper_mask = np.where(np.isnan(self.ineq_bounds[:, 1]), False, True)

        super().__init__(
            n_var=prob.num_dim,
            n_obj=prob.num_obj,
            n_eq_constr=self.num_h,
            n_ieq_constr=self.num_g,
            xl=self.l_bounds,
            xu=self.u_bounds,
        )

    def _evaluate(self, x, out, *args, **kwargs):

        num_pt = x.shape[0]

        x_scaled = (x - self.l_bounds) / (self.u_bounds - self.l_bounds)

        # sample objectives
        out["F"] = np.full((num_pt, self.n_obj), np.nan)
        for o_idx in range(self.n_obj):
            out["F"][:, o_idx] = self.f_callables[o_idx](x_scaled).ravel()

        # sample equality constraints
        if self.n_eq_constr > 0:
            out["H"] = np.empty((num_pt, self.n_eq_constr))
            for h_idx in range(self.n_eq_constr):
                out["H"][:, h_idx] = self.h_callables[h_idx](x_scaled).ravel()

        if self.n_ieq_constr > 0:
            g_lower = np.empty((num_pt, self.ineq_bounds.shape[0]))
            g_upper = np.empty((num_pt, self.ineq_bounds.shape[0]))

            # sample inequality constraints
            if self.n_ieq_constr > 0:
                for g_idx, g_call in enumerate(self.g_callables):
                    g_vals = self.g_callables[g_idx](x_scaled).ravel()

                    if not np.isnan(self.ineq_bounds[g_idx, 0]):
                        g_lower[:, g_idx] = self.ineq_bounds[g_idx, 0] - g_vals

                    elif not np.isnan(self.ineq_bounds[g_idx, 1]):
                        g_upper[:, g_idx] = g_vals - self.ineq_bounds[g_idx, 1]

                out["G"] = np.hstack(
                    (g_lower[:, self.g_lower_mask], g_upper[:, self.g_upper_mask])
                )

1	import numpy as np	2✔
2	from scipy.spatial.distance import cdist	2✔
3
4	from moocore import hypervolume as moocore_hv	2✔
5	from pymoo.core.problem import Problem as PymooProblem	2✔
6
7
8	def get_pareto_mask(Y: np.ndarray) -> np.ndarray:	2✔
9
10	n = Y.shape[0]	2✔
11	is_pareto = np.ones(n, dtype=bool)	2✔
12
13	for i in range(n):	2✔
14	if not is_pareto[i]:	2✔
NEW 15	continue	×
16
17	# A point is dominated if another point is <= in all objectives
18	# and strictly < in at least one
19	dominates = np.all(Y <= Y[i], axis=1) & np.any(Y < Y[i], axis=1)	2✔
20
21	# If any point dominates i -> i is not Pareto
22	if np.any(dominates):	2✔
23	is_pareto[i] = False	2✔
24
25	return is_pareto	2✔
26
27
28	def get_pareto_front(Y: np.ndarray) -> np.ndarray:	2✔
29	"""
30	Return the non-dominated objective vectors from ``Y``.
31
32	Parameters
33	----------
34	Y : np.ndarray
35	Array of shape ``(n_samples, n_objectives)`` containing objective
36	values for each sample.
37
38	Returns
39	-------
40	np.ndarray
41	Array of shape ``(n_pareto, n_objectives)`` containing the
42	non-dominated objective vectors (the Pareto front).
43
44	Notes
45	-----
46	Assumes a minimization problem for all objectives and no constraints.
47	"""
48
NEW 49	pareto_mask = get_pareto_mask(Y)	×
NEW 50	pareto = Y[pareto_mask]	×
NEW 51	return pareto	×
52
53
54	def get_pf_from_dataset(	2✔
55	dataset, ctol: float = 1e-4, fid: int = -1, return_dict: bool = False
56	) -> np.ndarray \| dict:
57
58	data_dict = dataset.export_as_dict()	2✔
59	obj = data_dict["obj"]	2✔
60	rscv = data_dict["rscv"]	2✔
61	fidelity = data_dict["fidelity"]	2✔
62
63	feas_mask = rscv <= ctol	2✔
64
65	if fid == -1:	2✔
66	fid = np.max(fidelity)	2✔
67
68	fid_mask = fidelity == fid	2✔
69
70	fid_feas_mask = fid_mask & feas_mask	2✔
71
72	pareto_mask = get_pareto_mask(obj[fid_feas_mask])	2✔
73
74	pareto_front = obj[fid_feas_mask][pareto_mask]	2✔
75
76	if return_dict:	2✔
77	data = {	2✔
78	"x": data_dict["x"][fid_feas_mask][pareto_mask],
79	"obj": pareto_front,
80	"cstr": data_dict["cstr"][fid_feas_mask][pareto_mask],
81	"rscv": rscv[fid_feas_mask][pareto_mask],
82	"fidelity": fidelity[fid_feas_mask][pareto_mask],
83	}
84
85	return data	2✔
86
87	return pareto_front	2✔
88
89
90	def hypervolume_2d(pf: np.ndarray, ref: np.ndarray) -> float:	2✔
91	"""
92	Compute the 2D hypervolume indicator he hypervolume of the Pareto front.
93
94	Parameters
95	----------
96	pf: np.ndarray of shape (num_points, 2)
97	Pareto front.
98
99	ref: np.ndarray of shape (2, )
100	Reference objective values.
101
102	Returns
103	-------
104	float
105	Hypervolume indicator value.
106
107	Notes:
108	Assume both objective are minimized.
109	"""
110
NEW 111	if pf.shape[1] != 2 or ref.shape[0] != 2:	×
NEW 112	raise Exception(	×
113	"Current hypervolume implementation is only for bi-objective optimization."
114	)
115
NEW 116	sorted_idx = np.argsort(pf[:, 0])	×
NEW 117	sorted_pf = pf[sorted_idx]	×
118
NEW 119	hv = 0.0	×
NEW 120	prev_f2 = ref[1]	×
121
NEW 122	for idx in range(sorted_pf.shape[0]):	×
NEW 123	f1 = sorted_pf[idx, 0]	×
NEW 124	f2 = sorted_pf[idx, 1]	×
125
NEW 126	width = ref[0] - f1	×
NEW 127	height = prev_f2 - f2	×
128
NEW 129	if width > 0 and height > 0:	×
NEW 130	hv += width * height	×
131
NEW 132	prev_f2 = min(prev_f2, f2)	×
133
NEW 134	return hv	×
135
136
137	def hypervolume(pf: np.ndarray, ref: np.ndarray) -> float:	2✔
138	"""
139	Compute the hypervolume indicator of the Pareto front.
140
141	Uses the `moocore` implementation:
142	https://multi-objective.github.io/moocore/python/reference/generated/moocore.hypervolume.html
143
144	Parameters
145	----------
146	pf: np.ndarray of shape (num_points, 2)
147	Pareto front.
148
149	ref: np.ndarray of shape (2, )
150	Reference objective values.
151
152	Returns
153	-------
154	float
155	Hypervolume indicator value.
156
157	Notes:
158	Assume both objective are minimized.
159	"""
160
NEW 161	return moocore_hv(pf, ref=ref)	×
162
163
164	def spacing(pf: np.ndarray) -> float:	2✔
165	"""
166	Compute the spacing indicator of the Pareto front (Schott, 1995). A lower value is better.
167
168	Parameters
169	----------
170	pf: np.ndarray of shape (num_points, num_objectives)
171	Pareto front.
172
173	Returns
174	-------
175	float
176	Spacing indicator value.
177
178	Notes:
179	Assume both objective are minimized.
180	"""
181
NEW 182	num_pf = pf.shape[0]	×
183
NEW 184	if num_pf <= 1:	×
NEW 185	return np.nan	×
186
NEW 187	distances = cdist(pf, pf, "cityblock")	×
NEW 188	np.fill_diagonal(distances, np.inf)	×
189
NEW 190	d1 = np.min(distances, axis=1)	×
NEW 191	d1_mean = np.mean(d1)	×
192
NEW 193	value = np.sqrt(1 / (num_pf - 1) * np.sum((d1_mean - d1) ** 2))	×
194
NEW 195	return value	×
196
197
198	class PymooStateWrapper(PymooProblem):	2✔
199	def __init__(self, state, scaled: bool = False, train: bool = True):	2✔
200	"""
201	Create a Pymoo problem instance for use with Pymoo algorithms, such as NSGA-II for multi-objective optimization.
202
203	The quantities of interest are modeled using the surrogate models.
204
205	Parameters
206	----------
207	state : State
208	scaled : bool
209	Normalize the input variables and standardize the quantities of interest
210	(before training the surrogate models, if applicable).
211	train : bool
212	Train the surrogate models.
213
214	Notes
215	-----
216	When using `PymooStateWrapper` after a Bayesian optimization (BO) iteration,
217	`train` should be set to `True`, as the surrogate models are not trained on
218	the most recently evaluated infill point.
219	"""
220
221	self.state = state	2✔
222	self.scaled = scaled	2✔
223	self.train = train	2✔
224
225	if not self.state.problem.design_space.is_all_cont:	2✔
NEW 226	raise ValueError(	×
227	"PymooStateWrapper currently requires a continuous optimization problem."
228	)
229
230	if self.train:	2✔
NEW 231	self.state.scale_dataset(self.scaled)	×
NEW 232	self.state.build_models()	×
233
234	prob = self.state.problem	2✔
235
236	if not self.scaled:	2✔
NEW 237	l_bounds = []	×
NEW 238	u_bounds = []	×
239
NEW 240	for idx, var in enumerate(prob.design_space.design_variables):	×
NEW 241	l_bounds.append(var.lower)	×
NEW 242	u_bounds.append(var.upper)	×
243
NEW 244	self.l_bounds = np.array(l_bounds)	×
NEW 245	self.u_bounds = np.array(u_bounds)	×
246
247	else:
248	self.l_bounds = np.zeros(prob.num_dim)	2✔
249	self.u_bounds = np.ones(prob.num_dim)	2✔
250
251	self.f_callables = []	2✔
252	self.g_callables = []	2✔
253	self.num_g = 0	2✔
254	self.h_callables = []	2✔
255	self.num_h = 0	2✔
256
257	ineq_bounds = []	2✔
258	self.num_g_lower = 0	2✔
259	self.num_g_upper = 0	2✔
260
261	for o_id, o_config in enumerate(prob.obj_configs):	2✔
262	self.f_callables.append(	2✔
263	lambda x, f=self.state.obj_models[o_id].predict_values: f(x).ravel()
264	)
265
266	for c_id, c_config in enumerate(prob.cstr_configs):	2✔
267	if c_config.equal is not None:	2✔
NEW 268	self.h_callables.append(	×
269	lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.equal: (
270	f(x).ravel() - self.state.cstr_equal[c_id]
271	)
272	)
NEW 273	self.num_h += 1	×
274
275	else:
276	self.g_callables.append(	2✔
277	lambda x, f=self.state.cstr_models[c_id].predict_values: f(
278	x
279	).ravel()
280	)
281
282	ineq_bounds.append(np.full(2, np.nan))	2✔
283
284	if c_config.lower is not None:	2✔
NEW 285	ineq_bounds[-1][0] = self.state.cstr_lower[c_id]	×
NEW 286	self.num_g += 1	×
287
288	if c_config.upper is not None:	2✔
289	ineq_bounds[-1][1] = self.state.cstr_upper[c_id]	2✔
290	self.num_g += 1	2✔
291
292	if len(self.g_callables):	2✔
293	self.ineq_bounds = np.array(ineq_bounds)	2✔
294
295	self.g_lower_mask = np.where(np.isnan(self.ineq_bounds[:, 0]), False, True)	2✔
296	self.g_upper_mask = np.where(np.isnan(self.ineq_bounds[:, 1]), False, True)	2✔
297
298	super().__init__(	2✔
299	n_var=prob.num_dim,
300	n_obj=prob.num_obj,
301	n_eq_constr=self.num_h,
302	n_ieq_constr=self.num_g,
303	xl=self.l_bounds,
304	xu=self.u_bounds,
305	)
306
307	def _evaluate(self, x, out, args, *kwargs):	2✔
308
309	num_pt = x.shape[0]	2✔
310
311	x_scaled = (x - self.l_bounds) / (self.u_bounds - self.l_bounds)	2✔
312
313	# sample objectives
314	out["F"] = np.full((num_pt, self.n_obj), np.nan)	2✔
315	for o_idx in range(self.n_obj):	2✔
316	out["F"][:, o_idx] = self.f_callables[o_idx](x_scaled).ravel()	2✔
317
318	# sample equality constraints
319	if self.n_eq_constr > 0:	2✔
NEW 320	out["H"] = np.empty((num_pt, self.n_eq_constr))	×
NEW 321	for h_idx in range(self.n_eq_constr):	×
NEW 322	out["H"][:, h_idx] = self.h_callables[h_idx](x_scaled).ravel()	×
323
324	if self.n_ieq_constr > 0:	2✔
325	g_lower = np.empty((num_pt, self.ineq_bounds.shape[0]))	2✔
326	g_upper = np.empty((num_pt, self.ineq_bounds.shape[0]))	2✔
327
328	# sample inequality constraints
329	if self.n_ieq_constr > 0:	2✔
330	for g_idx, g_call in enumerate(self.g_callables):	2✔
331	g_vals = self.g_callables[g_idx](x_scaled).ravel()	2✔
332
333	if not np.isnan(self.ineq_bounds[g_idx, 0]):	2✔
NEW 334	g_lower[:, g_idx] = self.ineq_bounds[g_idx, 0] - g_vals	×
335
336	elif not np.isnan(self.ineq_bounds[g_idx, 1]):	2✔
337	g_upper[:, g_idx] = g_vals - self.ineq_bounds[g_idx, 1]	2✔
338
339	out["G"] = np.hstack(	2✔
340	(g_lower[:, self.g_lower_mask], g_upper[:, self.g_upper_mask])
341	)

SMTorg / smt-optim / 28954437142

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous