mamut_routing_tools.campaign.select

Pick the subset that spans the feature space, subject to coverage quotas.

Two things have to be true of a benchmark that claims diversity, and they pull against each other:

Spread -- no two instances should be near-duplicates in feature space, and the extremes should be represented. That is a max-min (farthest-point) criterion on the z-scored descriptors.

Coverage -- every level of every design axis must appear often enough to be analysable. Seven demand distributions are useless if one of them shows up twice. Max-min alone will not deliver that: it chases outliers, and a whole level can sit in the middle of the cloud and never get picked.

So the greedy runs under residual-feasibility guards. Before each pick it asks, per quota, whether the remaining budget can still cover the outstanding minimums; when a quota becomes binding, the pick is restricted to candidates that reduce its deficit. A swap pass then improves the spread without letting any quota slip.

QuotaInfeasibleError

Bases: ValueError

The pool cannot satisfy the quotas within the requested budget.

Quota dataclass

Quota(name: str, key: Callable[[EvaluatedCandidate], Hashable], minimum: Mapping[Hashable, int] = (), maximum: Mapping[Hashable, int] | None = None, per_level_maximum: int | None = None)

A coverage requirement over one design axis.

key maps a candidate to the level it occupies (its demand type, its city, its size bucket). minimum is how many of each level the selection must end up with; maximum caps a level -- {city: 1} style caps are written with :meth:per_level_maximum.

select_maxmin

select_maxmin(candidates: Sequence[EvaluatedCandidate], k: int, quotas: Sequence[Quota] = (), *, swap_passes: int = 2) -> list[EvaluatedCandidate]

k candidates spread over feature space and honouring every quota.

Source code in MAMUT-routing-tools/src/mamut_routing_tools/campaign/select.py

