28718428290

Committed 04 Jul 2026 08:17PM UTC coverage: 86.103% (-0.9%) from 87.05%

Build # 28718428290

Build Type

push

github

Committed by

web-flow

Commit Message

Feat/mo acquisition functions (#15)

Coverage Stats

270 of 363 new or added lines in 5 files covered. (74.38%)

3922 of 4555 relevant lines covered (86.1%)

1.72 hits per line

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

31.58

/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) -> np.ndarray:

    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[feas_mask] == fid

    pareto_front = get_pareto_front(obj[feas_mask][fid_mask])

    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):

        self.state = state

        self.state.scale_dataset(False)
        self.state.build_models()

        prob = self.state.problem

        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)

        self.f_callables = []
        self.g_callables = []
        self.h_callables = []

        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() - val
                    )
                )

            else:
                if c_config.lower is not None:
                    self.g_callables.append(
                        lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.lower: (
                            val - f(x).ravel()
                        )
                    )

                if c_config.upper is not None:
                    self.g_callables.append(
                        lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.upper: (
                            f(x).ravel() - val
                        )
                    )

        super().__init__(
            n_var=prob.num_dim,
            n_obj=prob.num_obj,
            n_eq_constr=len(self.h_callables),
            n_ieq_constr=len(self.g_callables),
            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)

        out["F"] = np.full((num_pt, self.n_obj), np.nan)

        if self.n_eq_constr > 0:
            out["H"] = np.empty((num_pt, self.n_eq_constr))

        if self.n_ieq_constr > 0:
            out["G"] = np.empty((num_pt, self.n_ieq_constr))

        for o_idx in range(self.n_obj):
            out["F"][:, o_idx] = self.f_callables[o_idx](x_scaled).ravel()

        for h_idx in range(self.n_eq_constr):
            out["H"][:, h_idx] = self.h_callables[h_idx](x_scaled).ravel()

        for g_idx in range(self.n_ieq_constr):
            out["G"][:, g_idx] = self.g_callables[g_idx](x_scaled).ravel()


def purity(pf: np.ndarray, ref_pf: np.ndarray) -> float:
    """
    Compute the purity indicator of the Pareto front.
    Purity is the ratio of points in the obtained Pareto front that are not
    strictly dominated by any point in the reference Pareto front.

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

    Returns
    -------
    float
        Purity indicator value (between 0.0 and 1.0).
    """
    if len(pf) == 0:
        return 0.0
    dominated = 0
    for p in pf:
        # A point is dominated if any point in ref_pf is <= in all dims and < in at least one
        is_dom = np.any(np.all(ref_pf <= p, axis=1) & np.any(ref_pf < p, axis=1))
        if is_dom:
            dominated += 1
    return 1.0 - dominated / len(pf)


