"""GARP (Generalized Axiom of Revealed Preference) detection algorithm."""
from __future__ import annotations
import time
from typing import Any, Callable, cast
import numpy as np
from numpy.typing import NDArray
from prefgraph.core.session import ConsumerSession
from prefgraph.core.result import (
GARPResult,
WARPResult,
SwapsIndexResult,
ObservationContributionResult,
MinimumCostIndexResult,
)
from prefgraph.core.types import Cycle
from prefgraph.core.exceptions import ComputationalLimitError
from prefgraph.algorithms._budget_axioms import (
check_budget_axiom_at_efficiency,
validate_efficiency_level,
)
from prefgraph.graph.transitive_closure import floyd_warshall_transitive_closure
from prefgraph._kernels import bfs_find_path_numba, find_violation_pairs_numba
def check_garp(
session: ConsumerSession,
tolerance: float = 1e-10,
efficiency: float = 1.0,
) -> GARPResult:
"""
Check if consumer data satisfies GARP using Warshall's algorithm.
GARP (Generalized Axiom of Revealed Preference) states that revealed
preferences must be acyclic when considering both weak and strict
preferences. A violation occurs when there exists a cycle in the
transitive closure that includes at least one strict preference.
The algorithm:
1. Compute direct revealed preference matrix R:
R[i,j] = True iff p_i @ x_i >= p_i @ x_j
(bundle j was affordable when i was chosen, so i is weakly preferred)
2. Compute strict revealed preference matrix P:
P[i,j] = True iff p_i @ x_i > p_i @ x_j
(bundle j was strictly cheaper, so i is strictly preferred)
3. Compute transitive closure R* of R using Floyd-Warshall
4. Check for violations: GARP is violated if exists i,j such that
R*[i,j] = True AND P[j,i] = True
(i is transitively preferred to j, but j is strictly preferred to i)
Args:
session: ConsumerSession with prices and quantities
tolerance: Numerical tolerance for floating-point comparisons
efficiency: Afriat-style budget efficiency level in [0, 1]. The default
1.0 checks exact GARP.
Returns:
GARPResult with consistency flag, violation cycles, and matrices
Example:
>>> import numpy as np
>>> from prefgraph import ConsumerSession, check_garp
>>> # Consistent data: when A is cheap, buy more A; when B is cheap, buy more B
>>> prices = np.array([[1.0, 2.0], [2.0, 1.0]])
>>> quantities = np.array([[4.0, 1.0], [1.0, 4.0]])
>>> session = ConsumerSession(prices=prices, quantities=quantities)
>>> result = check_garp(session)
>>> result.is_consistent
True
"""
start_time = time.perf_counter()
efficiency = validate_efficiency_level(efficiency)
if efficiency != 1.0:
axiom_result = check_budget_axiom_at_efficiency(
session,
axiom="garp",
efficiency=efficiency,
tolerance=tolerance,
)
computation_time = (time.perf_counter() - start_time) * 1000
return GARPResult(
is_consistent=axiom_result.is_consistent,
violations=axiom_result.violations, # type: ignore[arg-type]
direct_revealed_preference=axiom_result.direct_revealed_preference,
# axiom="garp" always populates transitive_closure (only the WARP
# branch leaves it None), so narrow away the Optional here.
transitive_closure=cast(
"NDArray[np.bool_]", axiom_result.transitive_closure
),
strict_revealed_preference=axiom_result.strict_revealed_preference,
computation_time_ms=computation_time,
)
# Rust backend uses garp_check_with_closure (Varian 1982, O(T³)
# Floyd-Warshall) because we need R* for violation cycle extraction.
# The O(T²) SCC-only test (Talla Nobibon et al. 2015 JOTA) correctly
# determines is_consistent but does not populate R*.
from prefgraph._rust_backend import HAS_RUST, _rust_build_preference_graph
if HAS_RUST:
try:
# _rust_build_preference_graph is Optional[Callable] for the
# HAS_RUST=False fallback; under this guard it is always set.
build_graph = cast("Callable[..., Any]", _rust_build_preference_graph)
g = build_graph(
np.ascontiguousarray(session.prices, dtype=np.float64),
np.ascontiguousarray(session.quantities, dtype=np.float64),
tolerance,
)
R = g["r"].astype(bool)
P = g["p"].astype(bool)
R_star = g["r_star"].astype(bool)
violation_matrix = R_star & P.T
is_consistent = bool(g["is_garp"])
violations: list[Cycle] = []
if not is_consistent:
violations = _find_violation_cycles(R, P, R_star, violation_matrix)
computation_time = (time.perf_counter() - start_time) * 1000
return GARPResult(
is_consistent=is_consistent,
violations=violations,
direct_revealed_preference=R,
transitive_closure=R_star,
strict_revealed_preference=P,
computation_time_ms=computation_time,
)
except Exception:
pass # Fall through to Python
# Python fallback
E = session.expenditure_matrix # T x T
own_exp = session.own_expenditures # Shape: (T,)
R = own_exp[:, np.newaxis] >= E - tolerance
P = own_exp[:, np.newaxis] > E + tolerance
np.fill_diagonal(P, False)
R_star = floyd_warshall_transitive_closure(R)
violation_matrix = R_star & P.T
is_consistent = not np.any(violation_matrix)
# `violations` is already declared as list[Cycle] in the HAS_RUST branch above;
# re-annotating here would be a redefinition, so just reassign.
violations = []
if not is_consistent:
violations = _find_violation_cycles(R, P, R_star, violation_matrix)
computation_time = (time.perf_counter() - start_time) * 1000
return GARPResult(
is_consistent=is_consistent,
violations=violations,
direct_revealed_preference=R,
transitive_closure=R_star,
strict_revealed_preference=P,
computation_time_ms=computation_time,
)
_MAX_REPORTED_VIOLATIONS = 1000
def _find_violation_cycles(
R: NDArray[np.bool_],
P: NDArray[np.bool_],
R_star: NDArray[np.bool_],
violation_matrix: NDArray[np.bool_],
) -> list[Cycle]:
"""
Find cycles that violate GARP.
A violation cycle is a sequence i1 -> i2 -> ... -> in -> i1 where:
- Each consecutive pair is connected by revealed preference (R)
- At least one edge is strict preference (P)
Optimized to scope BFS to within SCCs and cap reported violations.
Args:
R: Direct revealed preference matrix
P: Strict revealed preference matrix
R_star: Transitive closure of R
violation_matrix: R_star & P.T (pre-computed)
Returns:
List of violation cycles as tuples of observation indices
"""
from prefgraph.graph.scc import find_sccs
violations: list[Cycle] = []
seen_cycles: set[frozenset[int]] = set()
T = R.shape[0]
# For small graphs, use the direct approach
if T < 10:
return _find_violation_cycles_direct(R, P, R_star, violation_matrix)
# Find SCCs - violations only exist within SCCs
n_comp, labels = find_sccs(R)
scc_sizes = np.bincount(labels, minlength=n_comp)
for c in range(n_comp):
if scc_sizes[c] <= 1:
continue
if len(violations) >= _MAX_REPORTED_VIOLATIONS:
break
scc_mask = labels == c
scc_nodes = np.where(scc_mask)[0]
# Check if this SCC has violations
sub_violation = violation_matrix[np.ix_(scc_nodes, scc_nodes)]
if not np.any(sub_violation):
continue
# Extract local subgraphs for BFS
sub_R = np.ascontiguousarray(R[np.ix_(scc_nodes, scc_nodes)], dtype=np.bool_)
# Find violation pairs within this SCC
sub_R_star = np.ascontiguousarray(
R_star[np.ix_(scc_nodes, scc_nodes)], dtype=np.bool_
)
sub_P = np.ascontiguousarray(P[np.ix_(scc_nodes, scc_nodes)], dtype=np.bool_)
sub_pairs = find_violation_pairs_numba(sub_R_star, sub_P)
for idx in range(sub_pairs.shape[0]):
if len(violations) >= _MAX_REPORTED_VIOLATIONS:
break
li, lj = int(sub_pairs[idx, 0]), int(sub_pairs[idx, 1])
# BFS within the small SCC subgraph
path = _reconstruct_path_bfs(sub_R, li, lj)
if path is not None:
# Map local indices back to global
global_cycle = tuple(int(scc_nodes[k]) for k in path)
cycle_set = frozenset(global_cycle[:-1])
if cycle_set not in seen_cycles:
seen_cycles.add(cycle_set)
violations.append(global_cycle)
return violations
def _find_violation_cycles_direct(
R: NDArray[np.bool_],
P: NDArray[np.bool_],
R_star: NDArray[np.bool_],
violation_matrix: NDArray[np.bool_],
) -> list[Cycle]:
"""Direct violation cycle finding for small graphs (no SCC overhead)."""
violations: list[Cycle] = []
seen_cycles: set[frozenset[int]] = set()
R_star_c = np.ascontiguousarray(R_star, dtype=np.bool_)
P_c = np.ascontiguousarray(P, dtype=np.bool_)
violation_pairs = find_violation_pairs_numba(R_star_c, P_c)
R_c = np.ascontiguousarray(R, dtype=np.bool_)
for idx in range(violation_pairs.shape[0]):
i, j = int(violation_pairs[idx, 0]), int(violation_pairs[idx, 1])
path = _reconstruct_path_bfs(R_c, i, j)
if path is not None:
cycle = tuple(path)
cycle_set = frozenset(cycle[:-1])
if cycle_set not in seen_cycles:
seen_cycles.add(cycle_set)
violations.append(cycle)
return violations
def _reconstruct_path_bfs(
R: NDArray[np.bool_],
start: int,
end: int,
) -> list[int] | None:
"""
Reconstruct shortest path from start to end using BFS on R.
Uses Numba JIT for fast path finding.
Args:
R: Direct revealed preference adjacency matrix (must be contiguous)
start: Starting node index
end: Ending node index
Returns:
List of node indices forming the path (ending with start to complete cycle),
or None if no path exists
"""
path_arr = bfs_find_path_numba(R, np.int64(start), np.int64(end))
if len(path_arr) == 0 or path_arr[0] == -1:
return None
return list(path_arr)
def check_warp(
session: ConsumerSession,
tolerance: float = 1e-10,
efficiency: float = 1.0,
) -> WARPResult:
"""
Check if consumer data satisfies WARP (Weak Axiom of Revealed Preference).
WARP is a weaker condition than GARP. It only checks for direct (length-2)
violations: if x_i R x_j, then NOT x_j P x_i.
Args:
session: ConsumerSession with prices and quantities
tolerance: Numerical tolerance for comparisons
efficiency: Afriat-style budget efficiency level in [0, 1]. The default
1.0 checks exact WARP.
Returns:
WARPResult with is_consistent flag and list of violating pairs
Example:
>>> import numpy as np
>>> from prefgraph import BehaviorLog, check_warp
>>> prices = np.array([[1.0, 2.0], [2.0, 1.0]])
>>> quantities = np.array([[4.0, 1.0], [1.0, 4.0]])
>>> session = BehaviorLog(cost_vectors=prices, action_vectors=quantities)
>>> result = check_warp(session)
>>> result.is_consistent
True
"""
start_time = time.perf_counter()
efficiency = validate_efficiency_level(efficiency)
if efficiency != 1.0:
axiom_result = check_budget_axiom_at_efficiency(
session,
axiom="warp",
efficiency=efficiency,
tolerance=tolerance,
)
computation_time = (time.perf_counter() - start_time) * 1000
return WARPResult(
is_consistent=axiom_result.is_consistent,
violations=axiom_result.violations, # type: ignore[arg-type]
computation_time_ms=computation_time,
)
E = session.expenditure_matrix
own_exp = session.own_expenditures
R = own_exp[:, np.newaxis] >= E - tolerance
P = own_exp[:, np.newaxis] > E + tolerance
np.fill_diagonal(P, False)
# WARP violation: R[i,j] AND P[j,i]
violation_matrix = R & P.T
violations = [
(int(i), int(j))
for i, j in np.argwhere(violation_matrix)
if i < j # Avoid duplicates
]
computation_time = (time.perf_counter() - start_time) * 1000
return WARPResult(
is_consistent=len(violations) == 0,
violations=violations,
computation_time_ms=computation_time,
)
# =============================================================================
# SWAPS INDEX (Apesteguia & Ballester 2015 JPE)
# =============================================================================
def compute_swaps_index(
session: ConsumerSession,
method: str = "greedy",
tolerance: float = 1e-10,
) -> SwapsIndexResult:
"""
Compute the Swaps Index (Apesteguia & Ballester 2015 JPE).
The swaps index counts the minimum number of preference relations
that must be "swapped" (reversed) to make the data GARP-consistent.
This is more interpretable than AEI: "3 swaps needed" vs "AEI = 0.92".
Args:
session: ConsumerSession with prices and quantities
method: Algorithm to use:
- "greedy": Fast heuristic (default)
- "optimal": Exact solution via ILP (slower, for small n)
tolerance: Numerical tolerance for comparisons
Returns:
SwapsIndexResult with swap count and affected pairs
Example:
>>> from prefgraph import BehaviorLog, compute_swaps_index
>>> result = compute_swaps_index(log)
>>> print(f"Need {result.swaps_count} swaps for consistency")
>>> for obs_i, obs_j in result.swap_pairs:
... print(f" Swap preference between obs {obs_i} and {obs_j}")
References:
Apesteguia, J., & Ballester, M. A. (2015). A Measure of Rationality
and Welfare. Journal of Political Economy, 123(6), 1278-1310.
"""
start_time = time.perf_counter()
# First check GARP consistency
garp_result = check_garp(session, tolerance)
if garp_result.is_consistent:
computation_time = (time.perf_counter() - start_time) * 1000
return SwapsIndexResult(
swaps_count=0,
swaps_normalized=0.0,
max_possible_swaps=session.num_observations
* (session.num_observations - 1)
// 2,
swap_pairs=[],
is_consistent=True,
method=method,
computation_time_ms=computation_time,
)
# Get revealed preference structure
R = garp_result.direct_revealed_preference
P = garp_result.strict_revealed_preference
R_star = garp_result.transitive_closure
# Find cycles that need breaking
violations = garp_result.violations
if method == "greedy":
swap_pairs = _compute_swaps_greedy(R, P, R_star, violations)
elif method == "optimal":
raise ComputationalLimitError(
"Exact ILP-based swaps computation requires an ILP solver (e.g., CPLEX, Gurobi). "
"Use method='greedy' for a heuristic approximation."
)
else:
raise ValueError(f"Unknown method: {method}. Use 'greedy' or 'optimal'.")
# Compute max possible swaps
T = session.num_observations
max_possible = T * (T - 1) // 2
# Normalize
swaps_normalized = len(swap_pairs) / max(max_possible, 1)
computation_time = (time.perf_counter() - start_time) * 1000
return SwapsIndexResult(
swaps_count=len(swap_pairs),
swaps_normalized=swaps_normalized,
max_possible_swaps=max_possible,
swap_pairs=swap_pairs,
is_consistent=False,
method=method,
computation_time_ms=computation_time,
)
def _compute_swaps_greedy(
R: NDArray[np.bool_],
P: NDArray[np.bool_],
R_star: NDArray[np.bool_],
violations: list[Cycle],
) -> list[tuple[int, int]]:
"""
Compute minimum feedback arc set using greedy heuristic.
Repeatedly removes edges that participate in the most cycles.
"""
if not violations:
return []
# Count how many cycles each edge participates in
edge_counts: dict[tuple[int, int], int] = {}
for cycle in violations:
cycle_list = list(cycle)
if len(cycle_list) < 2:
continue
for i in range(len(cycle_list) - 1):
edge = (cycle_list[i], cycle_list[i + 1])
edge_counts[edge] = edge_counts.get(edge, 0) + 1
# Handle wrap-around
if cycle_list[0] != cycle_list[-1]:
edge = (cycle_list[-1], cycle_list[0])
edge_counts[edge] = edge_counts.get(edge, 0) + 1
# Greedy: pick edges in most cycles until all cycles broken
remaining_cycles = [set(c) for c in violations]
swap_pairs = []
while remaining_cycles:
# Find edge in most remaining cycles
best_edge = None
best_count = 0
for edge, _ in sorted(edge_counts.items(), key=lambda x: -x[1]):
count = sum(1 for c in remaining_cycles if edge[0] in c and edge[1] in c)
if count > best_count:
best_count = count
best_edge = edge
if best_edge is None or best_count == 0:
break
swap_pairs.append(best_edge)
# Remove cycles containing this edge
remaining_cycles = [
c for c in remaining_cycles if not (best_edge[0] in c and best_edge[1] in c)
]
return swap_pairs
# =============================================================================
# PER-OBSERVATION CONTRIBUTIONS (Varian 1990)
# =============================================================================
def compute_observation_contributions(
session: ConsumerSession,
method: str = "cycle_count",
tolerance: float = 1e-10,
) -> ObservationContributionResult:
"""
Compute per-observation contribution to GARP violations (Varian 1990).
Identifies which observations are responsible for inconsistency.
Useful for data cleaning or identifying outliers.
Args:
session: ConsumerSession with prices and quantities
method: Method for computing contributions:
- "cycle_count": Count cycles each observation appears in (fast)
- "removal": Leave-one-out AEI improvement (accurate but slower)
tolerance: Numerical tolerance for comparisons
Returns:
ObservationContributionResult with per-observation analysis
Example:
>>> from prefgraph import BehaviorLog, compute_observation_contributions
>>> result = compute_observation_contributions(log)
>>> print(f"Base AEI: {result.base_aei:.3f}")
>>> for obs_idx, contrib in result.worst_observations[:3]:
... print(f" Obs {obs_idx}: {contrib:.2%} contribution")
References:
Varian, H. R. (1990). Goodness-of-fit in optimizing models.
Journal of Econometrics, 46(1-2), 125-140.
"""
start_time = time.perf_counter()
T = session.num_observations
# Get GARP result and compute base AEI
garp_result = check_garp(session, tolerance)
# Compute base AEI
from prefgraph.algorithms.aei import compute_aei
base_aei_result = compute_aei(session, tolerance=tolerance)
base_aei = base_aei_result.efficiency_index
# Initialize contribution tracking
contributions = np.zeros(T)
cycle_participation: dict[int, int] = {i: 0 for i in range(T)}
removal_impact: dict[int, float] = {}
if garp_result.is_consistent:
# All consistent, no contributions
computation_time = (time.perf_counter() - start_time) * 1000
return ObservationContributionResult(
contributions=contributions,
worst_observations=[],
removal_impact=removal_impact,
cycle_participation=cycle_participation,
base_aei=base_aei,
method=method,
computation_time_ms=computation_time,
)
# Count cycle participation
for cycle in garp_result.violations:
for obs in cycle:
if obs < T:
cycle_participation[obs] += 1
if method == "cycle_count":
# Contribution = fraction of cycles containing this observation
total_cycles = len(garp_result.violations)
for i in range(T):
contributions[i] = cycle_participation[i] / max(total_cycles, 1)
elif method == "removal":
# Leave-one-out analysis
for i in range(T):
if cycle_participation[i] == 0:
contributions[i] = 0.0
removal_impact[i] = 0.0
continue
# Create subset excluding observation i
mask = np.ones(T, dtype=bool)
mask[i] = False
# Only do expensive computation if observation is in cycles.
# prices/quantities are always set after __post_init__; narrow.
subset_prices = cast(NDArray[np.float64], session.prices)[mask]
subset_quantities = cast(NDArray[np.float64], session.quantities)[mask]
if len(subset_prices) > 1:
subset_session = ConsumerSession(
prices=subset_prices, quantities=subset_quantities
)
subset_aei = compute_aei(subset_session, tolerance=tolerance)
improvement = subset_aei.efficiency_index - base_aei
removal_impact[i] = improvement
contributions[i] = improvement / (1.0 - base_aei + 1e-10)
else:
removal_impact[i] = 0.0
contributions[i] = 0.0
else:
# Default to cycle_count
total_cycles = len(garp_result.violations)
for i in range(T):
contributions[i] = cycle_participation[i] / max(total_cycles, 1)
# Normalize contributions to sum to 1 (if any)
total_contrib = np.sum(contributions)
if total_contrib > 0:
contributions = contributions / total_contrib
# Get worst observations
indexed_contrib = [(i, float(contributions[i])) for i in range(T)]
worst_observations = sorted(indexed_contrib, key=lambda x: -x[1])
# Filter out zero contributions
worst_observations = [(i, c) for i, c in worst_observations if c > 0]
computation_time = (time.perf_counter() - start_time) * 1000
return ObservationContributionResult(
contributions=contributions,
worst_observations=worst_observations,
removal_impact=removal_impact,
cycle_participation=cycle_participation,
base_aei=base_aei,
method=method,
computation_time_ms=computation_time,
)
# =============================================================================
# MINIMUM COST INDEX (Dean & Martin 2016)
# =============================================================================
def compute_minimum_cost_index(
session: "ConsumerSession",
solver: str = "highs",
tolerance: float = 1e-8,
) -> "MinimumCostIndexResult":
"""
Compute the Minimum Cost Index for GARP violations.
The MCI measures the minimum monetary cost needed to break all GARP
violation cycles. It provides an alternative to the Afriat Efficiency
Index (AEI) with a more direct economic interpretation: the dollar
amount that would need to be "wasted" to make behavior consistent.
For each observation t, we find adjustment e_t >= 0 such that:
- The adjusted expenditures p^t @ x^t - e_t eliminate all GARP violations
- The total adjustment sum(e_t) is minimized
The normalized MCI is sum(e_t) / sum(p^t @ x^t), giving a proportion
of total expenditure.
Args:
session: ConsumerSession with prices and quantities
solver: LP solver to use ("highs" recommended)
tolerance: Numerical tolerance
Returns:
MinimumCostIndexResult with MCI value, adjustments, and diagnostics
Example:
>>> from prefgraph import ConsumerSession, compute_minimum_cost_index
>>> result = compute_minimum_cost_index(session)
>>> print(f"MCI: {result.mci_normalized:.4f}")
>>> print(f"Violation cost: ${result.mci_value:.2f}")
References:
Dean, M., & Martin, D. (2016). Measuring rationality with the
minimum cost of revealed preference violations.
Review of Economics and Statistics, 98(3), 524-534.
"""
start_time = time.perf_counter()
T = session.num_records
# cost_vectors/action_vectors are Optional on the dataclass but always set
# after __post_init__ (see ConsumerSession); narrow as session.py does.
P = cast(NDArray[np.float64], session.cost_vectors) # prices
Q = cast(NDArray[np.float64], session.action_vectors) # quantities
# Own expenditures
own_exp = np.sum(P * Q, axis=1) # T-vector
total_expenditure = float(np.sum(own_exp))
# Expenditure matrix E[t,s] = p^t @ x^s
E = P @ Q.T # T x T matrix
# Check if already GARP consistent
garp_result = check_garp(session, tolerance=tolerance)
if garp_result.is_consistent:
computation_time = (time.perf_counter() - start_time) * 1000
return MinimumCostIndexResult(
mci_value=0.0,
mci_normalized=0.0,
adjustments={},
cycles_broken=0,
total_expenditure=total_expenditure,
is_consistent=True,
computation_time_ms=computation_time,
)
# Build LP to find minimum cost adjustments
# Variables: e_t >= 0 for t = 0, ..., T-1
# Objective: minimize sum(e_t)
# Constraints: For each violation pair (t, s) where t R s and s P t,
# we need to break one of the relations
# Get all violation pairs from GARP result
violations = garp_result.violations
num_violations = len(violations)
if num_violations == 0:
computation_time = (time.perf_counter() - start_time) * 1000
return MinimumCostIndexResult(
mci_value=0.0,
mci_normalized=0.0,
adjustments={},
cycles_broken=0,
total_expenditure=total_expenditure,
is_consistent=True,
computation_time_ms=computation_time,
)
# Greedy heuristic: for each violation cycle, find the edge with
# smallest slack (cheapest to break) and reduce that observation's
# expenditure by slack + tolerance. Dean & Martin (2016) describe
# the exact LP, but disjunctive cycle-breaking constraints make the
# full problem NP-hard. This greedy approach provides an upper bound.
# Count cycle participation
cycle_counts = np.zeros(T)
for cycle in violations:
for obs in cycle:
cycle_counts[obs] += 1
# Build LP:
# Variables: e_t for each observation
# For each cycle (t1, t2, ..., tk), at least one e must be positive
# enough to break the preference chain
# For each cycle, break the weakest edge by reducing own expenditure
# by at least slack + tolerance so the R relation no longer holds.
# For each cycle, estimate minimum adjustment needed
adjustments = {}
cycles_broken = 0
for cycle in violations:
if len(cycle) < 2:
continue
# Find the weakest link in the cycle.
# Cycles are stored as (v0, v1, ..., vk, v0) with the first node
# repeated, so iterate len-1 edges to avoid a self-loop.
min_slack = float("inf")
weak_obs = cycle[0]
n_edges = len(cycle) - 1 if cycle[0] == cycle[-1] else len(cycle)
for i in range(n_edges):
t = cycle[i]
s = cycle[(i + 1) % len(cycle)]
# Slack is how much t R s is "strong"
slack = own_exp[t] - E[t, s]
if slack < min_slack:
min_slack = slack
weak_obs = t
# Adjust expenditure at weak observation to break cycle
if weak_obs not in adjustments:
adjustments[weak_obs] = 0.0
# Amount needed to reduce expenditure enough to break this edge
adjustment_needed = max(0.0, min_slack + tolerance)
adjustments[weak_obs] = max(adjustments[weak_obs], adjustment_needed)
cycles_broken += 1
# Compute MCI
mci_value = sum(adjustments.values())
mci_normalized = mci_value / total_expenditure if total_expenditure > 0 else 0.0
computation_time = (time.perf_counter() - start_time) * 1000
return MinimumCostIndexResult(
mci_value=mci_value,
mci_normalized=mci_normalized,
adjustments=adjustments,
cycles_broken=cycles_broken,
total_expenditure=total_expenditure,
is_consistent=False,
computation_time_ms=computation_time,
)
# =============================================================================
# TECH-FRIENDLY ALIASES
# =============================================================================
# validate_consistency: Tech-friendly name for check_garp
validate_consistency = check_garp
"""
Validate that user behavior is internally consistent.
This is the tech-friendly alias for check_garp (GARP = Generalized Axiom
of Revealed Preference). Consistent behavior indicates:
- Single user (not a shared account)
- Not a bot (bots make random inconsistent choices)
- Not confused by the UI
Returns a ConsistencyResult with:
- is_consistent: True if behavior is consistent
- violations: List of detected violation cycles
Example:
>>> from prefgraph import BehaviorLog, validate_consistency
>>> result = validate_consistency(user_log)
>>> if not result.is_consistent:
... print(f"Found {result.num_violations} violations")
"""
# validate_consistency_weak: Tech-friendly name for check_warp
validate_consistency_weak = check_warp
"""
Weak consistency check (only checks direct contradictions).
Faster than full validate_consistency but may miss transitive inconsistencies.
"""