def select_maxmin(
    candidates: Sequence[EvaluatedCandidate],
    k: int,
    quotas: Sequence[Quota] = (),
    *,
    swap_passes: int = 2,
) -> list[EvaluatedCandidate]:
    """``k`` candidates spread over feature space and honouring every quota."""
    if k <= 0:
        raise ValueError("k must be positive")
    if len(candidates) < k:
        raise QuotaInfeasibleError(
            f"pool has {len(candidates)} usable candidates, fewer than the requested {k}"
        )

    features = _feature_matrix(candidates)
    distances = np.linalg.norm(features[:, None, :] - features[None, :, :], axis=2)
    order = sorted(range(len(candidates)), key=lambda i: _tie_break(candidates[i]))
    counts: dict[str, Counter] = {quota.name: Counter() for quota in quotas}
    chosen: list[int] = []
    chosen_set: set[int] = set()

    def level(quota: Quota, index: int) -> Hashable:
        return quota.key(candidates[index])

    def admissible(index: int) -> bool:
        return all(
            (cap := quota.cap_for(level(quota, index))) is None
            or counts[quota.name][level(quota, index)] < cap
            for quota in quotas
        )

    def deficits_reduced(index: int) -> int:
        """How many outstanding quota minimums taking ``index`` would advance."""
        return sum(
            1
            for quota in quotas
            if level(quota, index) in quota.deficit(counts[quota.name])
        )

    def outstanding() -> str:
        return ", ".join(
            f"{quota.name}={quota.deficit(counts[quota.name])}"
            for quota in quotas
            if quota.deficit(counts[quota.name])
        )

    def keeps_quotas_reachable(index: int) -> bool:
        """Would taking ``index`` leave every quota still satisfiable?

        A one-step lookahead rather than a "is this quota binding right now"
        test: it is what stops the greedy walking into a corner where the
        remaining budget cannot cover what the remaining candidates offer.
        """
        trial = {name: counter.copy() for name, counter in counts.items()}
        for quota in quotas:
            trial[quota.name][level(quota, index)] += 1
        return (
            _shortfall(quotas, candidates, trial, chosen_set | {index}, k - len(chosen) - 1)
            is None
        )

    def take(index: int) -> None:
        chosen.append(index)
        chosen_set.add(index)
        for quota in quotas:
            counts[quota.name][level(quota, index)] += 1

    # Fail before doing any work when the pool simply cannot meet the quotas.
    upfront = _shortfall(quotas, candidates, counts, chosen_set, k)
    if upfront is not None:
        raise QuotaInfeasibleError(upfront)

    # Seed with the most extreme admissible candidate: farthest from the pool's
    # centre of mass, so the first pick is an edge of the cloud rather than an
    # arbitrary point in the middle of it.
    from_centre = np.linalg.norm(features - features.mean(axis=0), axis=1)
    start = [i for i in order if admissible(i) and keeps_quotas_reachable(i)]
    if not start:
        raise QuotaInfeasibleError("no first pick leaves the quotas satisfiable")
    take(max(start, key=lambda i: (from_centre[i], _tie_break(candidates[i]))))

    while len(chosen) < k:
        open_picks = [i for i in order if i not in chosen_set and admissible(i)]
        if not open_picks:
            capped = ", ".join(
                quota.name for quota in quotas if quota.per_level_maximum or quota.maximum
            )
            raise QuotaInfeasibleError(
                f"only {len(chosen)} of {k} instances could be picked: every remaining "
                f"candidate is blocked by a per-level cap ({capped or 'none declared'})"
            )
        allowed = [i for i in open_picks if keeps_quotas_reachable(i)]
        if not allowed:
            # One-step lookahead cannot see joint tightness across several
            # quotas at once, so when it finds nothing safe, keep going with
            # whatever still advances a minimum.
            allowed = [i for i in open_picks if deficits_reduced(i) > 0]
        if not allowed:
            reason = _shortfall(quotas, candidates, counts, chosen_set, k - len(chosen))
            raise QuotaInfeasibleError(
                reason
                or f"no candidate advances the outstanding quotas after {len(chosen)} "
                f"pick(s): {outstanding()}"
            )

        # Coverage first, spread within it. While any minimum is outstanding the
        # primary key is how many deficits a pick advances; once they are all
        # met this is exactly pure max-min.
        #
        # The order matters and is not a detail. Chasing spread first and
        # repairing coverage at the end fails whenever the picks are already
        # nearly forced -- one instance per city over 102 cities with a budget of
        # 100 leaves the selector only *which* of a city's candidates to take,
        # and by the time the deficits are noticed the cities that could have
        # filled them are used up. Ties on deficit count are broken by max-min,
        # and early on most candidates advance the same number of quotas, so
        # spread still drives the choice.
        covering = any(quota.deficit(counts[quota.name]) for quota in quotas)
        nearest = distances[np.ix_(allowed, chosen)].min(axis=1)
        best = max(
            range(len(allowed)),
            key=lambda position: (
                deficits_reduced(allowed[position]) if covering else 0,
                nearest[position],
                _tie_break(candidates[allowed[position]]),
            ),
        )
        take(allowed[best])

    unmet = [quota for quota in quotas if quota.deficit(counts[quota.name])]
    if unmet:
        raise QuotaInfeasibleError(
            "the budget could not cover every minimum: "
            + ", ".join(f"{q.name}={q.deficit(counts[q.name])}" for q in unmet)
        )

    chosen = _improve_by_swapping(chosen, candidates, distances, quotas, swap_passes)
    return [candidates[index] for index in chosen]

diversity_report

diversity_report(selected: Sequence[EvaluatedCandidate], pool: Sequence[EvaluatedCandidate], quotas: Sequence[Quota] = ()) -> dict[str, Any]

Coverage, spread and confounding of a selected set, as publishable data.

feature_correlation is the part worth reading first. If two descriptors are strongly correlated across the selected set, the benchmark cannot separate their effects on a solver, however wide its spread looks.