def gamma_spread(pf: np.ndarray) -> float:
    """
    Compute the Gamma-spread (Maximum Spread) indicator of the Pareto front.
    It measures the maximum Euclidean distance between any two points in the
    Pareto front.

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

    Returns
    -------
    float
        Maximum spread indicator value.
    """
    if len(pf) <= 1:
        return np.nan

    from scipy.spatial.distance import pdist

    distances = pdist(pf)
    if len(distances) == 0:
        return np.nan
    return float(np.max(distances))

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
49	pareto_mask = get_pareto_mask(Y)	2✔
50	pareto = Y[pareto_mask]	2✔
51	return pareto	2✔
52
53
54	def get_pf_from_dataset(dataset, ctol: float = 1e-4, fid: int = -1) -> np.ndarray:	2✔
55
56	data_dict = dataset.export_as_dict()	2✔
57	obj = data_dict["obj"]	2✔
58	rscv = data_dict["rscv"]	2✔
59	fidelity = data_dict["fidelity"]	2✔
60
61	feas_mask = rscv <= ctol	2✔
62
63	if fid == -1:	2✔
64	fid = np.max(fidelity)	2✔
65
66	fid_mask = fidelity[feas_mask] == fid	2✔
67
68	pareto_front = get_pareto_front(obj[feas_mask][fid_mask])	2✔
69
70	return pareto_front	2✔
71
72
73	def hypervolume_2d(pf: np.ndarray, ref: np.ndarray) -> float:	2✔
74	"""
75	Compute the 2D hypervolume indicator he hypervolume of the Pareto front.
76
77	Parameters
78	----------
79	pf: np.ndarray of shape (num_points, 2)
80	Pareto front.
81
82	ref: np.ndarray of shape (2, )
83	Reference objective values.
84
85	Returns
86	-------
87	float
88	Hypervolume indicator value.
89
90	Notes:
91	Assume both objective are minimized.
92	"""
93
NEW 94	if pf.shape[1] != 2 or ref.shape[0] != 2:	×
NEW 95	raise Exception(	×
96	"Current hypervolume implementation is only for bi-objective optimization."
97	)
98
NEW 99	sorted_idx = np.argsort(pf[:, 0])	×
NEW 100	sorted_pf = pf[sorted_idx]	×
101
NEW 102	hv = 0.0	×
NEW 103	prev_f2 = ref[1]	×
104
NEW 105	for idx in range(sorted_pf.shape[0]):	×
NEW 106	f1 = sorted_pf[idx, 0]	×
NEW 107	f2 = sorted_pf[idx, 1]	×
108
NEW 109	width = ref[0] - f1	×
NEW 110	height = prev_f2 - f2	×
111
NEW 112	if width > 0 and height > 0:	×
NEW 113	hv += width * height	×
114
NEW 115	prev_f2 = min(prev_f2, f2)	×
116
NEW 117	return hv	×
118
119
120	def hypervolume(pf: np.ndarray, ref: np.ndarray) -> float:	2✔
121	"""
122	Compute the hypervolume indicator of the Pareto front.
123
124	Uses the `moocore` implementation:
125	https://multi-objective.github.io/moocore/python/reference/generated/moocore.hypervolume.html
126
127	Parameters
128	----------
129	pf: np.ndarray of shape (num_points, 2)
130	Pareto front.
131
132	ref: np.ndarray of shape (2, )
133	Reference objective values.
134
135	Returns
136	-------
137	float
138	Hypervolume indicator value.
139
140	Notes:
141	Assume both objective are minimized.
142	"""
143
NEW 144	return moocore_hv(pf, ref=ref)	×
145
146
147	def spacing(pf: np.ndarray) -> float:	2✔
148	"""
149	Compute the spacing indicator of the Pareto front (Schott, 1995). A lower value is better.
150
151	Parameters
152	----------
153	pf: np.ndarray of shape (num_points, num_objectives)
154	Pareto front.
155
156	Returns
157	-------
158	float
159	Spacing indicator value.
160
161	Notes:
162	Assume both objective are minimized.
163	"""
164
NEW 165	num_pf = pf.shape[0]	×
166
NEW 167	if num_pf <= 1:	×
NEW 168	return np.nan	×
169
NEW 170	distances = cdist(pf, pf, "cityblock")	×
NEW 171	np.fill_diagonal(distances, np.inf)	×
172
NEW 173	d1 = np.min(distances, axis=1)	×
NEW 174	d1_mean = np.mean(d1)	×
175
NEW 176	value = np.sqrt(1 / (num_pf - 1) * np.sum((d1_mean - d1) ** 2))	×
177
NEW 178	return value	×
179
180
181	class PymooStateWrapper(PymooProblem):	2✔
182	def __init__(self, state):	2✔
183
NEW 184	self.state = state	×
185
NEW 186	self.state.scale_dataset(False)	×
NEW 187	self.state.build_models()	×
188
NEW 189	prob = self.state.problem	×
190
NEW 191	l_bounds = []	×
NEW 192	u_bounds = []	×
NEW 193	for idx, var in enumerate(prob.design_space.design_variables):	×
NEW 194	l_bounds.append(var.lower)	×
NEW 195	u_bounds.append(var.upper)	×
196
NEW 197	self.l_bounds = np.array(l_bounds)	×
NEW 198	self.u_bounds = np.array(u_bounds)	×
199
NEW 200	self.f_callables = []	×
NEW 201	self.g_callables = []	×
NEW 202	self.h_callables = []	×
203
NEW 204	for o_id, o_config in enumerate(prob.obj_configs):	×
NEW 205	self.f_callables.append(	×
206	lambda x, f=self.state.obj_models[o_id].predict_values: f(x).ravel()
207	)
208
NEW 209	for c_id, c_config in enumerate(prob.cstr_configs):	×
NEW 210	if c_config.equal is not None:	×
NEW 211	self.h_callables.append(	×
212	lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.equal: (
213	f(x).ravel() - val
214	)
215	)
216
217	else:
NEW 218	if c_config.lower is not None:	×
NEW 219	self.g_callables.append(	×
220	lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.lower: (
221	val - f(x).ravel()
222	)
223	)
224
NEW 225	if c_config.upper is not None:	×
NEW 226	self.g_callables.append(	×
227	lambda x, f=self.state.cstr_models[c_id].predict_values, val=c_config.upper: (
228	f(x).ravel() - val
229	)
230	)
231
NEW 232	super().__init__(	×
233	n_var=prob.num_dim,
234	n_obj=prob.num_obj,
235	n_eq_constr=len(self.h_callables),
236	n_ieq_constr=len(self.g_callables),
237	xl=self.l_bounds,
238	xu=self.u_bounds,
239	)
240
241	def _evaluate(self, x, out, args, *kwargs):	2✔
242
NEW 243	num_pt = x.shape[0]	×
244
NEW 245	x_scaled = (x - self.l_bounds) / (self.u_bounds - self.l_bounds)	×
246
NEW 247	out["F"] = np.full((num_pt, self.n_obj), np.nan)	×
248
NEW 249	if self.n_eq_constr > 0:	×
NEW 250	out["H"] = np.empty((num_pt, self.n_eq_constr))	×
251
NEW 252	if self.n_ieq_constr > 0:	×
NEW 253	out["G"] = np.empty((num_pt, self.n_ieq_constr))	×
254
NEW 255	for o_idx in range(self.n_obj):	×
NEW 256	out["F"][:, o_idx] = self.f_callables[o_idx](x_scaled).ravel()	×
257
NEW 258	for h_idx in range(self.n_eq_constr):	×
NEW 259	out["H"][:, h_idx] = self.h_callables[h_idx](x_scaled).ravel()	×
260
NEW 261	for g_idx in range(self.n_ieq_constr):	×
NEW 262	out["G"][:, g_idx] = self.g_callables[g_idx](x_scaled).ravel()	×
263
264
265	def purity(pf: np.ndarray, ref_pf: np.ndarray) -> float:	2✔
266	"""
267	Compute the purity indicator of the Pareto front.
268	Purity is the ratio of points in the obtained Pareto front that are not
269	strictly dominated by any point in the reference Pareto front.
270
271	Parameters
272	----------
273	pf: np.ndarray of shape (num_points, num_objectives)
274	Pareto front approximation.
275	ref_pf: np.ndarray of shape (num_ref_points, num_objectives)
276	Reference Pareto front.
277
278	Returns
279	-------
280	float
281	Purity indicator value (between 0.0 and 1.0).
282	"""
NEW 283	if len(pf) == 0:	×
NEW 284	return 0.0	×
NEW 285	dominated = 0	×
NEW 286	for p in pf:	×
287	# A point is dominated if any point in ref_pf is <= in all dims and < in at least one
NEW 288	is_dom = np.any(np.all(ref_pf <= p, axis=1) & np.any(ref_pf < p, axis=1))	×
NEW 289	if is_dom:	×
NEW 290	dominated += 1	×
NEW 291	return 1.0 - dominated / len(pf)	×
292
293
294	def gamma_spread(pf: np.ndarray) -> float:	2✔
295	"""
296	Compute the Gamma-spread (Maximum Spread) indicator of the Pareto front.
297	It measures the maximum Euclidean distance between any two points in the
298	Pareto front.
299
300	Parameters
301	----------
302	pf: np.ndarray of shape (num_points, num_objectives)
303	Pareto front.
304
305	Returns
306	-------
307	float
308	Maximum spread indicator value.
309	"""
NEW 310	if len(pf) <= 1:	×
NEW 311	return np.nan	×
312
NEW 313	from scipy.spatial.distance import pdist	×
314
NEW 315	distances = pdist(pf)	×
NEW 316	if len(distances) == 0:	×
NEW 317	return np.nan	×
NEW 318	return float(np.max(distances))	×

SMTorg / smt-optim / 28718428290

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