Criteria#

mergen.criteria#

Optimisation criteria for space-filling design.

Each criterion exposes four methods:

evaluate(X, space)

Full computation on a normalised design matrix. Used once per algorithm restart to obtain the initial score.

incremental(X, i, new_pt, space, current_score)

O(n*d) update when the entire point i is replaced by new_pt. Used by single-point swap operators.

begin_1d(X, i, current_score)

Initialise a per-point cache for fast 1D coordinate-exchange trials.

try_1d(cache, axis, new_value, space)

O(n) update when only one coordinate changes. Inner-loop operation of the SCE algorithm.

All criteria operate on normalised coordinates in [0, 1]^d.

Available criteria#

‘umaxpro’ UMaxPro — Vorechovsky & Masek (2026) (default) ‘maxpro’ MaxPro — Joseph, Gul & Ba (2015) ‘phi_p’ PhiP (p=15) — Morris & Mitchell (1995) ‘cd2’ CD2 — Hickernell (1998) ‘stratified’ StratifiedL2 — Tian & Xu (2025)

Usage#

from mergen.criteria import get_criterion

crit = get_criterion(‘umaxpro’) score = crit.evaluate(X_norm, space)

class mergen.criteria.BaseCriterion[source]#

Bases: ABC

Abstract base for space-filling design criteria.

Subclasses implement evaluate() and incremental(). Fast 1D coordinate-exchange updates (begin_1d(), try_1d()) are provided in the base class as a fallback and may be overridden by criteria with separable structure for an O(n) inner loop.

Class attributes#

supports_nominalbool

Whether the criterion is mathematically well-defined on parameter spaces containing nominal (unordered categorical) factors. False for distance- and discrepancy-based criteria whose metric assumes ordered levels; Sampler refuses to run such a criterion when the space has any nominal factor. Subclasses designed for mixed factor spaces (e.g. MaxProQQ, QQD) set this to True.

supports_nominal: bool = False#
abstractmethod evaluate(X, space)[source]#

Compute the criterion value for design X.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in [0, 1]^d)

  • space (ParameterSpace — used for dimension metadata if needed)

Returns:

float — raw criterion score (not log-transformed)

Return type:

float

abstractmethod incremental(X, i, new_pt, space, current_score)[source]#

O(n·d) incremental update when point i is swapped for new_pt.

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point being swapped)

  • new_pt (np.ndarray, shape (d,) — candidate replacement (normalised))

  • space (ParameterSpace)

  • current_score (float — raw criterion score before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score). Negative → improvement; positive → degradation.

  • new_score (float) – Raw criterion score after the swap.

Return type:

Tuple[float, float]

begin_1d(X, i, current_score)[source]#

Initialise a per-point cache for fast 1D coordinate-exchange updates.

Called once before trying multiple coordinate changes on the same point i. The returned cache object is then passed to try_1d() for each candidate value.

Subclasses with separable structures (UMaxPro, MaxPro, PhiP) override this to precompute per-pair contributions, enabling O(n) per-try updates instead of O(n·d).

Default implementation stores (X, i, current_score) for the base try_1d() fallback, which builds the full new point and delegates to incremental().

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point to be modified)

  • current_score (float — raw criterion score before any change)

Returns:

cache (object) – Opaque structure used by try_1d(). Subclass-specific.

try_1d(cache, axis, new_value, space)[source]#

Try changing one coordinate of the cached point and return the delta.

For separable criteria, this runs in O(n). For non-separable criteria (CD2, StratifiedL2), the base implementation falls back to the full O(n·d) incremental() evaluation.

The returned cache reflects the trial state and may be reused for further trials on the same point. It does not modify the original design array.

Parameters:

  • cache (object — returned by begin_1d())

  • axis (int — coordinate index being changed (0 ≤ axis < d))

  • new_value (float — new normalised value at coordinate axis)

  • space (ParameterSpace)

Returns:

  • log_delta (float — log(new_score) − log(current_score))

  • new_score (float — raw criterion score after the change)

  • new_cache (object — cache updated to reflect the trial change)

Return type:

Tuple[float, float, object]

class mergen.criteria.UMaxPro[source]#

