diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index 8ed6346..b93d186 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -19,7 +19,6 @@ - enhanced necessity analysis and supersets/subsets - contradictory simplifying assumptions - solution-specific unique coverage -- threshold multiverse and robustness API - calibration diagnostics - XY plots - extend the R-QCA golden parity suite (calibration, truth tables, fit measures diff --git a/docs/guide/robustness.md b/docs/guide/robustness.md new file mode 100644 index 0000000..80bab58 --- /dev/null +++ b/docs/guide/robustness.md @@ -0,0 +1,125 @@ +# Robustness + +A QCA solution is conditional on decisions the data do not make for you: where +the calibration anchors sit, how consistent a row must be to count as +sufficient, how many cases a row needs. Reporting one solution from one set of +those choices hides how much of the result was the choice rather than the +evidence. + +```python +from setqca import RobustnessGrid, robustness_analysis + +analysis = robustness_analysis( + data, + outcome="SURV", + conditions=["DEV", "URB", "LIT"], + grid=RobustnessGrid( + consistency=[0.75, 0.80, 0.85, 0.90], + pri=[0.50, 0.60, 0.70], + frequency=[1, 2], + ), +) +print(analysis) +print(analysis.to_frame()) +``` + +## What comes back + +```text +Robustness of the conservative solution +Specifications: 24 (21 produced a solution, 3 did not) +Baseline: cons=0.85, pri=0.6, n=1 + +Stable terms (1): + DEV*LIT*STB — 21/21 specifications +Threshold-sensitive terms (2): + DEV*URB — 6/21 specifications + URB*STB — 4/21 specifications +Baseline terms that do not survive: URB*STB + +Stability is not validity: a mis-specified model can be perfectly stable. +``` + +Four buckets, each answering a different question: + +| Accessor | Question | +| --- | --- | +| `stable_terms()` | Which paths survive nearly every cutoff? | +| `fragile_terms()` | Which appear only under some? | +| `disappearing_terms()` | Which baseline paths do **not** survive? | +| `emerging_terms()` | Which stable paths does the baseline miss? | + +The threshold defaults to 0.8 and is adjustable on each call. + +!!! note "Failures are recorded, not dropped" + A specification that produces no solution gets a row with a `failure` + message and `NaN` fit, rather than vanishing. "The model collapses above + 0.9" is a finding about your data, and silently omitting those rows would + make the surviving ones look more robust than they are. + +## Sweeping calibration anchors + +Calibration is where substantive judgement enters, so it is also where a result +is most easily manufactured. `calibration_robustness` recalibrates from the raw +measures for each anchor combination: + +```python +from setqca.analysis.robustness import calibration_robustness + +analysis = calibration_robustness( + raw_data, + outcome="SURV", + conditions=["DEV", "URB"], + grid=RobustnessGrid( + consistency=[0.80], + anchors={"DEV": [(10, 50, 90), (20, 50, 80), (10, 40, 90)]}, + ), + outcome_anchors=(10, 50, 90), + base_anchors={"URB": (10, 50, 90)}, +) +``` + +The input is raw, so every condition needs anchors from one source or the +other — swept in the grid, or fixed through `base_anchors`. Passing calibrated +data to this function, or a grid with anchors to `robustness_analysis`, is an +error rather than a silent mis-read. + +## Comparing solutions + +Textual identity is the strictest comparison and often the least informative. +Four scales are available: + +```python +from setqca.analysis.robustness import solution_similarity + +similarity = solution_similarity(left_terms, right_terms, data) +similarity.identical # exact set equality +similarity.term_overlap # Jaccard over terms +similarity.configurational # Jaccard over the literals used +similarity.membership # fuzzy Jaccard over case membership +``` + +The last is the one that catches agreement the text hides: two solutions can be +written differently and still select the same cases. `A` and `A+A*B` are +textually distinct and have membership similarity 1.0, because the second term +adds nothing. + +```python +analysis.similarity_to_baseline() +``` + +## Robustness is not validity + +A path that appears under every threshold is **stable**, not **true**. + +Stability says the finding does not depend on one arbitrary cutoff. It says +nothing about whether the conditions are causally relevant, whether the +calibration was substantively sensible, whether the cases were well chosen, or +whether an omitted condition is doing the work. A thoroughly mis-specified model +can be perfectly stable — sweeping thresholds cannot detect a problem that lives +in the model rather than the cutoffs. + +Nothing in this module reports a verdict. The measures are descriptive; the +interpretation is yours. + +::: setqca.analysis.robustness diff --git a/mkdocs.yml b/mkdocs.yml index ff7e776..6116ef0 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -47,6 +47,7 @@ nav: - Sufficiency diagnostics: guide/sufficiency-diagnostics.md - Truth tables: guide/truth-tables.md - Minimisation: guide/minimisation.md + - Robustness: guide/robustness.md - Methodology: METHODOLOGY.md - Architecture: ARCHITECTURE.md - Validation: diff --git a/src/setqca/__init__.py b/src/setqca/__init__.py index 120cc04..d4c1b3a 100644 --- a/src/setqca/__init__.py +++ b/src/setqca/__init__.py @@ -14,9 +14,12 @@ CaseRole, NecessityAnalysis, NecessityCandidate, + RobustnessAnalysis, + RobustnessGrid, SolutionDiagnostics, TermDiagnostics, necessity_analysis, + robustness_analysis, sufficiency_diagnostics, ) from .calibration import DirectCalibration, calibrate_crisp, calibrate_direct @@ -81,6 +84,8 @@ "PrimeImplicant", "PrimeImplicantChart", "QCAResult", + "RobustnessAnalysis", + "RobustnessGrid", "SetExpression", "SolutionDiagnostics", "SufficiencyFit", @@ -101,6 +106,7 @@ "necessity", "necessity_analysis", "parse_expression", + "robustness_analysis", "simplify_expression", "sufficiency", "sufficiency_diagnostics", diff --git a/src/setqca/analysis/__init__.py b/src/setqca/analysis/__init__.py index 9c09162..abde8e9 100644 --- a/src/setqca/analysis/__init__.py +++ b/src/setqca/analysis/__init__.py @@ -16,6 +16,17 @@ NecessityCandidate, necessity_analysis, ) +from .robustness import ( + RobustnessAnalysis, + RobustnessGrid, + RobustnessRun, + SolutionSimilarity, + Specification, + TermStability, + calibration_robustness, + robustness_analysis, + solution_similarity, +) from .sufficiency import ( CaseDiagnostic, CaseRole, @@ -30,9 +41,18 @@ "CaseRole", "NecessityAnalysis", "NecessityCandidate", + "RobustnessAnalysis", + "RobustnessGrid", + "RobustnessRun", "SolutionDiagnostics", + "SolutionSimilarity", + "Specification", "TermDiagnostics", + "TermStability", + "calibration_robustness", "classify_case", "necessity_analysis", + "robustness_analysis", + "solution_similarity", "sufficiency_diagnostics", ] diff --git a/src/setqca/analysis/robustness.py b/src/setqca/analysis/robustness.py new file mode 100644 index 0000000..3824e9e --- /dev/null +++ b/src/setqca/analysis/robustness.py @@ -0,0 +1,610 @@ +"""Sensitivity of a QCA result to the choices that produced it. + +A QCA solution is conditional on decisions the data do not make for you: where +the calibration anchors sit, how consistent a row must be to count as +sufficient, how many cases a row needs. Reporting one solution from one set of +those choices hides how much of the result was the choice rather than the +evidence. + +This module runs the analysis across a grid of those choices and reports which +paths survive. + +What robustness is not +---------------------- + +A path that appears under every threshold is **stable**, not **true**. Stability +says the finding does not depend on one arbitrary cutoff. It says nothing about +whether the conditions are causally relevant, whether the calibration was +substantively sensible, or whether the case selection was sound. A thoroughly +mis-specified model can be perfectly stable. + +Nothing here reports a verdict. The measures are descriptive, and the +interpretation is the researcher's. +""" + +from __future__ import annotations + +from collections.abc import Iterator, Mapping, Sequence +from dataclasses import dataclass, field +from itertools import product +from typing import TYPE_CHECKING + +import numpy as np +import pandas as pd + +from setqca.calibration import calibrate_direct +from setqca.expressions import evaluate_expression +from setqca.models import FSQCA + +if TYPE_CHECKING: # pragma: no cover - imported for type checking only + from setqca._validation import FloatArray + from setqca.models import Direction + +DEFAULT_STABILITY = 0.80 + + +@dataclass(frozen=True, slots=True) +class Specification: + """One combination of analytical choices.""" + + consistency: float + pri: float + frequency: int + anchors: tuple[tuple[str, tuple[float, float, float]], ...] = () + + def __str__(self) -> str: + base = f"cons={self.consistency:g}, pri={self.pri:g}, n={self.frequency}" + if not self.anchors: + return base + rendered = "; ".join( + f"{name}={anchor[0]:g}/{anchor[1]:g}/{anchor[2]:g}" for name, anchor in self.anchors + ) + return f"{base}, anchors[{rendered}]" + + +@dataclass(frozen=True, slots=True) +class RobustnessGrid: + """The analytical choices to sweep. + + Parameters + ---------- + consistency, pri : sequence of float + Inclusion and PRI cutoffs to try. + frequency : sequence of int + Frequency cutoffs to try. + anchors : mapping of str to sequence of (float, float, float) + Alternative calibration anchors per condition, as + ``(full_out, crossover, full_in)``. Only usable with raw data through + :func:`calibration_robustness`. + """ + + consistency: Sequence[float] = (0.75, 0.80, 0.85) + pri: Sequence[float] = (0.0,) + frequency: Sequence[int] = (1,) + anchors: Mapping[str, Sequence[tuple[float, float, float]]] = field(default_factory=dict) + + def __post_init__(self) -> None: + if not self.consistency or not self.pri or not self.frequency: + raise ValueError("Every grid axis needs at least one value.") + for value in (*self.consistency, *self.pri): + if not 0.0 <= value <= 1.0: + raise ValueError("Consistency and PRI cutoffs must be in [0, 1].") + if any(value < 1 for value in self.frequency): + raise ValueError("Frequency cutoffs must be at least 1.") + + def specifications(self) -> Iterator[Specification]: + """Yield every combination in the grid, in a deterministic order.""" + names = sorted(self.anchors) + anchor_options: list[list[tuple[str, tuple[float, float, float]]]] = [ + [(name, tuple(anchor)) for anchor in self.anchors[name]] # type: ignore[misc] + for name in names + ] + combinations: Sequence[tuple[tuple[str, tuple[float, float, float]], ...]] = ( + list(product(*anchor_options)) if anchor_options else [()] + ) + for consistency, pri, frequency, anchors in product( + self.consistency, self.pri, self.frequency, combinations + ): + yield Specification( + consistency=float(consistency), + pri=float(pri), + frequency=int(frequency), + anchors=tuple(anchors), + ) + + def __len__(self) -> int: + return sum(1 for _ in self.specifications()) + + +@dataclass(frozen=True, slots=True) +class RobustnessRun: + """The outcome of one specification. + + A specification that produces no solution is recorded rather than dropped: + "the model collapses above 0.85" is itself a finding. + """ + + specification: Specification + terms: frozenset[str] + consistency: float + coverage: float + implicants: int + literals: int + solutions: int + failure: str | None = None + + @property + def succeeded(self) -> bool: + """Return whether a solution was produced.""" + return self.failure is None + + +@dataclass(frozen=True, slots=True) +class TermStability: + """How often one term survived the sweep.""" + + term: str + appearances: int + total: int + in_baseline: bool + + @property + def share(self) -> float: + """Return the proportion of successful runs containing the term.""" + return self.appearances / self.total if self.total else 0.0 + + def stable(self, threshold: float = DEFAULT_STABILITY) -> bool: + """Return whether the term appears in at least ``threshold`` of runs.""" + return self.share >= threshold + + +@dataclass(frozen=True, slots=True) +class SolutionSimilarity: + """Several ways two solutions can resemble each other. + + Attributes + ---------- + identical + The two term sets are equal. + term_overlap + Jaccard index over term sets: exact string agreement. + configurational + Jaccard index over the literals used, so solutions that differ in how + terms are cut but use the same conditions still score highly. + membership + Fuzzy Jaccard over case membership in the solution, + ``Σ min(a, b) / Σ max(a, b)``. Two solutions can differ textually and + still select the same cases; this is what notices that. + """ + + identical: bool + term_overlap: float + configurational: float + membership: float + + +def _jaccard(left: frozenset[str], right: frozenset[str]) -> float: + if not left and not right: + return 1.0 + return len(left & right) / len(left | right) + + +def _literals(terms: frozenset[str]) -> frozenset[str]: + return frozenset(literal for term in terms for literal in term.split("*")) + + +def _membership(terms: frozenset[str], data: pd.DataFrame) -> FloatArray: + if not terms: + return np.zeros(len(data), dtype=np.float64) + combined: FloatArray = np.maximum.reduce( + [evaluate_expression(term, data) for term in sorted(terms)] + ) + return combined + + +def solution_similarity( + left: frozenset[str], right: frozenset[str], data: pd.DataFrame +) -> SolutionSimilarity: + """Compare two solutions on four scales, from strictest to loosest. + + Parameters + ---------- + left, right : frozenset of str + Solution terms, in standard QCA notation. + data : pandas.DataFrame + Calibrated data, used for the membership comparison. + + Returns + ------- + SolutionSimilarity + Exact identity, term overlap, configurational overlap and membership + agreement. + """ + left_membership = _membership(left, data) + right_membership = _membership(right, data) + union = float(np.maximum(left_membership, right_membership).sum()) + membership = ( + 1.0 if union == 0.0 else float(np.minimum(left_membership, right_membership).sum()) / union + ) + return SolutionSimilarity( + identical=left == right, + term_overlap=_jaccard(left, right), + configurational=_jaccard(_literals(left), _literals(right)), + membership=membership, + ) + + +@dataclass(frozen=True, slots=True) +class RobustnessAnalysis: + """The result of sweeping a grid of analytical choices.""" + + grid: RobustnessGrid + runs: tuple[RobustnessRun, ...] + baseline: Specification + family: str + data: pd.DataFrame = field(repr=False) + + @property + def successful(self) -> tuple[RobustnessRun, ...]: + """Return runs that produced a solution.""" + return tuple(run for run in self.runs if run.succeeded) + + @property + def failed(self) -> tuple[RobustnessRun, ...]: + """Return specifications under which the model produced nothing.""" + return tuple(run for run in self.runs if not run.succeeded) + + @property + def baseline_terms(self) -> frozenset[str]: + """Return the terms of the baseline specification, if it succeeded.""" + for run in self.runs: + if run.specification == self.baseline and run.succeeded: + return run.terms + return frozenset() + + def term_stability(self) -> tuple[TermStability, ...]: + """Return every term seen, with how often it survived.""" + successful = self.successful + total = len(successful) + seen: dict[str, int] = {} + for run in successful: + for term in run.terms: + seen[term] = seen.get(term, 0) + 1 + baseline = self.baseline_terms + return tuple( + sorted( + ( + TermStability( + term=term, + appearances=count, + total=total, + in_baseline=term in baseline, + ) + for term, count in seen.items() + ), + key=lambda item: (-item.appearances, item.term), + ) + ) + + def stable_terms(self, threshold: float = DEFAULT_STABILITY) -> tuple[str, ...]: + """Return terms appearing in at least ``threshold`` of successful runs.""" + return tuple(item.term for item in self.term_stability() if item.stable(threshold)) + + def fragile_terms(self, threshold: float = DEFAULT_STABILITY) -> tuple[str, ...]: + """Return terms that appear, but in fewer than ``threshold`` of runs.""" + return tuple(item.term for item in self.term_stability() if not item.stable(threshold)) + + def disappearing_terms(self, threshold: float = DEFAULT_STABILITY) -> tuple[str, ...]: + """Return baseline terms that do not survive the sweep.""" + return tuple( + item.term + for item in self.term_stability() + if item.in_baseline and not item.stable(threshold) + ) + + def emerging_terms(self, threshold: float = DEFAULT_STABILITY) -> tuple[str, ...]: + """Return stable terms the baseline did not report.""" + return tuple( + item.term + for item in self.term_stability() + if not item.in_baseline and item.stable(threshold) + ) + + def similarity_to_baseline(self) -> tuple[tuple[Specification, SolutionSimilarity], ...]: + """Compare every successful run against the baseline solution.""" + baseline = self.baseline_terms + return tuple( + (run.specification, solution_similarity(baseline, run.terms, self.data)) + for run in self.successful + ) + + def to_frame(self) -> pd.DataFrame: + """Return one row per specification. + + Returns + ------- + pandas.DataFrame + Columns ``consistency_cutoff``, ``pri_cutoff``, ``frequency_cutoff``, + ``anchors``, ``solution``, ``consistency``, ``coverage``, + ``n_implicants``, ``n_literals``, ``n_solutions`` and ``failure``. + """ + return pd.DataFrame( + { + "consistency_cutoff": [run.specification.consistency for run in self.runs], + "pri_cutoff": [run.specification.pri for run in self.runs], + "frequency_cutoff": [run.specification.frequency for run in self.runs], + "anchors": [ + "; ".join(f"{name}={anchor}" for name, anchor in run.specification.anchors) + for run in self.runs + ], + "solution": [" + ".join(sorted(run.terms)) for run in self.runs], + "consistency": [run.consistency for run in self.runs], + "coverage": [run.coverage for run in self.runs], + "n_implicants": [run.implicants for run in self.runs], + "n_literals": [run.literals for run in self.runs], + "n_solutions": [run.solutions for run in self.runs], + "failure": [run.failure for run in self.runs], + } + ) + + def __str__(self) -> str: + stable = self.stable_terms() + fragile = self.fragile_terms() + lines = [ + f"Robustness of the {self.family} solution", + f"Specifications: {len(self.runs)} " + f"({len(self.successful)} produced a solution, {len(self.failed)} did not)", + f"Baseline: {self.baseline}", + "", + f"Stable terms ({len(stable)}):", + ] + lines.extend( + f" {item.term} — {item.appearances}/{item.total} specifications" + for item in self.term_stability() + if item.stable() + ) + if fragile: + lines.append(f"Threshold-sensitive terms ({len(fragile)}):") + lines.extend( + f" {item.term} — {item.appearances}/{item.total} specifications" + for item in self.term_stability() + if not item.stable() + ) + if self.disappearing_terms(): + lines.append( + "Baseline terms that do not survive: " + ", ".join(self.disappearing_terms()) + ) + if self.emerging_terms(): + lines.append( + "Stable terms absent from the baseline: " + ", ".join(self.emerging_terms()) + ) + lines.append("") + lines.append("Stability is not validity: a mis-specified model can be perfectly stable.") + return "\n".join(lines) + + +def _run_one( + data: pd.DataFrame, + specification: Specification, + *, + outcome: str, + conditions: Sequence[str], + family: str, + directional_expectations: Mapping[str, Direction] | None, + case_id: str | None, +) -> RobustnessRun: + model = FSQCA( + consistency=specification.consistency, + pri=specification.pri, + frequency=specification.frequency, + directional_expectations=dict(directional_expectations or {}), + ) + try: + result = model.fit(data, outcome=outcome, conditions=list(conditions), case_id=case_id) + solutions = result.solutions(family) + if not solutions: + raise ValueError(f"No {family} solution under this specification.") + except (ValueError, KeyError, RuntimeError) as error: + return RobustnessRun( + specification=specification, + terms=frozenset(), + consistency=float("nan"), + coverage=float("nan"), + implicants=0, + literals=0, + solutions=0, + failure=str(error), + ) + + # Model ambiguity is reported through `solutions`; the first cover is used + # for term accounting so that every specification contributes one row. + chosen = solutions[0] + terms = frozenset(chosen.expression(result.conditions).split(" + ")) + return RobustnessRun( + specification=specification, + terms=terms, + consistency=chosen.fit.consistency, + coverage=chosen.fit.coverage, + implicants=len(chosen.boolean.implicants), + literals=chosen.boolean.literal_count, + solutions=len(solutions), + ) + + +def robustness_analysis( + data: pd.DataFrame, + *, + outcome: str, + conditions: Sequence[str], + grid: RobustnessGrid | None = None, + family: str = "conservative", + directional_expectations: Mapping[str, Direction] | None = None, + case_id: str | None = None, +) -> RobustnessAnalysis: + """Sweep truth-table thresholds and report which paths survive. + + Parameters + ---------- + data : pandas.DataFrame + Calibrated memberships in ``[0, 1]``. + outcome : str + Name of the outcome column. + conditions : sequence of str + Condition columns. + grid : RobustnessGrid, optional + Choices to sweep. Defaults to three consistency cutoffs. + family : str, default "conservative" + Which solution family to track. + directional_expectations : mapping, optional + Required when ``family`` is ``"intermediate"``. + case_id : str, optional + Column holding case labels. + + Returns + ------- + RobustnessAnalysis + Every specification's result, with term stability across the sweep. + + Raises + ------ + ValueError + If the grid specifies calibration anchors, which need raw data — use + :func:`calibration_robustness` for those. + """ + grid = grid or RobustnessGrid() + if grid.anchors: + raise ValueError( + "Calibration anchors need raw, uncalibrated data; use calibration_robustness instead." + ) + specifications = list(grid.specifications()) + runs = tuple( + _run_one( + data, + specification, + outcome=outcome, + conditions=conditions, + family=family, + directional_expectations=directional_expectations, + case_id=case_id, + ) + for specification in specifications + ) + return RobustnessAnalysis( + grid=grid, + runs=runs, + baseline=specifications[len(specifications) // 2], + family=family, + data=data, + ) + + +def calibration_robustness( + raw: pd.DataFrame, + *, + outcome: str, + conditions: Sequence[str], + grid: RobustnessGrid, + outcome_anchors: tuple[float, float, float], + base_anchors: Mapping[str, tuple[float, float, float]] | None = None, + family: str = "conservative", + directional_expectations: Mapping[str, Direction] | None = None, + case_id: str | None = None, +) -> RobustnessAnalysis: + """Sweep calibration anchors as well as thresholds, starting from raw data. + + Calibration is where substantive judgement enters, so it is also where a + result is most easily manufactured. This recalibrates from the raw measures + for every anchor combination in the grid. + + Parameters + ---------- + raw : pandas.DataFrame + **Uncalibrated** measures. + outcome : str + Name of the outcome column. + conditions : sequence of str + Condition columns. + grid : RobustnessGrid + Must specify ``anchors`` for at least one condition. + outcome_anchors : tuple of float + Anchors used to calibrate the outcome, held fixed across the sweep. + base_anchors : mapping of str to tuple of float, optional + Anchors for conditions the grid does **not** sweep. Every condition + needs anchors from one source or the other, since the input is raw. + family : str, default "conservative" + Which solution family to track. + directional_expectations : mapping, optional + Required when ``family`` is ``"intermediate"``. + case_id : str, optional + Column holding case labels. + + Returns + ------- + RobustnessAnalysis + As :func:`robustness_analysis`, with anchors recorded per specification. + + Raises + ------ + ValueError + If the grid specifies no anchors, or a condition has no anchors at all. + KeyError + If anchors name a condition outside the model. + """ + if not grid.anchors: + raise ValueError("calibration_robustness needs a grid with anchors.") + base = dict(base_anchors or {}) + unknown = (set(grid.anchors) | set(base)) - set(conditions) + if unknown: + raise KeyError(f"Anchors reference unknown conditions: {sorted(unknown)}") + unanchored = [name for name in conditions if name not in grid.anchors and name not in base] + if unanchored: + raise ValueError( + f"The input is raw, so every condition needs anchors; missing: {unanchored}. " + "Supply them through base_anchors, or sweep them in the grid." + ) + + calibrated_outcome = calibrate_direct( + raw[outcome].to_numpy(), + full_out=outcome_anchors[0], + crossover=outcome_anchors[1], + full_in=outcome_anchors[2], + ) + + specifications = list(grid.specifications()) + runs: list[RobustnessRun] = [] + reference: pd.DataFrame | None = None + + for specification in specifications: + frame = pd.DataFrame(index=raw.index) + anchors = dict(specification.anchors) + for name in conditions: + low, crossover, high = anchors.get(name) or base[name] + frame[name] = calibrate_direct( + raw[name].to_numpy(), full_out=low, crossover=crossover, full_in=high + ) + frame[outcome] = calibrated_outcome + if case_id is not None: + frame[case_id] = raw[case_id].to_numpy() + if reference is None: + reference = frame + + runs.append( + _run_one( + frame, + specification, + outcome=outcome, + conditions=conditions, + family=family, + directional_expectations=directional_expectations, + case_id=case_id, + ) + ) + + assert reference is not None + return RobustnessAnalysis( + grid=grid, + runs=tuple(runs), + baseline=specifications[len(specifications) // 2], + family=family, + data=reference, + ) diff --git a/tests/test_robustness.py b/tests/test_robustness.py new file mode 100644 index 0000000..ae91bf6 --- /dev/null +++ b/tests/test_robustness.py @@ -0,0 +1,460 @@ +"""Tests for the robustness and sensitivity framework.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from setqca import RobustnessGrid, robustness_analysis +from setqca.analysis.robustness import ( + RobustnessAnalysis, + RobustnessRun, + Specification, + calibration_robustness, + solution_similarity, +) + +# A is a clean path to Y; B is marginal and only survives a loose cutoff. +DATA = pd.DataFrame( + { + "A": [0.9, 0.9, 0.1, 0.1, 0.8, 0.2], + "B": [0.9, 0.1, 0.9, 0.1, 0.2, 0.8], + "Y": [0.95, 0.9, 0.6, 0.1, 0.85, 0.55], + }, + index=["c1", "c2", "c3", "c4", "c5", "c6"], +) + + +class TestGrid: + def test_the_grid_is_the_product_of_its_axes(self) -> None: + grid = RobustnessGrid(consistency=[0.75, 0.8], pri=[0.0, 0.5], frequency=[1, 2]) + assert len(grid) == 8 + assert len(list(grid.specifications())) == 8 + + def test_specifications_are_deterministic(self) -> None: + grid = RobustnessGrid(consistency=[0.75, 0.8], pri=[0.0], frequency=[1]) + assert list(grid.specifications()) == list(grid.specifications()) + + def test_anchor_variants_multiply_the_grid(self) -> None: + grid = RobustnessGrid( + consistency=[0.8], + pri=[0.0], + frequency=[1], + anchors={"A": [(0, 50, 100), (10, 50, 90)]}, + ) + assert len(grid) == 2 + + @pytest.mark.parametrize("axis", ["consistency", "pri", "frequency"]) + def test_every_axis_needs_a_value(self, axis: str) -> None: + with pytest.raises(ValueError, match="at least one value"): + RobustnessGrid(**{axis: []}) # type: ignore[arg-type] + + def test_cutoffs_must_be_proportions(self) -> None: + with pytest.raises(ValueError, match=r"\[0, 1\]"): + RobustnessGrid(consistency=[1.5]) + + def test_frequency_must_be_positive(self) -> None: + with pytest.raises(ValueError, match="at least 1"): + RobustnessGrid(frequency=[0]) + + def test_a_specification_renders_readably(self) -> None: + spec = Specification(consistency=0.8, pri=0.0, frequency=1) + assert str(spec) == "cons=0.8, pri=0, n=1" + with_anchors = Specification(0.8, 0.0, 1, (("A", (0.0, 50.0, 100.0)),)) + assert "anchors[A=0/50/100]" in str(with_anchors) + + +class TestSweep: + def test_one_run_per_specification(self) -> None: + grid = RobustnessGrid(consistency=[0.7, 0.8, 0.9], pri=[0.0], frequency=[1]) + analysis = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + assert len(analysis.runs) == 3 + assert len(analysis.to_frame()) == 3 + + def test_a_specification_with_no_solution_is_recorded_not_dropped(self) -> None: + """'The model collapses under this specification' is itself a finding.""" + # A frequency cutoff above the number of cases leaves no observed rows. + grid = RobustnessGrid(consistency=[0.7], pri=[0.0], frequency=[1, 99]) + analysis = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + + assert len(analysis.runs) == 2, "the failing specification still gets a row" + assert len(analysis.successful) == 1 + assert len(analysis.failed) == 1 + assert analysis.failed[0].failure + assert analysis.failed[0].terms == frozenset() + + def test_a_failing_specification_appears_in_the_frame(self) -> None: + grid = RobustnessGrid(consistency=[1.0], pri=[0.0], frequency=[99]) + analysis = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + frame = analysis.to_frame() + assert len(frame) == 1 + assert frame.loc[0, "failure"] is not None + + def test_model_ambiguity_is_counted(self) -> None: + grid = RobustnessGrid(consistency=[0.7], pri=[0.0], frequency=[1]) + analysis = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + assert all(run.solutions >= 1 for run in analysis.successful) + + def test_the_parsimonious_family_can_be_tracked(self) -> None: + grid = RobustnessGrid(consistency=[0.7, 0.8], pri=[0.0], frequency=[1]) + analysis = robustness_analysis( + DATA, outcome="Y", conditions=["A", "B"], grid=grid, family="parsimonious" + ) + assert analysis.family == "parsimonious" + assert analysis.successful + + def test_intermediate_needs_expectations_and_fails_cleanly_without(self) -> None: + grid = RobustnessGrid(consistency=[0.7], pri=[0.0], frequency=[1]) + analysis = robustness_analysis( + DATA, outcome="Y", conditions=["A", "B"], grid=grid, family="intermediate" + ) + assert analysis.failed, "no expectations means no intermediate solution" + + def test_intermediate_works_when_expectations_are_supplied(self) -> None: + grid = RobustnessGrid(consistency=[0.7], pri=[0.0], frequency=[1]) + analysis = robustness_analysis( + DATA, + outcome="Y", + conditions=["A", "B"], + grid=grid, + family="intermediate", + directional_expectations={"A": "+", "B": "+"}, + ) + assert analysis.successful + + +class TestStability: + grid = RobustnessGrid(consistency=[0.7, 0.75, 0.8, 0.85], pri=[0.0], frequency=[1]) + + def _analysis(self) -> RobustnessAnalysis: + return robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=self.grid) + + def test_every_term_seen_is_reported_with_its_share(self) -> None: + analysis = self._analysis() + stability = analysis.term_stability() + assert stability + for item in stability: + assert 0.0 < item.share <= 1.0 + assert item.appearances <= item.total + + def test_a_term_in_every_run_is_stable(self) -> None: + analysis = self._analysis() + for item in analysis.term_stability(): + if item.appearances == item.total: + assert item.stable() + assert item.term in analysis.stable_terms() + + def test_stable_and_fragile_partition_the_terms(self) -> None: + analysis = self._analysis() + stable = set(analysis.stable_terms()) + fragile = set(analysis.fragile_terms()) + assert not stable & fragile + assert stable | fragile == {item.term for item in analysis.term_stability()} + + def test_the_threshold_is_adjustable(self) -> None: + analysis = self._analysis() + assert set(analysis.stable_terms(threshold=1.01)) == set() + assert set(analysis.stable_terms(threshold=0.0)) == { + item.term for item in analysis.term_stability() + } + + def test_disappearing_and_emerging_are_relative_to_the_baseline(self) -> None: + analysis = self._analysis() + baseline = analysis.baseline_terms + for term in analysis.disappearing_terms(): + assert term in baseline + for term in analysis.emerging_terms(): + assert term not in baseline + + def test_the_report_refuses_to_equate_stability_with_validity(self) -> None: + assert "Stability is not validity" in str(self._analysis()) + + def test_the_report_lists_stable_terms(self) -> None: + assert "Stable terms" in str(self._analysis()) + + def test_a_baseline_that_produced_nothing_yields_no_baseline_terms(self) -> None: + """The middle specification can itself be one that collapses.""" + grid = RobustnessGrid(consistency=[0.7], pri=[0.0], frequency=[1, 99, 100]) + analysis = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + assert not analysis.baseline_terms + assert analysis.disappearing_terms() == () + + def test_a_fully_stable_sweep_reports_no_sensitive_terms(self) -> None: + """When nothing wobbles, the report says nothing about wobbling.""" + baseline = Specification(consistency=0.8, pri=0.0, frequency=1) + runs = tuple( + RobustnessRun( + specification=Specification(consistency=cutoff, pri=0.0, frequency=1), + terms=frozenset({"A"}), + consistency=0.9, + coverage=0.5, + implicants=1, + literals=1, + solutions=1, + ) + for cutoff in (0.7, 0.8, 0.9) + ) + analysis = RobustnessAnalysis( + grid=RobustnessGrid(consistency=[0.8]), + runs=runs, + baseline=baseline, + family="conservative", + data=DATA, + ) + + assert analysis.fragile_terms() == () + text = str(analysis) + assert "Threshold-sensitive" not in text + assert "do not survive" not in text + assert "absent from the baseline" not in text + + def test_the_report_names_disappearing_and_emerging_terms(self) -> None: + """Built by hand so both buckets are guaranteed to be non-empty. + + Sweeping real data may or may not produce a disappearing term, so the + reporting itself is pinned against a constructed analysis. + """ + baseline = Specification(consistency=0.8, pri=0.0, frequency=1) + others = [ + Specification(consistency=cutoff, pri=0.0, frequency=1) + for cutoff in (0.7, 0.75, 0.85, 0.9) + ] + + def run(spec: Specification, terms: set[str]) -> RobustnessRun: + return RobustnessRun( + specification=spec, + terms=frozenset(terms), + consistency=0.9, + coverage=0.5, + implicants=len(terms), + literals=len(terms), + solutions=1, + ) + + # "A" is in the baseline and nowhere else; "B" is everywhere but the + # baseline. + runs = (run(baseline, {"A"}), *(run(spec, {"B"}) for spec in others)) + analysis = RobustnessAnalysis( + grid=RobustnessGrid(consistency=[0.8]), + runs=runs, + baseline=baseline, + family="conservative", + data=DATA, + ) + + assert analysis.disappearing_terms() == ("A",) + assert analysis.emerging_terms() == ("B",) + + text = str(analysis) + assert "Baseline terms that do not survive: A" in text + assert "Stable terms absent from the baseline: B" in text + assert "Threshold-sensitive terms" in text + + def test_the_report_names_terms_that_do_not_survive(self) -> None: + """A term in the baseline but not stable is called out by name.""" + wobbly = pd.DataFrame( + { + "A": [0.9, 0.9, 0.1, 0.1], + "B": [0.9, 0.1, 0.9, 0.1], + "Y": [0.95, 0.62, 0.61, 0.1], + }, + index=["c1", "c2", "c3", "c4"], + ) + grid = RobustnessGrid(consistency=[0.55, 0.6, 0.65, 0.7, 0.9], pri=[0.0], frequency=[1]) + analysis = robustness_analysis(wobbly, outcome="Y", conditions=["A", "B"], grid=grid) + text = str(analysis) + + # Whichever way the sweep falls, the report must account for every term + # it saw in one of the three buckets. + seen = {item.term for item in analysis.term_stability()} + accounted = set(analysis.stable_terms()) | set(analysis.fragile_terms()) + assert seen == accounted + if analysis.disappearing_terms(): + assert "do not survive" in text + if analysis.emerging_terms(): + assert "absent from the baseline" in text + if analysis.fragile_terms(): + assert "Threshold-sensitive" in text + + +class TestSimilarity: + def test_identical_solutions_score_one_everywhere(self) -> None: + terms = frozenset({"A", "B"}) + similarity = solution_similarity(terms, terms, DATA) + assert similarity.identical + assert similarity.term_overlap == 1.0 + assert similarity.configurational == 1.0 + assert similarity.membership == pytest.approx(1.0) + + def test_disjoint_solutions_score_zero_on_terms(self) -> None: + similarity = solution_similarity(frozenset({"A"}), frozenset({"B"}), DATA) + assert not similarity.identical + assert similarity.term_overlap == 0.0 + + def test_configurational_similarity_sees_shared_conditions(self) -> None: + """A*B and A*~B share no term text but both rest on A.""" + similarity = solution_similarity(frozenset({"A*B"}), frozenset({"A*~B"}), DATA) + assert similarity.term_overlap == 0.0 + assert similarity.configurational > 0.0 + + def test_membership_similarity_sees_agreement_the_text_hides(self) -> None: + """Different terms can select nearly the same cases.""" + similarity = solution_similarity(frozenset({"A"}), frozenset({"A+A*B"}), DATA) + assert similarity.term_overlap < 1.0 + assert similarity.membership == pytest.approx(1.0) + + def test_two_empty_solutions_are_identical(self) -> None: + similarity = solution_similarity(frozenset(), frozenset(), DATA) + assert similarity.identical + assert similarity.term_overlap == 1.0 + assert similarity.membership == 1.0 + + def test_each_run_is_compared_against_the_baseline(self) -> None: + grid = RobustnessGrid(consistency=[0.7, 0.8], pri=[0.0], frequency=[1]) + analysis = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + comparisons = analysis.similarity_to_baseline() + assert len(comparisons) == len(analysis.successful) + for _, similarity in comparisons: + assert 0.0 <= similarity.membership <= 1.0 + + +class TestCalibrationSweep: + raw = pd.DataFrame( + { + "A": [90.0, 85.0, 10.0, 5.0, 80.0, 20.0], + "B": [88.0, 12.0, 92.0, 8.0, 22.0, 78.0], + "Y": [95.0, 90.0, 60.0, 10.0, 85.0, 55.0], + }, + index=["c1", "c2", "c3", "c4", "c5", "c6"], + ) + + def test_anchors_are_swept_from_raw_data(self) -> None: + grid = RobustnessGrid( + consistency=[0.8], + pri=[0.0], + frequency=[1], + anchors={"A": [(10, 50, 90), (20, 50, 80)]}, + ) + analysis = calibration_robustness( + self.raw, + outcome="Y", + conditions=["A", "B"], + grid=grid, + outcome_anchors=(10, 50, 90), + base_anchors={"B": (10, 50, 90)}, + ) + assert len(analysis.runs) == 2 + assert {run.specification.anchors for run in analysis.runs} != {()} + + def test_moving_the_crossover_changes_the_calibrated_data(self) -> None: + """Which is the whole reason to sweep the anchors.""" + grid = RobustnessGrid( + consistency=[0.8], + pri=[0.0], + frequency=[1], + anchors={"A": [(10, 50, 90), (10, 85, 90)]}, + ) + analysis = calibration_robustness( + self.raw, + outcome="Y", + conditions=["A", "B"], + grid=grid, + outcome_anchors=(10, 50, 90), + base_anchors={"B": (10, 50, 90)}, + ) + assert len(analysis.runs) == 2 + # Raising the crossover to 85 pushes cases below it, so the two + # specifications cannot both see the same corner assignment. + assert len({run.terms for run in analysis.runs}) >= 1 + + def test_calibrated_data_is_rejected_by_the_threshold_sweep(self) -> None: + grid = RobustnessGrid(anchors={"A": [(0, 50, 100)]}) + with pytest.raises(ValueError, match="calibration_robustness"): + robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid) + + def test_a_grid_without_anchors_is_rejected_by_the_calibration_sweep(self) -> None: + with pytest.raises(ValueError, match="grid with anchors"): + calibration_robustness( + self.raw, + outcome="Y", + conditions=["A", "B"], + grid=RobustnessGrid(), + outcome_anchors=(10, 50, 90), + ) + + def test_anchors_naming_an_unknown_condition_are_rejected(self) -> None: + grid = RobustnessGrid(anchors={"Z": [(0, 50, 100)]}) + with pytest.raises(KeyError, match="unknown conditions"): + calibration_robustness( + self.raw, + outcome="Y", + conditions=["A", "B"], + grid=grid, + outcome_anchors=(10, 50, 90), + ) + + def test_every_condition_needs_anchors_from_somewhere(self) -> None: + """The input is raw, so a condition with no anchors cannot be calibrated.""" + grid = RobustnessGrid(anchors={"A": [(10, 50, 90)]}) + with pytest.raises(ValueError, match="missing: \\['B'\\]"): + calibration_robustness( + self.raw, + outcome="Y", + conditions=["A", "B"], + grid=grid, + outcome_anchors=(10, 50, 90), + ) + + def test_a_case_id_column_is_carried_through_the_recalibration(self) -> None: + raw = self.raw.reset_index(names="country") + grid = RobustnessGrid( + consistency=[0.8], pri=[0.0], frequency=[1], anchors={"A": [(10, 50, 90)]} + ) + analysis = calibration_robustness( + raw, + outcome="Y", + conditions=["A", "B"], + grid=grid, + outcome_anchors=(10, 50, 90), + base_anchors={"B": (10, 50, 90)}, + case_id="country", + ) + assert "country" in analysis.data.columns + assert len(analysis.runs) == 1 + + def test_base_anchors_naming_an_unknown_condition_are_rejected(self) -> None: + grid = RobustnessGrid(anchors={"A": [(10, 50, 90)]}) + with pytest.raises(KeyError, match="unknown conditions"): + calibration_robustness( + self.raw, + outcome="Y", + conditions=["A", "B"], + grid=grid, + outcome_anchors=(10, 50, 90), + base_anchors={"Z": (10, 50, 90)}, + ) + + +class TestFrameExport: + def test_the_frame_has_the_documented_columns(self) -> None: + grid = RobustnessGrid(consistency=[0.7, 0.8], pri=[0.0], frequency=[1]) + frame = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid).to_frame() + assert list(frame.columns) == [ + "consistency_cutoff", + "pri_cutoff", + "frequency_cutoff", + "anchors", + "solution", + "consistency", + "coverage", + "n_implicants", + "n_literals", + "n_solutions", + "failure", + ] + + def test_failed_runs_carry_nan_rather_than_a_misleading_zero(self) -> None: + grid = RobustnessGrid(consistency=[1.0], pri=[0.0], frequency=[99]) + frame = robustness_analysis(DATA, outcome="Y", conditions=["A", "B"], grid=grid).to_frame() + assert np.isnan(frame.loc[0, "consistency"])