"""Money Pump Index (MPI) computation for measuring exploitable inconsistency."""
from __future__ import annotations
import time
from typing import Any, cast
import numpy as np
from numpy.typing import NDArray
from prefgraph.core.session import ConsumerSession
from prefgraph.core.result import MPIResult, MPIBoundsResult, HoutmanMaksResult
from prefgraph.core.types import Cycle
# Hyperparameter: ILP vs greedy threshold for Houtman-Maks.
# The exact ILP (scipy.optimize.milp / HiGHS) is the only trustworthy method.
# The greedy fallback runs feedback-vertex-set on the transitive closure, which
# is a complete bidirectional digraph, so it keeps one observation per cyclic
# component and OVER-removes (1.5 to 3 times the optimum) at large T. It is a
# grossly loose upper bound on removals (at T=100 it removed 97 observations
# where the exact ILP removes 3), used only above this threshold where the exact
# MILP is intractable. Raised from 50 so per-user budget data up to T=100 uses
# the exact value; the exact ILP costs a few seconds at T=100 but correctness
# outweighs the greedy's wild over-removal.
HOUTMAN_MAKS_ILP_THRESHOLD = 100
def compute_mpi(
session: ConsumerSession,
tolerance: float = 1e-10,
method: str = "cycles",
) -> MPIResult:
"""
Compute Money Pump Index for the consumer data.
The MPI measures the percentage of total expenditure that could be
"pumped" from a consumer exhibiting cyclic preferences by an arbitrager.
For a violation cycle k1 -> k2 -> ... -> kn -> k1:
MPI = sum(p_ki @ x_ki - p_ki @ x_{ki+1}) / sum(p_ki @ x_ki)
Interpretation:
- MPI = 0.0: Consistent behavior (no money can be pumped)
- MPI = 0.10: 10% of budget could be extracted
- MPI = 1.0: Complete irrational behavior (theoretical maximum)
Args:
session: ConsumerSession with prices and quantities
tolerance: Numerical tolerance for GARP detection
method: Retained for backward compatibility. The headline MPI is now
always the exact maximum money-pump fraction (minimum
cost-to-budget cycle ratio) regardless of this value.
Returns:
MPIResult with MPI value, worst cycle, and all cycle costs
Example:
>>> import numpy as np
>>> from prefgraph import ConsumerSession, compute_mpi
>>> # Data with GARP violation
>>> prices = np.array([[1.0, 1.0], [1.0, 1.0]])
>>> quantities = np.array([[3.0, 1.0], [1.0, 3.0]])
>>> session = ConsumerSession(prices=prices, quantities=quantities)
>>> result = compute_mpi(session)
>>> print(f"MPI: {result.mpi_value:.4f}")
"""
start_time = time.perf_counter()
# --------------------------------------------------------------------------
# The Money Pump Index is the maximum, over revealed-preference cycles, of
# the money-pump cost as a fraction of the cycle's total budget:
#
# MPI = max_cycle sum_l (p^{k_l} . x^{k_l} - p^{k_l} . x^{k_{l+1}})
# / sum_l (p^{k_l} . x^{k_l})
#
# Echenique, Lee & Shum (2011) JPE 119(6) Eq. (2); Smeulders & Spieksma
# (2013) JPE 121(6) Theorem 2 show this maximum is the minimum cost-to-budget
# cycle ratio (Megiddo 1979), computed exactly by `compute_mpi_bounds`. The
# earlier code used Karp's minimum cycle MEAN (dividing by edge count, not
# summed budgets), which is a different objective and could report an MPI
# above the true maximum on heterogeneous-budget cycles.
# --------------------------------------------------------------------------
from prefgraph.algorithms.garp import check_garp
garp_result = check_garp(session, tolerance)
total_expenditure = float(session.own_expenditures.sum())
if garp_result.is_consistent:
computation_time = (time.perf_counter() - start_time) * 1000
return MPIResult(
mpi_value=0.0,
worst_cycle=None,
cycle_costs=[],
total_expenditure=total_expenditure,
computation_time_ms=computation_time,
)
# Headline value: the exact maximum money-pump fraction (minimum
# cost-to-budget cycle ratio). Use the fast Rust routine when available (it
# computes the same value after the objective fix) and fall back to the
# Python min-cycle-ratio otherwise. Both agree to binary-search precision.
from prefgraph._rust_backend import HAS_RUST, _rust_analyze_batch
mpi_val: float | None = None
if HAS_RUST:
try:
p = np.ascontiguousarray(session.prices, dtype=np.float64)
q = np.ascontiguousarray(session.quantities, dtype=np.float64)
# Args match engine.py _analyze_chunk_rust: (prices, quantities, ccei,
# mpi, harp, hm, utility, vei, vei_exact, network, tolerance).
results = cast(Any, _rust_analyze_batch)(
[p],
[q],
False,
True,
False,
False,
False,
False,
False,
False,
tolerance,
)
mpi_val = float(results[0]["mpi"])
except Exception:
mpi_val = None
if mpi_val is None:
mpi_val = compute_mpi_bounds(session, tolerance=tolerance).maximum_mpi
# Per-cycle breakdown (ratio-of-sums per GARP violation cycle) for
# worst_cycle and cycle_costs. These are a subset of all cycles, so the
# headline mpi_val can exceed the largest per-cycle value here.
E = session.expenditure_matrix
cycle_costs: list[tuple[Cycle, float]] = []
for cycle in garp_result.violations:
mc = _compute_cycle_mpi(cycle, E)
if mc > 0:
cycle_costs.append((cycle, mc))
worst_cycle: Cycle | None = None
if cycle_costs:
worst_cycle = max(cycle_costs, key=lambda x: x[1])[0]
elif garp_result.violations:
worst_cycle = garp_result.violations[0]
computation_time = (time.perf_counter() - start_time) * 1000
return MPIResult(
mpi_value=mpi_val,
worst_cycle=worst_cycle,
cycle_costs=cycle_costs,
total_expenditure=total_expenditure,
computation_time_ms=computation_time,
)
_RatioEdge = tuple[int, float, float]
[docs]
def compute_mpi_bounds(
session: ConsumerSession,
tolerance: float = 1e-10,
convergence_tolerance: float = 1e-10,
max_iterations: int = 100,
) -> MPIBoundsResult:
"""
Compute minimum and maximum Money Pump Index values.
This follows the cost/time graph construction used for the efficiently
computable MPI extremes: the minimum MPI is the minimum cost-to-budget
cycle ratio over revealed-preference violations, and the maximum MPI is
obtained by applying the same routine to negated savings.
Args:
session: ConsumerSession with prices and quantities
tolerance: Numerical tolerance for revealed-preference edges
convergence_tolerance: Binary-search tolerance for the cycle ratio
max_iterations: Maximum binary-search iterations
Returns:
MPIBoundsResult with minimum_mpi and maximum_mpi.
"""
start_time = time.perf_counter()
# prices/quantities are Optional on the dataclass but always populated after
# __post_init__; narrow to non-Optional for indexing below (no runtime effect).
prices = cast("NDArray[np.float64]", session.prices)
quantities = cast("NDArray[np.float64]", session.quantities)
expenditures = session.own_expenditures
total_expenditure = float(expenditures.sum())
positive_budget = expenditures > tolerance
if np.count_nonzero(positive_budget) < 2:
computation_time = (time.perf_counter() - start_time) * 1000
return MPIBoundsResult(
minimum_mpi=0.0,
maximum_mpi=0.0,
total_expenditure=total_expenditure,
computation_time_ms=computation_time,
)
prices = prices[positive_budget]
quantities = quantities[positive_budget]
expenditures = expenditures[positive_budget]
expenditure_matrix = prices @ quantities.T
min_graph = _build_mpi_ratio_graph(
expenditure_matrix,
quantities,
expenditures,
bound="minimum",
tolerance=tolerance,
)
max_graph = _build_mpi_ratio_graph(
expenditure_matrix,
quantities,
expenditures,
bound="maximum",
tolerance=tolerance,
)
minimum_mpi = _minimum_cost_time_ratio(
min_graph,
low=0.0,
high=1.0,
convergence_tolerance=convergence_tolerance,
max_iterations=max_iterations,
tolerance=tolerance,
)
maximum_ratio = _minimum_cost_time_ratio(
max_graph,
low=-1.0,
high=0.0,
convergence_tolerance=convergence_tolerance,
max_iterations=max_iterations,
tolerance=tolerance,
)
maximum_mpi = -maximum_ratio
minimum_mpi = float(np.clip(minimum_mpi, 0.0, 1.0))
maximum_mpi = float(np.clip(maximum_mpi, 0.0, 1.0))
if maximum_mpi <= tolerance:
minimum_mpi = 0.0
maximum_mpi = 0.0
elif minimum_mpi > maximum_mpi:
minimum_mpi = maximum_mpi
computation_time = (time.perf_counter() - start_time) * 1000
return MPIBoundsResult(
minimum_mpi=minimum_mpi,
maximum_mpi=maximum_mpi,
total_expenditure=total_expenditure,
computation_time_ms=computation_time,
)
def _build_mpi_ratio_graph(
expenditure_matrix: NDArray[np.float64],
quantities: NDArray[np.float64],
own_expenditures: NDArray[np.float64],
bound: str,
tolerance: float,
) -> list[list[_RatioEdge]]:
"""Build the cost/time graph used for MPI extrema."""
n_obs = expenditure_matrix.shape[0]
graph: list[list[_RatioEdge]] = [[] for _ in range(n_obs)]
for i in range(n_obs):
budget = float(own_expenditures[i])
if budget <= tolerance:
continue
for j in range(n_obs):
if i == j:
continue
if expenditure_matrix[i, j] > budget + tolerance:
continue
same_bundle = np.allclose(
quantities[i],
quantities[j],
rtol=tolerance,
atol=tolerance,
)
if bound == "minimum" and same_bundle:
continue
if bound == "minimum":
cost = budget - float(expenditure_matrix[i, j])
else:
cost = float(expenditure_matrix[i, j]) - budget
graph[i].append((j, cost, budget))
return graph
def _minimum_cost_time_ratio(
graph: list[list[_RatioEdge]],
low: float,
high: float,
convergence_tolerance: float,
max_iterations: int,
tolerance: float,
) -> float:
"""Binary search the minimum cycle cost/time ratio."""
if not _has_cycle(graph):
return 0.0
lo = float(low)
hi = float(high)
for _ in range(max_iterations):
if abs(hi - lo) <= convergence_tolerance:
break
mid = (lo + hi) / 2.0
reweighted = _reweight_graph(graph, mid)
has_negative_cycle, distances = _bellman_ford_negative_cycle(
reweighted, tolerance
)
if has_negative_cycle:
hi = mid
continue
zero_graph = _zero_residual_graph(reweighted, distances, tolerance)
if _has_cycle(zero_graph):
return mid
lo = mid
return (lo + hi) / 2.0
def _reweight_graph(
graph: list[list[_RatioEdge]],
ratio: float,
) -> list[list[_RatioEdge]]:
"""Apply cost <- cost - ratio * time to every edge."""
return [
[(target, cost - ratio * budget, budget) for target, cost, budget in edges]
for edges in graph
]
def _bellman_ford_negative_cycle(
graph: list[list[_RatioEdge]],
tolerance: float,
) -> tuple[bool, NDArray[np.float64]]:
"""Detect any negative cycle using a super-source initialization."""
n_obs = len(graph)
distances = np.zeros(n_obs, dtype=np.float64)
for _ in range(max(0, n_obs - 1)):
changed = False
for i, edges in enumerate(graph):
for j, cost, _ in edges:
candidate = distances[i] + cost
if candidate < distances[j] - tolerance:
distances[j] = candidate
changed = True
if not changed:
break
for i, edges in enumerate(graph):
for j, cost, _ in edges:
if distances[i] + cost < distances[j] - tolerance:
return True, distances
return False, distances
def _zero_residual_graph(
graph: list[list[_RatioEdge]],
distances: NDArray[np.float64],
tolerance: float,
) -> list[list[_RatioEdge]]:
"""Build the graph of zero reduced-cost edges."""
residual: list[list[_RatioEdge]] = [[] for _ in graph]
for i, edges in enumerate(graph):
for j, cost, budget in edges:
reduced_cost = cost + distances[i] - distances[j]
if abs(reduced_cost) <= max(tolerance, 1e-9):
residual[i].append((j, 0.0, budget))
return residual
def _has_cycle(graph: list[list[_RatioEdge]]) -> bool:
"""Return True if a directed graph contains a cycle."""
n_obs = len(graph)
indegree = [0] * n_obs
for edges in graph:
for j, _, _ in edges:
indegree[j] += 1
stack = [i for i, degree in enumerate(indegree) if degree == 0]
visited = 0
while stack:
node = stack.pop()
visited += 1
for j, _, _ in graph[node]:
indegree[j] -= 1
if indegree[j] == 0:
stack.append(j)
return visited < n_obs
def _karp_mpi(
E: NDArray[np.float64],
own_exp: NDArray[np.float64],
R: NDArray[np.bool_],
) -> tuple[float, Cycle | None]:
"""
Compute MPI using Karp's algorithm on negated weights.
MPI wants the MAX mean-weight cycle (worst exploitation). Karp's
algorithm finds MIN mean-weight. So we negate: w[i,j] = -(savings)
and the min of negated = negative of the max of original.
Args:
E: T x T expenditure matrix
own_exp: Own expenditures (diagonal of E)
R: Direct revealed preference matrix (adjacency)
Returns:
Tuple of (mpi_value, worst_cycle) or (0.0, None) if no cycle
"""
from prefgraph._kernels import karp_min_mean_cycle_numba
T = E.shape[0]
# Build NEGATED weight matrix so Karp's min-mean finds the max-mean cycle
# Original: w[i,j] = (own_exp[i] - E[i,j]) / own_exp[i] (money pump per step)
# Negated: w[i,j] = -(own_exp[i] - E[i,j]) / own_exp[i]
weights = np.full((T, T), 1e18, dtype=np.float64)
for i in range(T):
if own_exp[i] > 0:
for j in range(T):
if R[i, j] and i != j:
weights[i, j] = -(own_exp[i] - E[i, j]) / own_exp[i]
adjacency = np.ascontiguousarray(R, dtype=np.bool_)
np.fill_diagonal(adjacency, False)
mean_weight, cycle_arr = karp_min_mean_cycle_numba(
np.ascontiguousarray(weights, dtype=np.float64),
adjacency,
)
if cycle_arr[0] == -1 or mean_weight >= 1e17:
return 0.0, None
cycle = tuple(int(x) for x in cycle_arr)
# Negate back: MPI = -min_mean of negated weights = max_mean of original
mpi_val = max(0.0, -float(mean_weight))
return mpi_val, cycle
def _compute_cycle_mpi(
cycle: Cycle,
E: NDArray[np.float64],
) -> float:
"""
Compute MPI for a single cycle.
For cycle k1 -> k2 -> ... -> kn -> k1:
MPI = sum_{i=1}^{n}(E[ki, ki] - E[ki, k_{i+1}]) / sum_{i=1}^{n}(E[ki, ki])
The numerator is the total "savings" if the consumer had chosen the
next bundle in the cycle at each step. The denominator is total expenditure
in the cycle.
Args:
cycle: Tuple of observation indices forming the cycle
E: Expenditure matrix where E[i,j] = p_i @ q_j
Returns:
MPI value for this cycle (0 to 1)
"""
if len(cycle) < 2:
return 0.0
numerator = 0.0
denominator = 0.0
# cycle is (k1, k2, ..., kn, k1) where last element repeats first
for i in range(len(cycle) - 1):
ki = cycle[i]
ki_next = cycle[i + 1]
# E[ki, ki] - E[ki, ki_next] = savings from choosing ki_next instead
savings = E[ki, ki] - E[ki, ki_next]
numerator += savings
# E[ki, ki] = expenditure at observation ki
denominator += E[ki, ki]
if denominator <= 0:
return 0.0
mpi = numerator / denominator
# MPI should be non-negative; clamp to handle numerical issues
return max(0.0, mpi)
def _compute_simple_mpi(
session: ConsumerSession,
violations: list[Cycle],
) -> float:
"""
Compute a simple aggregate MPI measure.
This is a fallback when cycle-based MPI is not well-defined.
It computes the average "wasted" money across all violation pairs.
Args:
session: ConsumerSession
violations: List of violation cycles
Returns:
Simple MPI estimate
"""
if not violations:
return 0.0
E = session.expenditure_matrix
total_waste = 0.0
total_spend = 0.0
for cycle in violations:
for i in range(len(cycle) - 1):
ki = cycle[i]
ki_next = cycle[i + 1]
waste = E[ki, ki] - E[ki, ki_next]
if waste > 0:
total_waste += waste
total_spend += E[ki, ki]
if total_spend <= 0:
return 0.0
return total_waste / total_spend
def compute_houtman_maks_index(
session: ConsumerSession,
tolerance: float = 1e-10,
method: str = "auto",
) -> HoutmanMaksResult:
"""
Compute Houtman-Maks index: minimum observations to remove for consistency.
The Houtman-Maks index is the size of the smallest subset of observations
that, when removed, makes the remaining data satisfy GARP.
Supports two methods:
- "ilp": Exact solution via Integer Linear Programming (Big-M Afriat
formulation with scipy.optimize.milp). Optimal but slower for large T.
- "greedy": Fast SCC + greedy FVS approximation (2-approximation factor).
- "auto" (default): Uses "ilp" for T <= HOUTMAN_MAKS_ILP_THRESHOLD, else "greedy".
Args:
session: ConsumerSession
tolerance: Numerical tolerance
method: "auto", "ilp", or "greedy"
Returns:
HoutmanMaksResult with fraction and list of removed observation indices
"""
from prefgraph.algorithms.garp import check_garp
start_time = time.perf_counter()
T = session.num_observations
if T < 2:
computation_time = (time.perf_counter() - start_time) * 1000
return HoutmanMaksResult(
fraction=0.0,
removed_observations=[],
computation_time_ms=computation_time,
)
# Quick GARP check first
garp_result = check_garp(session, tolerance)
if garp_result.is_consistent:
computation_time = (time.perf_counter() - start_time) * 1000
return HoutmanMaksResult(
fraction=0.0,
removed_observations=[],
computation_time_ms=computation_time,
)
# Choose method
if method == "auto":
method = "ilp" if T <= HOUTMAN_MAKS_ILP_THRESHOLD else "greedy"
if method == "ilp":
removed = _houtman_maks_ilp(session, tolerance)
else:
removed = _houtman_maks_greedy(session, tolerance)
computation_time = (time.perf_counter() - start_time) * 1000
fraction = len(removed) / T
return HoutmanMaksResult(
fraction=fraction,
removed_observations=removed,
computation_time_ms=computation_time,
)
def _houtman_maks_ilp(
session: ConsumerSession,
tolerance: float,
) -> list[int]:
"""
Exact Houtman-Maks via Big-M Afriat ILP formulation.
Maximize sum(z_i) subject to Afriat's inequalities holding whenever
both observations i and j are kept (z_i = z_j = 1).
Variables: z_i (binary), U_i (utility), lambda_i (marginal utility)
Constraint: U_i - U_j - lambda_j*(E[j,i] - E[j,j]) <= M*(2 - z_i - z_j)
"""
from scipy.optimize import milp, LinearConstraint, Bounds
T = session.num_observations
E = session.expenditure_matrix
own_exp = session.own_expenditures
# Variables layout: [z_0..z_{T-1}, U_0..U_{T-1}, lambda_0..lambda_{T-1}]
n_vars = 3 * T
# Objective: maximize sum(z_i) = minimize -sum(z_i)
c = np.zeros(n_vars)
c[:T] = -1.0 # Minimize negative z = maximize z
# Variable bounds
max_exp = float(np.max(own_exp))
max_coeff = float(np.max(np.abs(E - own_exp[:, np.newaxis])))
U_max = max_exp * T
lambda_lb = 1e-3
lambda_ub = T * 1.0
# M must deactivate the constraint when at least one z=0.
# Max |LHS| = U_max + lambda_ub * max_coeff. M must exceed this.
M = U_max + lambda_ub * max_coeff + 1.0
# Build constraints: for each (i,j) pair with i != j:
# U_i - U_j - lambda_j * (E[j,i] - E[j,j]) + M*z_i + M*z_j <= 2*M
n_constraints = T * (T - 1)
A = np.zeros((n_constraints, n_vars))
b = np.full(n_constraints, 2.0 * M)
idx = 0
for i in range(T):
for j in range(T):
if i == j:
continue
A[idx, T + i] = 1.0 # U_i
A[idx, T + j] = -1.0 # -U_j
A[idx, 2 * T + j] = -(E[j, i] - own_exp[j]) # -lambda_j * coeff
A[idx, i] = M # M * z_i
A[idx, j] = M # M * z_j
b[idx] = 2.0 * M
idx += 1
constraints = LinearConstraint(A, ub=b)
lb = np.zeros(n_vars)
ub = np.full(n_vars, M)
ub[:T] = 1.0 # z_i in [0, 1]
ub[T : 2 * T] = U_max # U_i in [0, U_max]
ub[2 * T :] = lambda_ub # lambda_i in [lambda_lb, lambda_ub]
lb[2 * T :] = lambda_lb
bounds = Bounds(lb, ub)
# Integer constraints: z_i are binary
integrality = np.zeros(n_vars)
integrality[:T] = 1 # z variables are integer (binary with bounds [0,1])
try:
result = milp(
c,
constraints=constraints,
# scipy's milp accepts a float ndarray for integrality at runtime (it
# casts internally); the stub only permits integer arrays. Third-party
# stub strictness, values here are exactly 0.0/1.0.
integrality=integrality, # type: ignore[arg-type]
bounds=bounds,
)
if result.success:
z = result.x[:T]
removed = [i for i in range(T) if z[i] < 0.5]
return removed
except Exception:
pass
# Fallback to greedy if ILP fails
return _houtman_maks_greedy(session, tolerance)
def _houtman_maks_greedy(
session: ConsumerSession,
tolerance: float,
) -> list[int]:
"""Greedy FVS-based Houtman-Maks (fast approximation).
Houtman & Maks (1985) "Determining All Maximal Data Subsets Consistent
with Revealed Preference", Kwantitatieve Methoden 19, 89-104.
The HM index finds the largest subset of observations consistent with
GARP. Equivalently, it finds the minimum set of observations to remove
so that the remaining data satisfies GARP.
Heufer & Hjertstrand (2015) "Consistent Subsets", Section 2, define:
"The HM-index is the maximal fraction of non-zero elements in the
binary vector v such that GARP(v) holds."
GARP violation (Varian 1982, Def 2.1 in Smeulders et al. 2014):
"if q_i R q_j, then it is not the case that q_j P^0 q_i"
where R is the TRANSITIVE closure of R^0. Violations live in R_star,
not in the direct relation R.
Smeulders et al. (2014) Theorem 5.1 proves HI-GARP is NP-hard.
This greedy uses SCC decomposition + feedback vertex set (FVS) as
an approximation.
Key fix (2026-03-28): SCC decomposition must use R_star (transitive
closure), not R (direct relation). Purely transitive GARP violations
- where the violation path i ->* j goes through intermediates with no
direct back-edge - produce all-singleton SCCs in R. The greedy loop
then skips all components, returning removed=[] even when violations
exist. Using R_star ensures all observations participating in any
cycle (direct or transitive) are grouped into the same SCC.
"""
from prefgraph.graph.scc import find_sccs, greedy_feedback_vertex_set
from prefgraph.graph.transitive_closure import floyd_warshall_transitive_closure
E = session.expenditure_matrix
own_exp = session.own_expenditures
# Direct revealed-preference relations (Varian 1982):
# R[i,j] = True iff p_i · x_i >= p_i · x_j (i could afford j's bundle)
# P[i,j] = True iff p_i · x_i > p_i · x_j (strictly)
R = own_exp[:, np.newaxis] >= E - tolerance
P = own_exp[:, np.newaxis] > E + tolerance
np.fill_diagonal(P, False)
# R_star = transitive closure of R via Floyd-Warshall.
# GARP violation: R_star[i,j] AND P[j,i] - i.e., i transitively
# revealed-prefers j, yet j strictly prefers its own bundle over i's.
R_star = floyd_warshall_transitive_closure(R)
violation_matrix = R_star & P.T
if not np.any(violation_matrix):
return []
# SCC decomposition on R_star (transitive closure), NOT R (direct).
# In R_star, i and j are in the same SCC iff R_star[i,j] AND R_star[j,i],
# meaning they are mutually reachable through chains of revealed preference.
# This guarantees all observations involved in ANY violation cycle - even
# purely transitive ones - are grouped together for FVS removal.
n_comp, labels = find_sccs(R_star)
scc_sizes = np.bincount(labels, minlength=n_comp)
removed: list[int] = []
for c in range(n_comp):
if scc_sizes[c] <= 1:
continue
scc_nodes = np.where(labels == c)[0]
sub_violation = violation_matrix[np.ix_(scc_nodes, scc_nodes)]
if not np.any(sub_violation):
continue
# Greedy FVS on R_star within this SCC. Using R_star (not R) ensures
# the FVS algorithm sees all transitive edges as direct arcs, so it
# can identify which nodes to remove to break transitive cycles.
sub_R = R_star[np.ix_(scc_nodes, scc_nodes)].copy()
fvs_local = greedy_feedback_vertex_set(sub_R)
for local_idx in fvs_local:
removed.append(int(scc_nodes[local_idx]))
return removed
# =============================================================================
# TECH-FRIENDLY ALIASES
# =============================================================================
# compute_confusion_metric: Tech-friendly name for compute_mpi
compute_confusion_metric = compute_mpi
"""
Compute the confusion metric (how exploitable the user's decisions are).
This is the tech-friendly alias for compute_mpi (Money Pump Index).
The confusion metric measures how much value could be extracted from
a user making inconsistent decisions via preference cycling.
Example:
>>> from prefgraph import BehaviorLog, compute_confusion_metric
>>> result = compute_confusion_metric(user_log)
>>> if result.confusion_score > 0.15:
... alert_ux_team(user_id)
Returns:
ConfusionResult with confusion_score in [0, 1]
"""
# compute_minimal_outlier_fraction: Tech-friendly name for compute_houtman_maks_index
compute_minimal_outlier_fraction = compute_houtman_maks_index
"""
Compute the minimal fraction of observations to remove to achieve consistency.
Tech-friendly alias for compute_houtman_maks_index.
Returns the smallest fraction of user sessions that must be removed to
make the remaining behavior fully consistent. Useful for identifying
which specific transactions are problematic.
"""