Bases: BaseCriterion

Uniform MaxPro criterion with periodic (toroidal) distance.

Defines the per-axis squared periodic distance:

δ²_v(i, j) = min(|x_iv − x_jv|,  1 − |x_iv − x_jv|)²

The criterion to minimise is then:

uMaxPro(X) = Σ_{i<j}  1 / prod_v  δ²_v(i, j)

The toroidal wrap removes the boundary bias of plain MaxPro: corner regions are no longer artificially far from each other, so the resulting design is statistically uniform across the entire hypercube including its boundaries.

References

Vorechovsky, M. & Masek, J. (2026).

Computers & Structures. [uMaxPro]

Joseph, V. R., Gul, E. & Ba, S. (2015).

Biometrika, 102(2), 371–380. [original MaxPro]

evaluate(X, space)[source]#

Compute the criterion value for design X.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in [0, 1]^d)

  • space (ParameterSpace — used for dimension metadata if needed)

Returns:

float — raw criterion score (not log-transformed)

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n·d) incremental update when point i is swapped for new_pt.

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point being swapped)

  • new_pt (np.ndarray, shape (d,) — candidate replacement (normalised))

  • space (ParameterSpace)

  • current_score (float — raw criterion score before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score). Negative → improvement; positive → degradation.

  • new_score (float) – Raw criterion score after the swap.

begin_1d(X, i, current_score)[source]#

Precompute per-pair contributions and per-axis squared periodic distances for fast 1D coordinate-exchange trials.

Complexity#

O(n · d) once per point. Each subsequent try_1d() call then costs O(n), so the total cost for trying all d axes on the same point is O(n · d) — the same as a single full incremental(), but split across d independent trials.

try_1d(cache, axis, new_value, space)[source]#

O(n) update: change one coordinate and rescale the cached contributions by the ratio of squared periodic distances.

Mathematics#

For each j ≠ i, the pairwise contribution to uMaxPro is

c_j = 1 / prod_v δ_v(i, j)²

Changing only coordinate axis of point i replaces a single factor in the product:

c_j_new = c_j_old · ( δ_axis_old² / δ_axis_new² )

and the total score updates as

S_new = S_old − Σ_j c_j_old + Σ_j c_j_new.

class mergen.criteria.MaxPro[source]#

Bases: BaseCriterion

MaxPro criterion with Euclidean squared distance.

Criterion (minimise):

MaxPro = Σ_{i<j}  1 / prod_v  (x_iv - x_jv)²

Notes

Suffers from boundary bias: corner points are treated as maximally separated even when they are geometrically close under a periodic metric. Prefer UMaxPro for general use.

References

Joseph, V. R., Gul, E. & Ba, S. (2015).

Biometrika, 102(2), 371–380.

evaluate(X, space)[source]#

Compute the criterion value for design X.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in [0, 1]^d)

  • space (ParameterSpace — used for dimension metadata if needed)

Returns:

float — raw criterion score (not log-transformed)

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n·d) incremental update when point i is swapped for new_pt.

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point being swapped)

  • new_pt (np.ndarray, shape (d,) — candidate replacement (normalised))

  • space (ParameterSpace)

  • current_score (float — raw criterion score before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score). Negative → improvement; positive → degradation.

  • new_score (float) – Raw criterion score after the swap.

begin_1d(X, i, current_score)[source]#

Precompute per-pair contributions and per-axis squared distances.

try_1d(cache, axis, new_value, space)[source]#

O(n) rescaling update for MaxPro (Euclidean variant).

See UMaxPro.try_1d() for the algebra; the only difference is that the per-axis squared distance is the plain Euclidean (u - v)² rather than the periodic minimum form.

class mergen.criteria.MaxProQQ(delta=0.0)[source]#

Bases: BaseCriterion

Maximum Projection criterion for quantitative and qualitative factors (Joseph, Gul & Ba 2019).

Extends MaxPro to spaces that mix continuous, discrete-numeric / integer / ordinal, and nominal factors. When the space contains only continuous columns, the criterion reduces exactly to plain MaxPro.

Parameters:

delta (float, default 0.0) – Additive floor for the pairwise denominator. Guards against division by zero when two rows share the same value on a continuous column (which does not happen for a valid LHD but can arise in constrained or coarse-grid searches). Joseph et al. (2019) do not use a floor; the default of 0.0 reproduces the paper exactly. Set to a small positive value (e.g. 1e-12) for extra numerical safety.

Variables:

supports_nominal (bool) – True — the criterion is defined on spaces containing nominal factors.

Notes

The design matrix passed to evaluate() and incremental() is the same normalised X that every other Mergen criterion sees; nominal columns arrive as their integer level indices divided by L_h - 1 (so distinct levels stay distinct after normalisation). The indicator \(\mathbb{1}(v_i \neq v_j)\) is therefore evaluated by an np.isclose comparison rather than by looking at the raw labels.

References

Joseph, V. R., Gul, E. & Ba, S. (2020).

Journal of Quality Technology, 52(4), 343-354.

supports_nominal: bool = True#
evaluate(X, space)[source]#

Compute the MaxProQQ score on the normalised design X.

To stay consistent with Mergen’s MaxPro — which reports the inner pairwise sum rather than the outer \((\cdot)^{1/d}\) root — this method also returns the raw sum

\[\text{score}(D) = \sum_{i<j} \frac{1}{\prod_{k} c_k(x_{ik}, x_{jk}) + \delta},\]

where the per-column term \(c_k\) follows Joseph, Gul & Ba (2019) Eq. 9. The Mergen optimisers rank designs by log(score), so the monotone \(1/d\) power is redundant and would only shrink dynamic range for the SA/SCE acceptance rules. Users who want the paper’s \(\psi(D)\) can recover it as (score / (n*(n-1)/2)) ** (1/d).

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in) – \([0, 1]^d\).

  • space (ParameterSpace — supplies the per-column factor type) – and level count.

Returns:

float — the raw pairwise score (smaller is better).

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n*d) update when point i is replaced by new_pt.

The MaxProQQ objective is a sum over pairs, so replacing a single row changes only the n-1 pairs that involve that row. We subtract those pairs’ old contributions and add the new ones without recomputing the full \(O(n^2)\) sum.

Parameters:

  • X (np.ndarray, shape (n, d))

  • i (int — row being replaced)

  • new_pt (np.ndarray, shape (d,))

  • space (ParameterSpace)

  • current_score (float — raw pairwise sum before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score); negative values mark an improvement.

  • new_score (float) – Raw pairwise sum after the swap.

Return type:

Tuple[float, float]

class mergen.criteria.PhiP(p=15)[source]#

Bases: BaseCriterion

Φ_p (phi-p) space-filling criterion.

Criterion (minimise):

Φ_p(X) = ( Σ_{i<j}  ||x_i − x_j||^{−p} )^{1/p}

As p → ∞ this converges to the maximin criterion (maximise the minimum pairwise distance). p = 15 is the standard default introduced by Morris & Mitchell (1995) as a smooth, differentiable maximin proxy that is well behaved under local search.

Parameters:

p (int or float, default 15)

References

Morris, M. D. & Mitchell, T. J. (1995).

Journal of Statistical Planning and Inference, 43, 381–402.

evaluate(X, space)[source]#

Compute the criterion value for design X.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in [0, 1]^d)

  • space (ParameterSpace — used for dimension metadata if needed)

Returns:

float — raw criterion score (not log-transformed)

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n·d) incremental update when point i is swapped for new_pt.

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point being swapped)

  • new_pt (np.ndarray, shape (d,) — candidate replacement (normalised))

  • space (ParameterSpace)

  • current_score (float — raw criterion score before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score). Negative → improvement; positive → degradation.

  • new_score (float) – Raw criterion score after the swap.

begin_1d(X, i, current_score)[source]#

Precompute per-pair squared distances (per axis and total) plus the per-pair contribution d^(-p), supporting O(n) updates when only one coordinate of point i changes.

try_1d(cache, axis, new_value, space)[source]#

O(n) update for Φ_p when only one coordinate of point i changes.

Mathematics#

The squared distance to each neighbour updates as

d²_new(i, j) = d²_old(i, j) − δ_axis_old² + δ_axis_new²,

