"""Abstract Choice Theory algorithms for menu-based preference analysis.
This module implements revealed preference axioms for abstract choice data
(choices from menus without prices). Based on Chapters 1-2 of
Chambers & Echenique (2016) "Revealed Preference Theory".
Tech-Friendly Names (Primary):
- validate_menu_warp(): Check WARP for menu choices
- validate_menu_sarp(): Check SARP for menu choices
- validate_menu_consistency(): Check full rationalizability (Congruence)
- compute_menu_efficiency(): Houtman-Maks efficiency index
- fit_menu_preferences(): Recover ordinal preference ranking
Economics Names (Legacy Aliases):
- check_abstract_warp() -> validate_menu_warp()
- check_abstract_sarp() -> validate_menu_sarp()
- check_congruence() -> validate_menu_consistency()
- compute_abstract_efficiency() -> compute_menu_efficiency()
- recover_ordinal_utility() -> fit_menu_preferences()
"""
from __future__ import annotations
import time
from typing import TYPE_CHECKING
import numpy as np
from numpy.typing import NDArray
from prefgraph.core.result import (
AbstractWARPResult,
AbstractSARPResult,
CongruenceResult,
HoutmanMaksAbstractResult,
OrdinalUtilityResult,
)
from prefgraph.core.types import Cycle
from prefgraph.graph.transitive_closure import floyd_warshall_transitive_closure
from prefgraph._kernels import (
bfs_find_path_numba,
find_symmetric_pairs_bool_numba,
topological_sort_numba,
)
if TYPE_CHECKING:
from prefgraph.core.session import MenuChoiceLog
def _find_cycle_from_pair(
R: NDArray[np.bool_],
start: int,
end: int,
) -> Cycle | None:
"""
Find a cycle between two nodes using BFS (numba-accelerated).
Args:
R: Revealed preference adjacency matrix
start: Starting node
end: Ending node (should be reachable from start and vice versa)
Returns:
Tuple of node indices forming the cycle, or None if not found
"""
# BFS from start to end using numba kernel
path_to_end = bfs_find_path_numba(R, np.int64(start), np.int64(end))
if path_to_end[0] == -1:
return None
# BFS from end back to start using numba kernel
path_back = bfs_find_path_numba(R, np.int64(end), np.int64(start))
if path_back[0] == -1:
return None
# Combine paths to form cycle
# path_to_end: start -> ... -> end -> start (numba returns with start at end)
# path_back: end -> ... -> start -> end (numba returns with end at end)
# We need: start -> ... -> end -> ... -> start
# path_to_end already goes start->...->end->start, just need start->end path
# path_back goes end->...->start->end, we need end->start path
# Extract path from start to end (exclude the cycling back)
path_to_end_list = list(path_to_end[:-1]) # Remove the repeat at end
path_back_list = list(path_back[:-1]) # Remove the repeat at end
# Cycle: start -> ... -> end -> ... -> start
cycle = path_to_end_list[:-1] + path_back_list
return tuple(cycle)
MENU_HM_ILP_THRESHOLD = 60
def _menu_houtman_maks_exact(
obs_edges: list[list[tuple[int, int]]],
n_items: int,
n_obs: int,
) -> list[int]:
"""Exact menu Houtman-Maks over observations.
Maximise the number of kept observations subject to the kept item graph
being acyclic, encoded as a ranking ILP: a binary keep-flag z_o per
observation and a real rank r_i per item, with r[chosen] >= r[other] + 1 for
every edge of a kept observation (the constraint is inactive when z_o = 0).
A consistent ranking exists iff the kept graph is acyclic, i.e. SARP holds.
"""
from scipy.optimize import Bounds, LinearConstraint, milp
big_m = float(n_items + 1)
n_vars = n_obs + n_items
rows: list[list[float]] = []
upper: list[float] = []
for o, edges in enumerate(obs_edges):
for chosen, other in edges:
row = [0.0] * n_vars
row[n_obs + chosen] -= 1.0
row[n_obs + other] += 1.0
row[o] += big_m
rows.append(row)
upper.append(big_m - 1.0)
if not rows:
return []
c = np.zeros(n_vars)
c[:n_obs] = -1.0 # maximize sum(z)
integrality = np.zeros(n_vars, dtype=int)
integrality[:n_obs] = 1
lb = np.zeros(n_vars)
ub = np.concatenate([np.ones(n_obs), np.full(n_items, float(n_items))])
res = milp(
c,
constraints=LinearConstraint(np.array(rows), -np.inf, np.array(upper)),
integrality=integrality,
bounds=Bounds(lb, ub),
)
if not getattr(res, "success", False) or res.x is None:
return _menu_houtman_maks_greedy(obs_edges, n_items, n_obs)
z = res.x[:n_obs]
return [o for o in range(n_obs) if z[o] < 0.5]
def _menu_houtman_maks_greedy(
obs_edges: list[list[tuple[int, int]]],
n_items: int,
n_obs: int,
) -> list[int]:
"""Greedy upper bound on removals: drop the observation touching the most
SARP-violating item pairs until the remaining item graph is acyclic. Used
only above the exact-ILP threshold; it over-removes relative to the optimum."""
kept = set(range(n_obs))
removed: list[int] = []
while True:
graph = np.zeros((n_items, n_items), dtype=np.bool_)
for o in kept:
for c, k in obs_edges[o]:
graph[c, k] = True
closure = floyd_warshall_transitive_closure(graph)
viol = closure & closure.T
np.fill_diagonal(viol, False)
if not np.any(viol):
break
best_o, best_score = -1, -1
for o in kept:
score = sum(1 for c, k in obs_edges[o] if viol[c, k] or viol[k, c])
if score > best_score:
best_o, best_score = o, score
if best_o < 0:
break
kept.discard(best_o)
removed.append(best_o)
return removed
# =============================================================================
# LEGACY ALIASES (Economics terminology)
# =============================================================================
check_abstract_warp = validate_menu_warp
"""Legacy alias: use validate_menu_warp instead."""
check_abstract_sarp = validate_menu_sarp
"""Legacy alias: use validate_menu_sarp instead."""
check_sarp = validate_menu_sarp
"""Compatibility alias: use validate_menu_sarp instead."""
menu_sarp_check = validate_menu_sarp
"""Technical alias: use validate_menu_sarp instead."""
menu_warp_check = validate_menu_warp
"""Technical alias: use validate_menu_warp instead."""
check_congruence = validate_menu_consistency
"""Legacy alias: use validate_menu_consistency instead."""
compute_abstract_efficiency = compute_menu_efficiency
"""Legacy alias: use compute_menu_efficiency instead."""
recover_ordinal_utility = fit_menu_preferences
"""Legacy alias: use fit_menu_preferences instead."""