Source code in MAMUT-routing-tools/src/mamut_routing_tools/campaign/select.py

def diversity_report(
    selected: Sequence[EvaluatedCandidate],
    pool: Sequence[EvaluatedCandidate],
    quotas: Sequence[Quota] = (),
) -> dict[str, Any]:
    """Coverage, spread and confounding of a selected set, as publishable data.

    ``feature_correlation`` is the part worth reading first. If two descriptors
    are strongly correlated *across the selected set*, the benchmark cannot
    separate their effects on a solver, however wide its spread looks.
    """
    if not selected:
        return {"selected": 0, "pool": len(pool)}

    # Descriptor units, not the selector's transformed space: this block is read
    # by people, and a route size reported as 1.097 when the instance holds three
    # customers per route is a misreport, not a rounding.
    selected_features = np.asarray(
        [descriptor_values(c.descriptors) for c in selected], dtype=float
    )
    pool_features = np.asarray(
        [descriptor_values(c.descriptors) for c in pool], dtype=float
    )

    features: dict[str, Any] = {}
    for position, name in enumerate(SELECTION_FEATURES):
        column, reference = selected_features[:, position], pool_features[:, position]
        pool_range = float(reference.max() - reference.min())
        selected_range = float(column.max() - column.min())
        features[name] = {
            "min": float(column.min()),
            "median": float(np.median(column)),
            "max": float(column.max()),
            "std": float(column.std()),
            # How much of what the pool offered on this axis the selection kept.
            "range_fraction_of_pool": selected_range / pool_range if pool_range > 0 else None,
        }

    # A feature that happens to be constant across the selection has no
    # correlation with anything; np.corrcoef would divide by its zero spread and
    # fill the row with NaN and a warning. Say "undefined" deliberately instead.
    varying = selected_features.std(axis=0) > 0.0
    correlation = np.full((len(SELECTION_FEATURES), len(SELECTION_FEATURES)), np.nan)
    if varying.sum() >= 2:
        block = np.corrcoef(selected_features[:, varying], rowvar=False)
        correlation[np.ix_(varying, varying)] = block
    worst_pair, worst_value = None, 0.0
    for i in range(len(SELECTION_FEATURES)):
        for j in range(i + 1, len(SELECTION_FEATURES)):
            value = correlation[i, j]
            if np.isfinite(value) and abs(value) > abs(worst_value):
                worst_pair = (SELECTION_FEATURES[i], SELECTION_FEATURES[j])
                worst_value = float(value)

    coverage = {
        quota.name: {
            str(level): count
            for level, count in sorted(
                Counter(quota.key(c) for c in selected).items(), key=lambda kv: str(kv[0])
            )
        }
        for quota in quotas
    }
    unmet = {
        quota.name: {str(level): short for level, short in quota.deficit(
            Counter(quota.key(c) for c in selected)
        ).items()}
        for quota in quotas
        if quota.deficit(Counter(quota.key(c) for c in selected))
    }

    standardized = _feature_matrix(list(pool))
    index_of = {id(c): i for i, c in enumerate(pool)}
    chosen = [index_of[id(c)] for c in selected if id(c) in index_of]
    spread = _spread_score(
        np.linalg.norm(standardized[:, None, :] - standardized[None, :, :], axis=2), chosen
    ) if len(chosen) == len(selected) else ()

    return {
        "selected": len(selected),
        "pool": len(pool),
        "coverage": coverage,
        "unmet_quotas": unmet,
        "features": features,
        "feature_correlation": {
            "names": list(SELECTION_FEATURES),
            # JSON has no NaN; an undefined correlation is null.
            "matrix": [
                [float(v) if np.isfinite(v) else None for v in row] for row in correlation
            ],
            "strongest_pair": list(worst_pair) if worst_pair else None,
            "strongest_value": worst_value,
        },
        "min_pairwise_distance": float(spread[0]) if spread else None,
    }