and the per-pair contribution becomes

c_j_new = max(√d²_new, ε)^(−p).

The total inner sum is then

S_new = S_old − Σ_j c_j_old + Σ_j c_j_new,

from which the new Φ_p score is S_new^(1/p).

class mergen.criteria.CD2[source]#

Bases: BaseCriterion

Centred L2 discrepancy.

Measures deviation of the empirical distribution from the uniform distribution on [0, 1]^d. Lower values indicate better uniformity.

Criterion (minimise):

CD2 = sqrt(
    (13/12)^d
    - (2/n) Σ_i  prod_v (1 + 0.5|x_iv - 0.5| - 0.5(x_iv - 0.5)²)
    + (1/n²) Σ_{i,k} prod_v (1 + 0.5|x_iv - 0.5|
                                + 0.5|x_kv - 0.5|
                                - 0.5|x_iv - x_kv|)
)

References

Hickernell, F. J. (1998).

Mathematics of Computation, 67(221), 299–322.

evaluate(X, space)[source]#

Compute the criterion value for design X.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in [0, 1]^d)

  • space (ParameterSpace — used for dimension metadata if needed)

Returns:

float — raw criterion score (not log-transformed)

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n·d) incremental update when point i is swapped for new_pt.

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point being swapped)

  • new_pt (np.ndarray, shape (d,) — candidate replacement (normalised))

  • space (ParameterSpace)

  • current_score (float — raw criterion score before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score). Negative → improvement; positive → degradation.

  • new_score (float) – Raw criterion score after the swap.

class mergen.criteria.QQD[source]#

Bases: BaseCriterion

Qualitative-Quantitative Discrepancy (Zhang, Yang & Zhou 2021).

Discrepancy-style companion to MaxProQQ for designs that mix quantitative (continuous, discrete, integer, ordinal) and nominal factors. QQD generalises Hickernell’s wrap-around L2 discrepancy: on all-quantitative spaces it reduces to the WD of the quantitative sub-design.

Variables:

supports_nominal (bool) – True — the criterion is well-defined on spaces containing nominal factors and is one of the two entry points (with MaxProQQ) that Mergen exposes for such spaces.

Notes

evaluate() returns \(\text{QQD}^{2}(D)\) (the value reported throughout Zhang, Yang & Zhou 2021). The discrepancy is strictly positive on any non-degenerate design; a small floor at _EPS guards log(score) calls inside the SA / SCE acceptance rules against harmless numerical drift below zero.

Compared with the pair-sum objective of MaxProQQ, QQD is naturally computed as an n x n matrix sum because the paper writes it over both symmetric and diagonal terms. The diagonal contribution is \(n \cdot (3/2)^{p+q}\) — independent of the design — so a design swap changes only the off-diagonal contributions of the affected row and column, giving an \(O(n\,d)\) incremental update.

References

Zhang, M., Yang, F. & Zhou, Y.-D. (2021). Statistics,

arXiv:2101.02416, Theorem 1.

supports_nominal: bool = True#
evaluate(X, space)[source]#

Compute \(\text{QQD}^{2}(D)\) on the normalised design.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in) – \([0, 1]^d\). Nominal columns carry their level indices divided by L_h - 1.

  • space (ParameterSpace — supplies the per-column factor type) – and level count.

Returns:

float — \(\text{QQD}^{2}(D)\), smaller is better.

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n*d) update when point i is replaced by new_pt.

The QQD closed form sums a symmetric n x n kernel matrix. Replacing row i changes only row i and column i; by symmetry those two contributions are equal and can be subtracted / added as 2 * row_i after excluding the (i, i) diagonal entry, which is design-independent.

Parameters:

  • X (np.ndarray, shape (n, d))

  • i (int — row being replaced)

  • new_pt (np.ndarray, shape (d,))

  • space (ParameterSpace)

  • current_score (float — \(\text{QQD}^{2}(D)\) before swap)

Returns:

  • log_delta (float — log(new_score) - log(current_score).)

  • new_score (float — \(\text{QQD}^{2}(D)\) after the swap.)

Return type:

Tuple[float, float]

class mergen.criteria.StratifiedL2(s=2, p=None, weights='auto')[source]#

Bases: BaseCriterion

Stratified L2-discrepancy (Tian & Xu 2025).

Measures uniformity by aggregating the L2-norm of the local discrepancy over a hierarchical family of stratified regions. For a given base s and maximum depth p, the unit hypercube \([0,1)^m\) is stratified into s^u equal intervals on each axis for \(u = 0, 1, \ldots, p\), producing \((p+1)^m\) possible stratifications. The criterion sums the squared L2 deviations from uniformity across all of these stratifications, weighted by w(u).

Closed-form computational expression (Theorem 1, Eq. 13):

\[SD(P)^2 = -\left[\sum_{i=0}^p w(i)\, s^{-2i}\right]^m + \frac{1}{n^2} \sum_{a,b=1}^n \prod_{j=1}^m \left[\sum_{i=0}^p w(i)\, s^{-i}\, \delta_i(x_{aj}, x_{bj})\right]\]

where \(\delta_i(t, z) = \mathbf{1}\{\lfloor s^i t \rfloor = \lfloor s^i z \rfloor\}\). Note that \(\delta_0(t, z) \equiv 1\).

Parameters:

  • s (int, default 2) – Base number of strata per dimension. Common choices are s = 2 (binary partitions) or s = 3 (ternary).

  • p (int or None, default None) – Maximum stratification depth. None selects floor(log_s(n)) at evaluation time so that s^p ≤ n (Tian & Xu’s recommendation to avoid over-stratification).

  • weights ({'constant', 'exponential'} or array-like, default 'auto') –

    • 'constant': w(i) = 1 for all i, suitable for low to moderate dimension.

    • 'exponential': w(i) = y^i with y = 2/(m+1), which bounds SD(P)^2 < e regardless of dimension (Corollary 1).

    • 'auto' (default): exponential when m ≥ 8 (curse-of- dimensionality territory), constant otherwise.

    • array-like: explicit weights [w(0), w(1), …, w(p)].

Notes

Computational complexity is \(O(n^2 m p)\) per evaluation. The base incremental() falls back to full re-evaluation, which keeps SA’s cost at \(O(n^2 m p)\) per move; this matches CD2 and is acceptable for n \le 200.

References

Tian, Y. & Xu, H. (2025). A stratified L2-discrepancy with

application to space-filling designs. Journal of the Royal Statistical Society, Series B, 88(2).

evaluate(X, space)[source]#

Compute the criterion value for design X.

Parameters:

  • X (np.ndarray, shape (n, d) — normalised coordinates in [0, 1]^d)

  • space (ParameterSpace — used for dimension metadata if needed)

Returns:

float — raw criterion score (not log-transformed)

Return type:

float

incremental(X, i, new_pt, space, current_score)[source]#

O(n·d) incremental update when point i is swapped for new_pt.

Parameters:

  • X (np.ndarray, shape (n, d) — current normalised design)

  • i (int — index of the point being swapped)

  • new_pt (np.ndarray, shape (d,) — candidate replacement (normalised))

  • space (ParameterSpace)

  • current_score (float — raw criterion score before the swap)

Returns:

  • log_delta (float) – log(new_score) - log(current_score). Negative → improvement; positive → degradation.

  • new_score (float) – Raw criterion score after the swap.

mergen.criteria.get_criterion(name)[source]#

Instantiate a criterion by name.

Parameters:

name (str) – One of: 'umaxpro', 'maxpro', 'phi_p', 'cd2', 'stratified' (or 'stratified_l2'). Case-insensitive.

Returns:

BaseCriterion instance

Raises:

ValueError if name is not recognised. –

Return type:

BaseCriterion

Examples

>>> crit = get_criterion('umaxpro')
>>> type(crit).__name__
'UMaxPro'
mergen.criteria.list_criteria()[source]#

Return the list of available criterion names (canonical, no aliases).

Return type:

list

mergen.criteria.criterion_latex(name)[source]#

Return the LaTeX math label for a criterion, or the name itself.

Aliases resolve to their canonical label. Plain-text output targets (terminal, CSV) should use the raw name instead of this label.

Parameters:

name (str)

Return type:

str