diff --git a/README.md b/README.md index c16ef50..84bd34a 100644 --- a/README.md +++ b/README.md @@ -197,7 +197,6 @@ contract. ## Non-goals for 0.1 - claiming complete parity with R `QCA`; -- mvQCA; - tQCA; - CCubes/eQMC performance parity. diff --git a/docs/ROADMAP.md b/docs/ROADMAP.md index b93d186..d4be7a7 100644 --- a/docs/ROADMAP.md +++ b/docs/ROADMAP.md @@ -27,10 +27,8 @@ - optional R-compatible calibration snapping, so extreme memberships can be reported exactly as R does when replicating an existing analysis -## 0.3 — multi-value and performance +## 0.3 — performance -- mvQCA -- categorical-set expressions - faster bitset/cube minimiser - prime-implicant consistency filters - row dominance diff --git a/docs/guide/multivalue.md b/docs/guide/multivalue.md new file mode 100644 index 0000000..dd2e541 --- /dev/null +++ b/docs/guide/multivalue.md @@ -0,0 +1,101 @@ +# Multi-value QCA + +A multi-value condition takes one of several unordered categories — regime type, +welfare regime, sector — rather than being present or absent. Forcing such a +condition into a binary set either loses information or invents a dichotomy the +concept does not have. + +```python +from setqca.multivalue import MVQCA + +result = MVQCA(consistency=0.8).fit( + data, outcome="Y", conditions=["regime", "wealth"] +) +print(result) +print(result.truth_table.to_frame()) +print(result.summary_frame("parsimonious")) +``` + +Conditions hold integer category codes from `0`; the outcome is a membership in +`[0, 1]`. The workflow deliberately mirrors `FSQCA` and `CSQCA` — moving between +them is a change of estimator, not a change of method. + +## Notation + +`regime{0,2}*wealth{1}` reads "regime is 0 or 2, and wealth is 1". A condition +allowing *every* level constrains nothing and is omitted from the expression, so +the binary case reduces to familiar QCA notation. + +## The property space + +```python +result.domain # regime{0,1,2}, wealth{0,1} +result.domain.size # 6 logically possible configurations +``` + +Configurations are indexed in mixed radix, which generalises the binary minterm +and reduces to it exactly when every condition has two levels. + +!!! warning "Declare levels that have no cases" + Levels are inferred from the data, which understates a category that is + theoretically possible but happens to be unobserved. That matters: an + unobserved level is a *remainder*, and remainders change the parsimonious + solution. + + ```python + MVQCA(levels={"regime": 4, "wealth": 2}).fit(...) + ``` + + Declaring fewer levels than the data contain is an error. + +## Why not Boolean dummies + +The obvious shortcut is to encode `A{0,1,2}` as three binary indicators and +reuse the binary minimiser. **That transformation does not preserve the +semantics.** + +The binary space contains points such as `A_0 = A_1 = 1` — a case that is +simultaneously in two mutually exclusive categories, which corresponds to no +configuration at all. The minimiser is free to build implicants across those +points, producing terms that look valid and describe nothing. Recovering a +multi-value expression afterwards requires exactly the mutual-exclusivity +constraints the encoding threw away. + +So the cube algebra is implemented directly. A cube allows a **set** of levels +per condition, and merging generalises the binary rule: + +> two cubes that agree on every condition but one merge into a single cube whose +> set at that condition is the union of the two. + +Because the two cubes agree everywhere else, the merged cube covers exactly +their union and nothing more — the same property the binary rule relies on. A +test asserts precisely that, and another asserts every cube covers only real +configurations. + +The exact cover is then solved by the **same verified solver the binary engine +uses**, so both inherit one exactness guarantee rather than two implementations. +Minimisation is checked against exhaustive enumeration for four different level +combinations, and against the binary minimiser for every three-condition +problem. + +## Agreement with R + +R `QCA` supports multi-value and writes literals as `regime[2]`. The truth table +and the parsimonious solution match exactly on the benchmarks in +`validation/fixtures/r_qca.json`. + +The conservative solution can differ in *representation*: + +```text +R: regime[2] + regime[1]*wealth[1] +setqca: regime{2} + regime{1,2}*wealth{1} +``` + +Both cover the same configurations and both cost two terms and three literals, +so both are minimal. The difference is that R writes single-value literals only, +while `setqca` also forms subset literals — and here R's `regime[1]*wealth[1]` +is a **proper subset** of `regime{1,2}*wealth{1}`, so R's term is not a prime +implicant. The parity tests therefore compare cost and coverage rather than +text, which is the comparison that carries meaning. + +::: setqca.multivalue diff --git a/docs/guide/truth-tables.md b/docs/guide/truth-tables.md index 84bb326..096eabc 100644 --- a/docs/guide/truth-tables.md +++ b/docs/guide/truth-tables.md @@ -87,11 +87,11 @@ recorded. A row's outcome code alone conflates situations that call for different responses: ```python -table.positive_rows() # coded "1" -table.negative_rows() # coded "0" +table.positive_rows() # coded "1" +table.negative_rows() # coded "0" table.contradictions() # coded "C" -table.remainders() # coded "R" -table.excluded_rows() # kept out by a *threshold*, not by the evidence +table.remainders() # coded "R" +table.excluded_rows() # kept out by a *threshold*, not by the evidence print(table.summary()) ``` @@ -124,8 +124,8 @@ stored and re-minimised without recalibrating or rebuilding: text = table.to_json() restored = TruthTable.from_json(text) -restored.minimize() # conservative -restored.minimize(include_remainders=True) # parsimonious +restored.minimize() # conservative +restored.minimize(include_remainders=True) # parsimonious ``` Both agree with the estimator exactly — there are tests asserting so. Only the diff --git a/mkdocs.yml b/mkdocs.yml index 6116ef0..348a306 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 + - Multi-value QCA: guide/multivalue.md - Robustness: guide/robustness.md - Methodology: METHODOLOGY.md - Architecture: ARCHITECTURE.md diff --git a/src/setqca/__init__.py b/src/setqca/__init__.py index 62ea289..f615fbc 100644 --- a/src/setqca/__init__.py +++ b/src/setqca/__init__.py @@ -65,6 +65,7 @@ minimize_chart, ) from .models import CSQCA, FSQCA, Direction +from .multivalue import MVQCA, MultiValueDomain, MultiValueResult, MultiValueTruthTable from .results import FittedSolution, QCAResult from .sets import Condition, Intersection, Negation, SetExpression, Union from .truth_table import TruthCode, TruthTable, TruthTableRow, build_truth_table @@ -77,6 +78,7 @@ __all__ = [ "CSQCA", "FSQCA", + "MVQCA", "AnchorSuggestion", "BooleanSolution", "CalibrationDiagnostics", @@ -98,6 +100,9 @@ "Intersection", "MinimalCover", "MinimizationResult", + "MultiValueDomain", + "MultiValueResult", + "MultiValueTruthTable", "NecessityAnalysis", "NecessityCandidate", "NecessityFit", diff --git a/src/setqca/minimize/qmc.py b/src/setqca/minimize/qmc.py index 4c98927..f319b76 100644 --- a/src/setqca/minimize/qmc.py +++ b/src/setqca/minimize/qmc.py @@ -3,6 +3,7 @@ from __future__ import annotations from collections import defaultdict +from collections.abc import Sequence from dataclasses import dataclass from .implicant import Implicant, minterm_to_implicant @@ -119,32 +120,67 @@ def exact_minimum_covers( """ if not on_set: return (BooleanSolution(()),) - cover_map = {m: tuple(i for i, p in enumerate(primes) if p.covers(m)) for m in on_set} - if any(not choices for choices in cover_map.values()): - raise RuntimeError("Prime-implicant chart cannot cover every positive row.") + covered = [frozenset(m for m in on_set if prime.covers(m)) for prime in primes] + literals = [prime.literals for prime in primes] + choices = solve_minimum_cover(covered, literals, on_set, max_solutions=max_solutions) + return tuple(BooleanSolution(tuple(primes[i] for i in indices)) for indices in choices) + + +def solve_minimum_cover( + covered: Sequence[frozenset[int]], + literals: Sequence[int], + on_set: set[int], + *, + max_solutions: int = 256, +) -> tuple[tuple[int, ...], ...]: + """Solve a covering problem exactly, independently of what is being covered. + + The exactness guarantee lives here, so every minimiser in the package — + binary and multi-value alike — shares one verified implementation rather + than repeating the search. + + Parameters + ---------- + covered : sequence of frozenset of int + For each candidate, the elements of ``on_set`` it covers. + literals : sequence of int + Cost of each candidate, used as the secondary objective. + on_set : set of int + Elements that must be covered. + max_solutions : int, default 256 + Upper bound on the number of tied minimum covers returned. + + Returns + ------- + tuple of tuple of int + Candidate indices, one tuple per tied minimum cover. - covered_by = { - index: frozenset(m for m in on_set if primes[index].covers(m)) - for index in {i for choices in cover_map.values() for i in choices} - } + Raises + ------ + RuntimeError + If some element is covered by no candidate. + """ + cover_map = {m: tuple(i for i, reach in enumerate(covered) if m in reach) for m in on_set} + if any(not options for options in cover_map.values()): + raise RuntimeError("Prime-implicant chart cannot cover every positive row.") - # A minterm with a single candidate forces that prime into every cover. - essential = frozenset(choices[0] for choices in cover_map.values() if len(choices) == 1) + # A minterm with a single candidate forces that candidate into every cover. + essential = frozenset(options[0] for options in cover_map.values() if len(options) == 1) start_uncovered = ( - frozenset(on_set).difference(*(covered_by[i] for i in essential)) - if (essential) + frozenset(on_set).difference(*(covered[i] for i in essential)) + if essential else frozenset(on_set) ) def lower_bound(uncovered: frozenset[int]) -> int: - """Return a lower bound on the number of further primes required.""" + """Return a lower bound on the number of further candidates required.""" blocked: set[int] = set() bound = 0 for minterm in sorted(uncovered, key=lambda m: len(cover_map[m])): - choices = cover_map[minterm] - if blocked.isdisjoint(choices): + options = cover_map[minterm] + if blocked.isdisjoint(options): bound += 1 - blocked.update(choices) + blocked.update(options) return bound best_cost: tuple[int, int] | None = None @@ -153,7 +189,7 @@ def lower_bound(uncovered: frozenset[int]) -> int: def search(chosen: frozenset[int], uncovered: frozenset[int]) -> None: nonlocal best_cost - current_cost = (len(chosen), sum(primes[i].literals for i in chosen)) + current_cost = (len(chosen), sum(literals[i] for i in chosen)) if not uncovered: if best_cost is None or current_cost < best_cost: best_cost = current_cost @@ -164,8 +200,8 @@ def search(chosen: frozenset[int], uncovered: frozenset[int]) -> None: if best_cost is not None: if current_cost >= best_cost: return - # Any completion needs at least `lower_bound` further primes, and - # primes never reduce the literal count. + # Any completion needs at least `lower_bound` further candidates, + # and candidates never reduce the cost. if (current_cost[0] + lower_bound(uncovered), current_cost[1]) > best_cost: return previous = seen.get(uncovered) @@ -174,18 +210,18 @@ def search(chosen: frozenset[int], uncovered: frozenset[int]) -> None: if previous is None or current_cost < previous: seen[uncovered] = current_cost - # Branch on the minterm with the fewest candidates: every cover must + # Branch on the element with the fewest candidates: every cover must # contain one of them, so this is a complete and narrow branching rule. target = min(uncovered, key=lambda m: len(cover_map[m])) - candidates = sorted( + ordered = sorted( cover_map[target], - key=lambda i: (-len(covered_by[i] & uncovered), primes[i].literals, i), + key=lambda i: (-len(covered[i] & uncovered), literals[i], i), ) - for index in candidates: - search(chosen | {index}, uncovered - covered_by[index]) + for index in ordered: + search(chosen | {index}, uncovered - covered[index]) search(essential, start_uncovered) - return tuple(BooleanSolution(tuple(primes[i] for i in indices)) for indices in sorted(best)) + return tuple(sorted(best)) def minimize( diff --git a/src/setqca/multivalue/__init__.py b/src/setqca/multivalue/__init__.py new file mode 100644 index 0000000..458031c --- /dev/null +++ b/src/setqca/multivalue/__init__.py @@ -0,0 +1,52 @@ +"""Multi-value QCA: categorical conditions with more than two levels. + +A multi-value condition takes one of several unordered categories — regime type, +welfare regime, sector — rather than being present or absent. Forcing such a +condition into a binary set either loses information or invents a dichotomy the +concept does not have. + +The cube algebra is implemented directly rather than by encoding categories as +Boolean indicators; see :mod:`setqca.multivalue._cube` for why that encoding is +unsound. The exact cover is solved by the same verified solver the binary +engine uses, so both inherit one exactness guarantee. + +Examples +-------- +>>> import pandas as pd +>>> from setqca.multivalue import MVQCA +>>> data = pd.DataFrame( +... {"regime": [0, 1, 2, 1], "wealth": [0, 1, 1, 0], "Y": [0.1, 0.9, 0.9, 0.2]} +... ) +>>> result = MVQCA(consistency=0.8).fit(data, outcome="Y", conditions=["regime", "wealth"]) +>>> print(result.summary_frame()) # doctest: +SKIP +""" + +from __future__ import annotations + +from ._cube import ( + MultiValueCube, + MultiValueSolution, + minimize_multivalue, + prime_cubes, +) +from ._domain import MultiValueDomain +from ._model import ( + MVQCA, + MultiValueResult, + MultiValueRow, + MultiValueTruthTable, + build_multivalue_truth_table, +) + +__all__ = [ + "MVQCA", + "MultiValueCube", + "MultiValueDomain", + "MultiValueResult", + "MultiValueRow", + "MultiValueSolution", + "MultiValueTruthTable", + "build_multivalue_truth_table", + "minimize_multivalue", + "prime_cubes", +] diff --git a/src/setqca/multivalue/_cube.py b/src/setqca/multivalue/_cube.py new file mode 100644 index 0000000..772a860 --- /dev/null +++ b/src/setqca/multivalue/_cube.py @@ -0,0 +1,226 @@ +"""Multi-value cubes and their exact minimisation. + +A cube in a multi-value property space allows a **set** of levels for each +condition, rather than a single value or a don't-care. ``A{0,2}*B{1}`` is a +legitimate term, and a condition whose set contains every level is simply +absent from the expression. + +Why not Boolean dummies +----------------------- + +The obvious shortcut is to encode ``A{0,1,2}`` as three binary indicators and +reuse the binary minimiser. That transformation does **not** preserve the +semantics. The binary space contains points such as ``A_0 = A_1 = 1``, which +correspond to no configuration at all, and the minimiser is free to build +implicants across them — producing terms that look valid and describe nothing. +Recovering a multi-value expression afterwards requires exactly the +mutual-exclusivity constraints the encoding discarded. + +So the cube algebra is implemented directly. Merging is the generalisation of +the binary rule: + + two cubes that agree on every condition but one merge into a single cube + whose set at that condition is the union of the two. + +Because the two cubes agree everywhere else, the merged cube covers exactly +their union and nothing more, so merging can never introduce coverage of a +configuration that was not already covered. That is the property the binary +rule relies on, and it holds unchanged for sets. + +The exact cover is then solved by the same verified solver the binary engine +uses. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import TYPE_CHECKING + +from setqca.minimize.qmc import solve_minimum_cover + +if TYPE_CHECKING: # pragma: no cover - imported for type checking only + from ._domain import MultiValueDomain + + +@dataclass(frozen=True, slots=True) +class MultiValueCube: + """A conjunction allowing a set of levels for each condition.""" + + pattern: tuple[frozenset[int], ...] + + @classmethod + def from_configuration(cls, values: tuple[int, ...]) -> MultiValueCube: + """Build the cube covering exactly one configuration.""" + return cls(tuple(frozenset({value}) for value in values)) + + def literals(self, domain: MultiValueDomain) -> int: + """Return the number of conditions the cube actually constrains.""" + return sum( + 1 + for allowed, count in zip(self.pattern, domain.levels, strict=True) + if len(allowed) < count + ) + + def is_tautology(self, domain: MultiValueDomain) -> bool: + """Return whether the cube constrains nothing.""" + return self.literals(domain) == 0 + + def covers_values(self, values: tuple[int, ...]) -> bool: + """Return whether a configuration falls inside the cube.""" + return all(value in allowed for value, allowed in zip(values, self.pattern, strict=True)) + + def covers(self, index: int, domain: MultiValueDomain) -> bool: + """Return whether the configuration at an index falls inside the cube.""" + return self.covers_values(domain.values_of(index)) + + def contains(self, other: MultiValueCube) -> bool: + """Return whether this cube covers everything ``other`` covers.""" + return all(theirs <= mine for mine, theirs in zip(self.pattern, other.pattern, strict=True)) + + def merge(self, other: MultiValueCube) -> MultiValueCube | None: + """Merge two cubes differing at exactly one condition. + + Returns ``None`` when they differ at none or several, in which case the + union would cover configurations neither cube covers. + """ + differing = [ + index + for index, (mine, theirs) in enumerate(zip(self.pattern, other.pattern, strict=True)) + if mine != theirs + ] + if len(differing) != 1: + return None + position = differing[0] + pattern = list(self.pattern) + pattern[position] = self.pattern[position] | other.pattern[position] + return MultiValueCube(tuple(pattern)) + + def as_expression(self, domain: MultiValueDomain) -> str: + """Render in multi-value QCA notation, for example ``A{0,2}*B{1}``. + + Conditions allowing every level are omitted, since they constrain + nothing. A cube constraining nothing renders as ``1``. + """ + parts = [ + f"{name}{{{','.join(str(value) for value in sorted(allowed))}}}" + for name, allowed, count in zip( + domain.conditions, self.pattern, domain.levels, strict=True + ) + if len(allowed) < count + ] + return "*".join(parts) if parts else "1" + + +@dataclass(frozen=True, slots=True) +class MultiValueSolution: + """A minimal cover of multi-value configurations.""" + + cubes: tuple[MultiValueCube, ...] + + def literal_count(self, domain: MultiValueDomain) -> int: + """Return the total number of constrained conditions.""" + return sum(cube.literals(domain) for cube in self.cubes) + + def as_expression(self, domain: MultiValueDomain) -> str: + """Render the whole cover, for example ``A{0}*B{1} + A{2}``.""" + return " + ".join(cube.as_expression(domain) for cube in self.cubes) + + def covers(self, index: int, domain: MultiValueDomain) -> bool: + """Return whether any cube covers a configuration.""" + return any(cube.covers(index, domain) for cube in self.cubes) + + +def prime_cubes( + on_set: set[int], dont_cares: set[int], domain: MultiValueDomain +) -> tuple[MultiValueCube, ...]: + """Generate every prime cube for a multi-value problem. + + Parameters + ---------- + on_set : set of int + Configuration indices that must be covered. + dont_cares : set of int + Configurations usable but not required. + domain : MultiValueDomain + The property space. + + Returns + ------- + tuple of MultiValueCube + Prime cubes, ordered by literal count then rendered form, and filtered + to those covering at least one required configuration. + + Raises + ------ + ValueError + If the two sets overlap. + """ + if on_set & dont_cares: + raise ValueError("on_set and dont_cares must be disjoint.") + universe = on_set | dont_cares + if not universe: + return () + + generated = {MultiValueCube.from_configuration(domain.values_of(index)) for index in universe} + frontier = set(generated) + + while frontier: + produced: set[MultiValueCube] = set() + ordered = sorted(frontier, key=lambda cube: tuple(sorted(sorted(s) for s in cube.pattern))) + for position, left in enumerate(ordered): + for right in ordered[position + 1 :]: + merged = left.merge(right) + if merged is not None and merged not in generated: + produced.add(merged) + generated |= produced + frontier = produced + + # A cube contained in another is not prime. Equality is excluded so that + # two identical cubes do not eliminate each other. + primes = [ + cube + for cube in generated + if not any(other != cube and other.contains(cube) for other in generated) + ] + useful = [cube for cube in primes if any(cube.covers(index, domain) for index in on_set)] + return tuple( + sorted(useful, key=lambda cube: (cube.literals(domain), cube.as_expression(domain))) + ) + + +def minimize_multivalue( + on_set: set[int], + *, + domain: MultiValueDomain, + dont_cares: set[int] | None = None, + max_solutions: int = 256, +) -> tuple[MultiValueSolution, ...]: + """Return every exact minimum cover of a multi-value problem. + + Parameters + ---------- + on_set : set of int + Configuration indices that must be covered. + domain : MultiValueDomain + The property space. + dont_cares : set of int, optional + Configurations usable but not required, typically logical remainders. + max_solutions : int, default 256 + Upper bound on tied minimum covers. + + Returns + ------- + tuple of MultiValueSolution + Every cover of provably minimal cost. + """ + required = set(on_set) + if not required: + return (MultiValueSolution(()),) + + primes = prime_cubes(required, set(dont_cares or ()), domain) + covered = [ + frozenset(index for index in required if cube.covers(index, domain)) for cube in primes + ] + literals = [cube.literals(domain) for cube in primes] + choices = solve_minimum_cover(covered, literals, required, max_solutions=max_solutions) + return tuple(MultiValueSolution(tuple(primes[i] for i in indices)) for indices in choices) diff --git a/src/setqca/multivalue/_domain.py b/src/setqca/multivalue/_domain.py new file mode 100644 index 0000000..06676dc --- /dev/null +++ b/src/setqca/multivalue/_domain.py @@ -0,0 +1,104 @@ +"""The property space of a multi-value model. + +A multi-value condition takes one of several unordered categories rather than +being present or absent. ``A{0,1,2}`` is a three-level condition, and a case +takes exactly one of those levels. + +Configurations are indexed in mixed radix, generalising the binary minterm. +With levels ``(2, 3)`` the index of ``(1, 2)`` is ``1 * 3 + 2 = 5``. Big-endian, +so the first condition is the most significant, matching the binary case. +""" + +from __future__ import annotations + +from collections.abc import Iterator, Mapping, Sequence +from dataclasses import dataclass +from itertools import product + + +@dataclass(frozen=True, slots=True) +class MultiValueDomain: + """Condition names and how many levels each takes. + + Parameters + ---------- + conditions : tuple of str + Condition names, in the order used for indexing. + levels : tuple of int + Number of categories per condition. Level values are ``0..levels-1``. + """ + + conditions: tuple[str, ...] + levels: tuple[int, ...] + + def __post_init__(self) -> None: + if not self.conditions: + raise ValueError("At least one condition is required.") + if len(self.conditions) != len(self.levels): + raise ValueError("conditions and levels must have the same length.") + if len(set(self.conditions)) != len(self.conditions): + raise ValueError("Condition names must be unique.") + if any(count < 2 for count in self.levels): + raise ValueError("Every condition needs at least two levels.") + + @classmethod + def from_mapping(cls, levels: Mapping[str, int]) -> MultiValueDomain: + """Build a domain from a ``{condition: levels}`` mapping.""" + return cls(tuple(levels), tuple(levels.values())) + + @property + def size(self) -> int: + """Return the number of logically possible configurations.""" + total = 1 + for count in self.levels: + total *= count + return total + + @property + def width(self) -> int: + """Return the number of conditions.""" + return len(self.conditions) + + def index_of(self, values: Sequence[int]) -> int: + """Return the mixed-radix index of one configuration. + + Raises + ------ + ValueError + If the length is wrong or a value is outside its condition's range. + """ + if len(values) != self.width: + raise ValueError(f"Expected {self.width} values, got {len(values)}.") + index = 0 + for value, count in zip(values, self.levels, strict=True): + if not 0 <= value < count: + raise ValueError(f"Value {value} is outside the range 0..{count - 1}.") + index = index * count + value + return index + + def values_of(self, index: int) -> tuple[int, ...]: + """Return the configuration at a mixed-radix index. + + Raises + ------ + ValueError + If the index is outside the property space. + """ + if not 0 <= index < self.size: + raise ValueError(f"Index {index} is outside the property space of size {self.size}.") + values: list[int] = [] + remaining = index + for count in reversed(self.levels): + values.append(remaining % count) + remaining //= count + return tuple(reversed(values)) + + def configurations(self) -> Iterator[tuple[int, ...]]: + """Yield every configuration, in index order.""" + yield from product(*(range(count) for count in self.levels)) + + def __str__(self) -> str: + return ", ".join( + f"{name}{{{','.join(str(value) for value in range(count))}}}" + for name, count in zip(self.conditions, self.levels, strict=True) + ) diff --git a/src/setqca/multivalue/_model.py b/src/setqca/multivalue/_model.py new file mode 100644 index 0000000..24cf172 --- /dev/null +++ b/src/setqca/multivalue/_model.py @@ -0,0 +1,369 @@ +"""Multi-value truth tables and the mvQCA estimator. + +The public shape mirrors :mod:`setqca.models` so that moving between csQCA, +fsQCA and mvQCA is a change of estimator rather than a change of workflow. + +Conditions are categorical and each case falls in exactly one configuration, so +configuration membership is crisp. The outcome may still be fuzzy: sufficiency +of a crisp configuration for a fuzzy outcome is well defined, and reduces to +the proportion of cases showing the outcome when the outcome is crisp too. +""" + +from __future__ import annotations + +from dataclasses import dataclass, field +from typing import TYPE_CHECKING, Any, Literal + +import numpy as np +import pandas as pd + +from setqca._validation import validate_columns, validate_membership +from setqca.metrics import SufficiencyFit, sufficiency + +from ._cube import MultiValueSolution, minimize_multivalue +from ._domain import MultiValueDomain + +if TYPE_CHECKING: # pragma: no cover - imported for type checking only + from collections.abc import Mapping + +MultiValueCode = Literal["1", "0", "C", "R"] + + +@dataclass(frozen=True, slots=True) +class MultiValueRow: + """One configuration of the multi-value property space.""" + + index: int + configuration: tuple[int, ...] + frequency: int + consistency: float + pri: float + outcome: MultiValueCode + cases: tuple[str, ...] + exclusion_reason: str | None = None + + @property + def observed(self) -> bool: + """Return whether the configuration passed the frequency cutoff.""" + return self.outcome != "R" + + +@dataclass(frozen=True, slots=True) +class MultiValueTruthTable: + """A complete multi-value truth table.""" + + domain: MultiValueDomain + outcome_name: str + rows: tuple[MultiValueRow, ...] + inclusion_cutoff: float + frequency_cutoff: int + + @property + def positive_indices(self) -> set[int]: + """Return configurations coded sufficient.""" + return {row.index for row in self.rows if row.outcome == "1"} + + @property + def remainder_indices(self) -> set[int]: + """Return configurations with too few cases to judge.""" + return {row.index for row in self.rows if row.outcome == "R"} + + def rows_with(self, code: MultiValueCode) -> tuple[MultiValueRow, ...]: + """Return rows carrying one outcome code.""" + return tuple(row for row in self.rows if row.outcome == code) + + def to_frame(self) -> pd.DataFrame: + """Return a tidy representation, one row per configuration.""" + records: list[dict[str, Any]] = [] + for row in self.rows: + record: dict[str, Any] = dict( + zip(self.domain.conditions, row.configuration, strict=True) + ) + record.update( + { + "index": row.index, + "n": row.frequency, + "consistency": row.consistency, + "PRI": row.pri, + "OUT": row.outcome, + "cases": ", ".join(row.cases), + "excluded_because": row.exclusion_reason or "", + } + ) + records.append(record) + return pd.DataFrame.from_records(records) + + def minimize( + self, *, include_remainders: bool = False, max_solutions: int = 256 + ) -> tuple[MultiValueSolution, ...]: + """Minimise directly from the table. + + Raises + ------ + ValueError + If no configuration is coded sufficient. + """ + on_set = self.positive_indices + if not on_set: + raise ValueError("No configuration is sufficient under the chosen thresholds.") + return minimize_multivalue( + on_set, + domain=self.domain, + dont_cares=self.remainder_indices if include_remainders else None, + max_solutions=max_solutions, + ) + + +@dataclass(frozen=True, slots=True) +class MultiValueResult: + """A fitted mvQCA analysis.""" + + domain: MultiValueDomain + outcome: str + truth_table: MultiValueTruthTable + conservative: tuple[MultiValueSolution, ...] + parsimonious: tuple[MultiValueSolution, ...] + fits: dict[str, SufficiencyFit] = field(default_factory=dict) + + def summary_frame(self, solution: str = "conservative") -> pd.DataFrame: + """Return one row per minimal solution of a family.""" + if solution not in ("conservative", "parsimonious"): + raise ValueError( + f"Unknown solution kind {solution!r}; expected 'conservative' or 'parsimonious'." + ) + solutions: tuple[MultiValueSolution, ...] = getattr(self, solution) + expressions = [item.as_expression(self.domain) for item in solutions] + return pd.DataFrame( + { + "solution": expressions, + "n_cubes": [len(item.cubes) for item in solutions], + "n_literals": [item.literal_count(self.domain) for item in solutions], + "consistency": [ + self.fits[expression].consistency if expression in self.fits else float("nan") + for expression in expressions + ], + "coverage": [ + self.fits[expression].coverage if expression in self.fits else float("nan") + for expression in expressions + ], + } + ) + + def __str__(self) -> str: + lines = [ + "Multi-value Qualitative Comparative Analysis", + f"Outcome: {self.outcome}", + f"Property space: {self.domain}", + f"Sufficient configurations: {len(self.truth_table.positive_indices)}", + "", + "Conservative solution(s):", + ] + lines.extend(f" {item.as_expression(self.domain)}" for item in self.conservative) + lines.append("Parsimonious solution(s):") + lines.extend(f" {item.as_expression(self.domain)}" for item in self.parsimonious) + return "\n".join(lines) + + +def _levels(data: pd.DataFrame, conditions: list[str]) -> tuple[int, ...]: + counts: list[int] = [] + for name in conditions: + values = data[name].to_numpy() + if not np.all(np.equal(np.mod(values, 1), 0)): + raise ValueError( + f"Condition {name!r} must hold integer category codes; mvQCA " + "conditions are categorical, not fuzzy." + ) + integers = values.astype(np.int64) + if integers.min() < 0: + raise ValueError(f"Condition {name!r} has a negative category code.") + counts.append(int(integers.max()) + 1) + return tuple(counts) + + +def build_multivalue_truth_table( + data: pd.DataFrame, + *, + outcome: str, + conditions: list[str] | tuple[str, ...], + levels: Mapping[str, int] | None = None, + inclusion_cutoff: float = 0.8, + frequency_cutoff: int = 1, + case_id: str | None = None, +) -> MultiValueTruthTable: + """Build a complete multi-value truth table. + + Parameters + ---------- + data : pandas.DataFrame + Condition columns holding integer category codes from ``0``, and an + outcome column holding memberships in ``[0, 1]``. + outcome : str + Name of the outcome column. + conditions : list of str or tuple of str + Condition columns. + levels : mapping of str to int, optional + Number of categories per condition. Inferred from the data when + omitted, which can understate a level that no case happens to take. + inclusion_cutoff : float, default 0.8 + Minimum sufficiency consistency for a configuration to count. + frequency_cutoff : int, default 1 + Minimum number of cases for a configuration to be observed. + case_id : str, optional + Column holding case labels. Defaults to the frame index. + + Returns + ------- + MultiValueTruthTable + One row per logically possible configuration. + + Raises + ------ + ValueError + If a condition is not categorical, a level count is too small, or a + cutoff is out of range. + """ + if not isinstance(data, pd.DataFrame): + raise TypeError("data must be a pandas DataFrame.") + names = validate_columns(data, conditions) + if not names: + raise ValueError("At least one condition is required.") + validate_columns(data, [outcome]) + if not 0.0 <= inclusion_cutoff <= 1.0: + raise ValueError("inclusion_cutoff must be in [0, 1].") + if frequency_cutoff < 1: + raise ValueError("frequency_cutoff must be at least 1.") + + counts = _levels(data, names) + if levels is not None: + missing = [name for name in names if name not in levels] + if missing: + raise KeyError(f"Levels missing for conditions: {missing}") + declared = tuple(int(levels[name]) for name in names) + for name, observed, stated in zip(names, counts, declared, strict=True): + if stated < observed: + raise ValueError( + f"Condition {name!r} declares {stated} levels but the data use {observed}." + ) + counts = declared + + domain = MultiValueDomain(tuple(names), counts) + y = validate_membership(data[outcome].to_numpy(), name=outcome) + codes = data[names].to_numpy().astype(np.int64) + + if case_id is None: + labels = np.asarray([str(index) for index in data.index], dtype=object) + else: + validate_columns(data, [case_id]) + labels = data[case_id].astype(str).to_numpy(dtype=object) + + rows: list[MultiValueRow] = [] + for configuration in domain.configurations(): + selector = np.all(codes == np.asarray(configuration), axis=1) + n = int(selector.sum()) + # Each case belongs to exactly one configuration, so membership is crisp. + membership = selector.astype(np.float64) + fit = sufficiency(membership, y) + + reason: str | None = None + if n < frequency_cutoff: + code: MultiValueCode = "R" + reason = f"frequency {n} below the cutoff of {frequency_cutoff}" + elif fit.consistency >= inclusion_cutoff: + code = "1" + else: + code = "0" + reason = ( + f"consistency {fit.consistency:.3f} below the inclusion cutoff " + f"of {inclusion_cutoff}" + ) + + rows.append( + MultiValueRow( + index=domain.index_of(configuration), + configuration=configuration, + frequency=n, + consistency=fit.consistency, + pri=fit.pri, + outcome=code, + cases=tuple(str(value) for value in labels[selector]), + exclusion_reason=reason, + ) + ) + + return MultiValueTruthTable( + domain=domain, + outcome_name=outcome, + rows=tuple(rows), + inclusion_cutoff=inclusion_cutoff, + frequency_cutoff=frequency_cutoff, + ) + + +@dataclass(slots=True) +class MVQCA: + """Multi-value Qualitative Comparative Analysis estimator. + + Parameters + ---------- + consistency : float, default 0.8 + Inclusion cutoff on sufficiency consistency. + frequency : int, default 1 + Minimum number of cases for a configuration to be observed. + max_solutions : int, default 256 + Upper bound on tied minimal covers. + levels : mapping of str to int, optional + Declared number of categories per condition. Supply this when a level + is theoretically possible but happens to have no cases, since it + changes the property space and therefore the remainders. + """ + + consistency: float = 0.8 + frequency: int = 1 + max_solutions: int = 256 + levels: Mapping[str, int] | None = None + + def fit( + self, + data: pd.DataFrame, + *, + outcome: str, + conditions: list[str] | tuple[str, ...], + case_id: str | None = None, + ) -> MultiValueResult: + """Fit mvQCA to categorical conditions and a calibrated outcome.""" + table = build_multivalue_truth_table( + data, + outcome=outcome, + conditions=conditions, + levels=self.levels, + inclusion_cutoff=self.consistency, + frequency_cutoff=self.frequency, + case_id=case_id, + ) + conservative = table.minimize(max_solutions=self.max_solutions) + parsimonious = table.minimize(include_remainders=True, max_solutions=self.max_solutions) + + y = validate_membership(data[outcome].to_numpy(), name=outcome) + codes = data[list(table.domain.conditions)].to_numpy().astype(np.int64) + fits: dict[str, SufficiencyFit] = {} + for solution in (*conservative, *parsimonious): + expression = solution.as_expression(table.domain) + if expression in fits: + continue + membership = np.asarray( + [ + 1.0 if solution.covers(table.domain.index_of(tuple(row)), table.domain) else 0.0 + for row in codes + ], + dtype=np.float64, + ) + fits[expression] = sufficiency(membership, y) + + return MultiValueResult( + domain=table.domain, + outcome=outcome, + truth_table=table, + conservative=conservative, + parsimonious=parsimonious, + fits=fits, + ) diff --git a/tests/test_multivalue.py b/tests/test_multivalue.py new file mode 100644 index 0000000..b3dc815 --- /dev/null +++ b/tests/test_multivalue.py @@ -0,0 +1,442 @@ +"""Tests for multi-value QCA. + +The central claim under test is that the multi-value engine minimises cubes +directly rather than reducing to Boolean indicators, and that its answers are +exact. +""" + +from __future__ import annotations + +from itertools import combinations, product + +import pandas as pd +import pytest + +from setqca.multivalue import ( + MVQCA, + MultiValueCube, + MultiValueDomain, + build_multivalue_truth_table, + minimize_multivalue, + prime_cubes, +) + +# Two conditions: a three-level regime type and a binary wealth indicator. +DOMAIN = MultiValueDomain(("regime", "wealth"), (3, 2)) + + +class TestDomain: + def test_the_property_space_is_the_product_of_the_levels(self) -> None: + assert DOMAIN.size == 6 + assert len(list(DOMAIN.configurations())) == 6 + + def test_indexing_is_mixed_radix_and_big_endian(self) -> None: + assert DOMAIN.index_of((0, 0)) == 0 + assert DOMAIN.index_of((0, 1)) == 1 + assert DOMAIN.index_of((1, 0)) == 2 + assert DOMAIN.index_of((2, 1)) == 5 + + def test_indices_round_trip(self) -> None: + for configuration in DOMAIN.configurations(): + assert DOMAIN.values_of(DOMAIN.index_of(configuration)) == configuration + + def test_a_binary_domain_reproduces_the_minterm_encoding(self) -> None: + """With two levels everywhere, mixed radix is the ordinary minterm.""" + binary = MultiValueDomain(("A", "B", "C"), (2, 2, 2)) + assert binary.index_of((1, 1, 0)) == 6 + assert binary.values_of(5) == (1, 0, 1) + + def test_a_value_outside_its_range_is_rejected(self) -> None: + with pytest.raises(ValueError, match="outside the range"): + DOMAIN.index_of((3, 0)) + + def test_an_index_outside_the_space_is_rejected(self) -> None: + with pytest.raises(ValueError, match="outside the property space"): + DOMAIN.values_of(6) + + def test_the_wrong_number_of_values_is_rejected(self) -> None: + with pytest.raises(ValueError, match="Expected 2 values"): + DOMAIN.index_of((0,)) + + @pytest.mark.parametrize( + ("conditions", "levels", "message"), + [ + ((), (), "At least one condition"), + (("A",), (2, 3), "same length"), + (("A", "A"), (2, 2), "unique"), + (("A",), (1,), "at least two levels"), + ], + ) + def test_malformed_domains_are_rejected( + self, conditions: tuple[str, ...], levels: tuple[int, ...], message: str + ) -> None: + with pytest.raises(ValueError, match=message): + MultiValueDomain(conditions, levels) + + def test_a_domain_renders_in_multi_value_notation(self) -> None: + assert str(DOMAIN) == "regime{0,1,2}, wealth{0,1}" + + def test_a_domain_can_be_built_from_a_mapping(self) -> None: + assert MultiValueDomain.from_mapping({"regime": 3, "wealth": 2}) == DOMAIN + + +class TestCubes: + def test_a_configuration_cube_covers_only_itself(self) -> None: + cube = MultiValueCube.from_configuration((1, 0)) + assert cube.covers(DOMAIN.index_of((1, 0)), DOMAIN) + assert not cube.covers(DOMAIN.index_of((2, 0)), DOMAIN) + + def test_merging_takes_the_union_at_the_one_differing_condition(self) -> None: + left = MultiValueCube.from_configuration((0, 1)) + right = MultiValueCube.from_configuration((2, 1)) + merged = left.merge(right) + assert merged is not None + assert merged.pattern[0] == frozenset({0, 2}) + assert merged.pattern[1] == frozenset({1}) + + def test_a_merged_cube_covers_exactly_the_union(self) -> None: + """The property the merge rule rests on: no extra coverage appears.""" + left = MultiValueCube.from_configuration((0, 1)) + right = MultiValueCube.from_configuration((2, 1)) + merged = left.merge(right) + assert merged is not None + covered = {i for i in range(DOMAIN.size) if merged.covers(i, DOMAIN)} + assert covered == {DOMAIN.index_of((0, 1)), DOMAIN.index_of((2, 1))} + + def test_cubes_differing_at_two_conditions_do_not_merge(self) -> None: + left = MultiValueCube.from_configuration((0, 0)) + right = MultiValueCube.from_configuration((1, 1)) + assert left.merge(right) is None + + def test_identical_cubes_do_not_merge(self) -> None: + cube = MultiValueCube.from_configuration((0, 0)) + assert cube.merge(cube) is None + + def test_a_condition_allowing_every_level_is_not_a_literal(self) -> None: + cube = MultiValueCube((frozenset({0, 1, 2}), frozenset({1}))) + assert cube.literals(DOMAIN) == 1 + assert cube.as_expression(DOMAIN) == "wealth{1}" + + def test_a_cube_constraining_nothing_renders_as_one(self) -> None: + cube = MultiValueCube((frozenset({0, 1, 2}), frozenset({0, 1}))) + assert cube.is_tautology(DOMAIN) + assert cube.as_expression(DOMAIN) == "1" + + def test_containment_is_per_condition(self) -> None: + wide = MultiValueCube((frozenset({0, 1}), frozenset({0, 1}))) + narrow = MultiValueCube((frozenset({0}), frozenset({1}))) + assert wide.contains(narrow) + assert not narrow.contains(wide) + + +class TestMinimisation: + def test_all_levels_of_a_condition_collapse_it(self) -> None: + """Regime taking every level, with wealth fixed, eliminates regime.""" + on_set = {DOMAIN.index_of((level, 1)) for level in range(3)} + solutions = minimize_multivalue(on_set, domain=DOMAIN) + assert len(solutions) == 1 + assert solutions[0].as_expression(DOMAIN) == "wealth{1}" + + def test_a_partial_set_of_levels_survives_as_a_subset_literal(self) -> None: + """Two of three levels cannot eliminate the condition, but do combine.""" + on_set = {DOMAIN.index_of((0, 1)), DOMAIN.index_of((2, 1))} + solutions = minimize_multivalue(on_set, domain=DOMAIN) + assert len(solutions) == 1 + assert solutions[0].as_expression(DOMAIN) == "regime{0,2}*wealth{1}" + + def test_a_single_configuration_is_its_own_solution(self) -> None: + on_set = {DOMAIN.index_of((1, 0))} + solutions = minimize_multivalue(on_set, domain=DOMAIN) + assert solutions[0].as_expression(DOMAIN) == "regime{1}*wealth{0}" + + def test_an_empty_problem_yields_the_empty_cover(self) -> None: + assert minimize_multivalue(set(), domain=DOMAIN)[0].cubes == () + + def test_a_problem_with_nothing_at_all_yields_no_primes(self) -> None: + assert prime_cubes(set(), set(), DOMAIN) == () + + def test_rows_know_whether_they_were_observed(self) -> None: + table = build_multivalue_truth_table( + pd.DataFrame({"regime": [0, 1], "wealth": [0, 1], "Y": [0.1, 0.9]}), + outcome="Y", + conditions=["regime", "wealth"], + ) + observed = [row for row in table.rows if row.observed] + assert len(observed) == 2 + assert all(row.frequency > 0 for row in observed) + + def test_remainders_permit_a_simpler_cover(self) -> None: + on_set = {DOMAIN.index_of((0, 1)), DOMAIN.index_of((1, 1))} + remainder = {DOMAIN.index_of((2, 1))} + conservative = minimize_multivalue(on_set, domain=DOMAIN)[0] + parsimonious = minimize_multivalue(on_set, domain=DOMAIN, dont_cares=remainder)[0] + + assert parsimonious.literal_count(DOMAIN) < conservative.literal_count(DOMAIN) + assert parsimonious.as_expression(DOMAIN) == "wealth{1}" + + def test_overlapping_sets_are_rejected(self) -> None: + with pytest.raises(ValueError, match="disjoint"): + prime_cubes({0, 1}, {1}, DOMAIN) + + def test_every_returned_cover_is_valid_and_tied_on_cost(self) -> None: + on_set = {DOMAIN.index_of((0, 0)), DOMAIN.index_of((1, 1)), DOMAIN.index_of((2, 0))} + off_set = set(range(DOMAIN.size)) - on_set + solutions = minimize_multivalue(on_set, domain=DOMAIN) + + costs = {(len(item.cubes), item.literal_count(DOMAIN)) for item in solutions} + assert len(costs) == 1 + for solution in solutions: + covered = {i for i in range(DOMAIN.size) if solution.covers(i, DOMAIN)} + assert on_set <= covered + assert not covered & off_set + + +class TestExactness: + """Verified against exhaustive search, as the binary engine is.""" + + @pytest.mark.parametrize("levels", [(2, 2), (3, 2), (2, 3), (3, 3)]) + def test_minimisation_matches_brute_force(self, levels: tuple[int, ...]) -> None: + domain = MultiValueDomain(tuple("AB"[: len(levels)]), levels) + universe = list(range(domain.size)) + + # Every legal cube: a non-empty subset of levels per condition. + options = [ + [ + frozenset(subset) + for size in range(1, count + 1) + for subset in combinations(range(count), size) + ] + for count in levels + ] + all_cubes = [MultiValueCube(tuple(pattern)) for pattern in product(*options)] + + for size in (1, 2, 3): + for on_tuple in combinations(universe, size): + on_set = set(on_tuple) + off_set = set(universe) - on_set + + legal = [ + cube + for cube in all_cubes + if not any(cube.covers(i, domain) for i in off_set) + and any(cube.covers(i, domain) for i in on_set) + ] + expected: tuple[int, int] | None = None + for count in range(1, len(legal) + 1): + candidates = [ + (count, sum(cube.literals(domain) for cube in choice)) + for choice in combinations(legal, count) + if on_set + <= {i for i in universe for cube in choice if cube.covers(i, domain)} + ] + if candidates: + expected = min(candidates) + break + + solution = minimize_multivalue(on_set, domain=domain)[0] + obtained = (len(solution.cubes), solution.literal_count(domain)) + assert obtained == expected, f"levels={levels} on_set={sorted(on_set)}" + + +class TestNotBooleanDummies: + """The encoding shortcut this implementation deliberately avoids.""" + + def test_a_three_level_condition_is_not_split_into_indicators(self) -> None: + on_set = {DOMAIN.index_of((0, 1)), DOMAIN.index_of((2, 1))} + expression = minimize_multivalue(on_set, domain=DOMAIN)[0].as_expression(DOMAIN) + + # A dummy encoding would express this through indicator variables such + # as regime_0 and regime_2; the multi-value form keeps one literal. + assert "regime{0,2}" in expression + assert "regime_0" not in expression + + def test_no_cube_can_describe_an_impossible_case(self) -> None: + """Every cube covers only real configurations, by construction. + + A Boolean dummy encoding admits points where two indicators for the + same condition are both true, which correspond to no case at all. + """ + on_set = set(range(DOMAIN.size)) + for cube in prime_cubes(on_set, set(), DOMAIN): + covered = [i for i in range(DOMAIN.size) if cube.covers(i, DOMAIN)] + for index in covered: + values = DOMAIN.values_of(index) + assert len(values) == DOMAIN.width + for value, count in zip(values, DOMAIN.levels, strict=True): + assert 0 <= value < count + + +class TestTruthTable: + data = pd.DataFrame( + { + "regime": [0, 0, 1, 1, 2, 2], + "wealth": [0, 0, 1, 1, 1, 1], + "Y": [0.1, 0.2, 0.9, 0.95, 0.85, 0.9], + }, + index=["a", "b", "c", "d", "e", "f"], + ) + + def test_the_table_covers_the_whole_property_space(self) -> None: + table = build_multivalue_truth_table( + self.data, outcome="Y", conditions=["regime", "wealth"] + ) + assert len(table.rows) == 6 + assert table.domain.levels == (3, 2) + + def test_unobserved_configurations_are_remainders(self) -> None: + table = build_multivalue_truth_table( + self.data, outcome="Y", conditions=["regime", "wealth"] + ) + assert table.remainder_indices + for row in table.rows_with("R"): + assert row.frequency == 0 + assert "frequency" in (row.exclusion_reason or "") + + def test_consistency_is_the_share_of_the_outcome_in_the_configuration(self) -> None: + table = build_multivalue_truth_table( + self.data, outcome="Y", conditions=["regime", "wealth"] + ) + row = next(row for row in table.rows if row.configuration == (1, 1)) + assert row.frequency == 2 + assert row.consistency == pytest.approx((0.9 + 0.95) / 2) + + def test_declared_levels_enlarge_the_property_space(self) -> None: + """A level with no cases still exists, and becomes a remainder.""" + table = build_multivalue_truth_table( + self.data, + outcome="Y", + conditions=["regime", "wealth"], + levels={"regime": 4, "wealth": 2}, + ) + assert table.domain.levels == (4, 2) + assert len(table.rows) == 8 + + def test_declared_levels_may_not_contradict_the_data(self) -> None: + with pytest.raises(ValueError, match="declares 2 levels"): + build_multivalue_truth_table( + self.data, + outcome="Y", + conditions=["regime", "wealth"], + levels={"regime": 2, "wealth": 2}, + ) + + def test_missing_declared_levels_are_reported(self) -> None: + with pytest.raises(KeyError, match="Levels missing"): + build_multivalue_truth_table( + self.data, + outcome="Y", + conditions=["regime", "wealth"], + levels={"regime": 3}, + ) + + def test_fuzzy_conditions_are_rejected(self) -> None: + frame = pd.DataFrame({"regime": [0.5, 1.0], "Y": [0.9, 0.1]}) + with pytest.raises(ValueError, match="categorical, not fuzzy"): + build_multivalue_truth_table(frame, outcome="Y", conditions=["regime"]) + + def test_negative_category_codes_are_rejected(self) -> None: + frame = pd.DataFrame({"regime": [-1, 1], "Y": [0.9, 0.1]}) + with pytest.raises(ValueError, match="negative category code"): + build_multivalue_truth_table(frame, outcome="Y", conditions=["regime"]) + + def test_the_frame_export_names_the_conditions(self) -> None: + frame = build_multivalue_truth_table( + self.data, outcome="Y", conditions=["regime", "wealth"] + ).to_frame() + assert list(frame.columns)[:2] == ["regime", "wealth"] + assert "OUT" in frame.columns + + def test_a_table_with_nothing_sufficient_refuses_to_minimise(self) -> None: + table = build_multivalue_truth_table( + self.data, outcome="Y", conditions=["regime", "wealth"], inclusion_cutoff=1.0 + ) + with pytest.raises(ValueError, match="No configuration is sufficient"): + table.minimize() + + @pytest.mark.parametrize( + ("kwargs", "message"), + [ + ({"inclusion_cutoff": 1.5}, "inclusion_cutoff"), + ({"frequency_cutoff": 0}, "frequency_cutoff"), + ({"conditions": []}, "At least one condition"), + ], + ) + def test_guards(self, kwargs: dict[str, object], message: str) -> None: + arguments: dict[str, object] = { + "outcome": "Y", + "conditions": ["regime", "wealth"], + } + arguments.update(kwargs) + with pytest.raises(ValueError, match=message): + build_multivalue_truth_table(self.data, **arguments) # type: ignore[arg-type] + + def test_data_must_be_a_frame(self) -> None: + with pytest.raises(TypeError, match="pandas DataFrame"): + build_multivalue_truth_table({"a": [1]}, outcome="Y", conditions=["a"]) # type: ignore[arg-type] + + +class TestEstimator: + data = TestTruthTable.data + + def test_the_estimator_returns_both_families(self) -> None: + result = MVQCA(consistency=0.8).fit(self.data, outcome="Y", conditions=["regime", "wealth"]) + assert result.conservative + assert result.parsimonious + + def test_the_parsimonious_solution_is_no_more_complex(self) -> None: + result = MVQCA(consistency=0.8).fit(self.data, outcome="Y", conditions=["regime", "wealth"]) + assert min(item.literal_count(result.domain) for item in result.parsimonious) <= min( + item.literal_count(result.domain) for item in result.conservative + ) + + def test_the_summary_frame_carries_fit(self) -> None: + result = MVQCA(consistency=0.8).fit(self.data, outcome="Y", conditions=["regime", "wealth"]) + frame = result.summary_frame("parsimonious") + assert list(frame.columns) == [ + "solution", + "n_cubes", + "n_literals", + "consistency", + "coverage", + ] + assert (frame["consistency"] >= 0.8).all() + + def test_an_unknown_family_is_rejected(self) -> None: + result = MVQCA().fit(self.data, outcome="Y", conditions=["regime", "wealth"]) + with pytest.raises(ValueError, match="Unknown solution kind"): + result.summary_frame("intermediate") + + def test_the_report_names_the_property_space(self) -> None: + result = MVQCA().fit(self.data, outcome="Y", conditions=["regime", "wealth"]) + text = str(result) + assert "Multi-value" in text + assert "regime{0,1,2}" in text + + def test_case_labels_can_come_from_a_column(self) -> None: + frame = self.data.reset_index(names="country") + result = MVQCA().fit(frame, outcome="Y", conditions=["regime", "wealth"], case_id="country") + assert any(row.cases for row in result.truth_table.rows) + + def test_declared_levels_flow_through_the_estimator(self) -> None: + result = MVQCA(levels={"regime": 4, "wealth": 2}).fit( + self.data, outcome="Y", conditions=["regime", "wealth"] + ) + assert result.domain.levels == (4, 2) + + +class TestBinaryAgreement: + """With two levels everywhere, mvQCA must reproduce csQCA.""" + + def test_a_binary_problem_gives_the_same_minimal_cost(self) -> None: + from setqca.minimize import minimize + + domain = MultiValueDomain(("A", "B", "C"), (2, 2, 2)) + for on_tuple in combinations(range(8), 3): + on_set = set(on_tuple) + binary = minimize(on_set, width=3)[0] + multi = minimize_multivalue(on_set, domain=domain)[0] + + assert (len(multi.cubes), multi.literal_count(domain)) == ( + len(binary.implicants), + binary.literal_count, + ) diff --git a/tests/test_parity.py b/tests/test_parity.py index ca20a67..f32dbf1 100644 --- a/tests/test_parity.py +++ b/tests/test_parity.py @@ -33,6 +33,7 @@ sufficiency, sufficiency_diagnostics, ) +from setqca.multivalue import MVQCA, MultiValueDomain, build_multivalue_truth_table pytestmark = pytest.mark.parity @@ -295,6 +296,87 @@ def test_parsimonious_solution_matches_r(analysis: dict[str, Any]) -> None: assert obtained == _canonical_solutions(analysis["parsimonious"]) +# --------------------------------------------------------------------------- +# Multi-value QCA +# --------------------------------------------------------------------------- + +_MULTIVALUE = FIXTURE.get("multivalue", []) + + +def _r_term_coverage(term: str, domain: MultiValueDomain) -> set[int]: + """Return the configurations an R multi-value term such as ``A[1]*B[2]`` covers.""" + fixed: dict[str, int] = {} + for literal in term.split("*"): + name, _, rest = literal.partition("[") + fixed[name.strip()] = int(rest.rstrip("]")) + return { + domain.index_of(configuration) + for configuration in domain.configurations() + if all( + configuration[domain.conditions.index(name)] == value for name, value in fixed.items() + ) + } + + +@pytest.mark.parametrize("case", _MULTIVALUE, ids=_ids(_MULTIVALUE)) +def test_multivalue_truth_table_matches_r(case: dict[str, Any]) -> None: + frame = pd.DataFrame(case["data"]) + table = build_multivalue_truth_table( + frame, + outcome=case["outcome"], + conditions=list(case["conditions"]), + inclusion_cutoff=float(case["incl_cut"]), + ) + obtained = {row.configuration: row for row in table.rows} + + assert len(obtained) == len(case["truth_table"]) + for expected in case["truth_table"]: + row = obtained[tuple(expected["configuration"])] + assert row.frequency == expected["n"] + assert row.outcome == expected["out"] + if expected["consistency"] is not None: + assert row.consistency == pytest.approx(expected["consistency"], abs=TOLERANCE) + + +@pytest.mark.parametrize("case", _MULTIVALUE, ids=_ids(_MULTIVALUE)) +@pytest.mark.parametrize("family", ["conservative", "parsimonious"]) +def test_multivalue_solutions_match_r_in_cost_and_coverage( + case: dict[str, Any], family: str +) -> None: + """Compared by cost and coverage, not by text. + + R writes single-value literals only (``regime[1]``); setqca also forms + subset literals (``regime{1,2}``) when they are prime. On the regime/wealth + benchmark R's conservative term ``regime[1]*wealth[1]`` is a proper subset + of setqca's ``regime{1,2}*wealth{1}``, so R's is not a prime implicant. Both + covers are minimal and cover the same configurations, which is what the + comparison asserts. + """ + expected_solutions = case[family] + if not expected_solutions: + pytest.skip("R produced no solution for this family") + + frame = pd.DataFrame(case["data"]) + result = MVQCA(consistency=float(case["incl_cut"])).fit( + frame, outcome=case["outcome"], conditions=list(case["conditions"]) + ) + domain = result.domain + obtained: tuple[Any, ...] = getattr(result, family) + + r_terms = _as_list(expected_solutions[0]) + r_coverage: set[int] = set() + for term in r_terms: + r_coverage |= _r_term_coverage(term, domain) + r_literals = sum(len(term.split("*")) for term in r_terms) + + ours = obtained[0] + our_coverage = {i for i in range(domain.size) if ours.covers(i, domain)} + + assert our_coverage == r_coverage, "the two solutions must select the same configurations" + assert len(ours.cubes) == len(r_terms), "the number of terms must agree" + assert ours.literal_count(domain) == r_literals, "the literal count must agree" + + # --------------------------------------------------------------------------- # Per-term fit, including unique coverage # --------------------------------------------------------------------------- diff --git a/validation/fixtures/r_qca.json b/validation/fixtures/r_qca.json index 2f6c2dd..2617f13 100644 --- a/validation/fixtures/r_qca.json +++ b/validation/fixtures/r_qca.json @@ -1313,5 +1313,135 @@ } ] } + ], + "multivalue": [ + { + "id": "regime-wealth", + "outcome": "Y", + "conditions": ["regime", "wealth"], + "incl_cut": 0.8, + "data": { + "regime": [0, 0, 1, 1, 2, 2, 1, 2], + "wealth": [0, 0, 1, 1, 1, 1, 0, 0], + "Y": [0, 0, 1, 1, 1, 1, 0, 1] + }, + "truth_table": [ + { + "configuration": [0, 0], + "n": 2, + "consistency": 0, + "out": "0" + }, + { + "configuration": [0, 1], + "n": 0, + "consistency": null, + "out": "R" + }, + { + "configuration": [1, 0], + "n": 1, + "consistency": 0, + "out": "0" + }, + { + "configuration": [1, 1], + "n": 2, + "consistency": 1, + "out": "1" + }, + { + "configuration": [2, 0], + "n": 1, + "consistency": 1, + "out": "1" + }, + { + "configuration": [2, 1], + "n": 2, + "consistency": 1, + "out": "1" + } + ], + "conservative": [ + ["regime[2]", "regime[1]*wealth[1]"] + ], + "parsimonious": [ + ["regime[2]", "wealth[1]"] + ] + }, + { + "id": "three-by-three", + "outcome": "Y", + "conditions": ["A", "B"], + "incl_cut": 0.8, + "data": { + "A": [0, 1, 2, 0, 1, 2, 0, 1, 2], + "B": [0, 0, 0, 1, 1, 1, 2, 2, 2], + "Y": [0, 0, 1, 0, 1, 1, 1, 1, 1] + }, + "truth_table": [ + { + "configuration": [0, 0], + "n": 1, + "consistency": 0, + "out": "0" + }, + { + "configuration": [0, 1], + "n": 1, + "consistency": 0, + "out": "0" + }, + { + "configuration": [0, 2], + "n": 1, + "consistency": 1, + "out": "1" + }, + { + "configuration": [1, 0], + "n": 1, + "consistency": 0, + "out": "0" + }, + { + "configuration": [1, 1], + "n": 1, + "consistency": 1, + "out": "1" + }, + { + "configuration": [1, 2], + "n": 1, + "consistency": 1, + "out": "1" + }, + { + "configuration": [2, 0], + "n": 1, + "consistency": 1, + "out": "1" + }, + { + "configuration": [2, 1], + "n": 1, + "consistency": 1, + "out": "1" + }, + { + "configuration": [2, 2], + "n": 1, + "consistency": 1, + "out": "1" + } + ], + "conservative": [ + ["A[2]", "B[2]", "A[1]*B[1]"] + ], + "parsimonious": [ + ["A[2]", "B[2]", "A[1]*B[1]"] + ] + } ] } diff --git a/validation/r/generate_fixtures.R b/validation/r/generate_fixtures.R index 2d969c1..ec0bb2d 100644 --- a/validation/r/generate_fixtures.R +++ b/validation/r/generate_fixtures.R @@ -352,6 +352,61 @@ calibration_cases <- list( ) ) +# --------------------------------------------------------------------------- +# Multi-value QCA +# --------------------------------------------------------------------------- + +# R writes multi-value literals as `regime[2]`; setqca writes `regime{2}` and +# additionally forms subset literals such as `regime{1,2}`, which R does not. +# The truth table and the parsimonious solution are directly comparable. +build_multivalue <- function(id, frame, outcome, conditions, incl_cut) { + tt <- truthTable(frame, + outcome = outcome, conditions = conditions, + incl.cut = incl_cut, show.cases = TRUE + ) + rows <- lapply(seq_len(nrow(tt$tt)), function(i) { + row <- tt$tt[i, ] + incl <- suppressWarnings(as.numeric(as.character(row[["incl"]]))) + list( + configuration = as.integer(unlist(row[conditions])), + n = as.integer(row[["n"]]), + consistency = if (is.na(incl)) NULL else incl, + out = recode_out(row[["OUT"]]) + ) + }) + list( + id = id, + outcome = outcome, + conditions = conditions, + incl_cut = incl_cut, + data = as.list(frame), + truth_table = rows, + conservative = safe_minimize(tt, include = ""), + parsimonious = safe_minimize(tt, include = "?") + ) +} + +multivalue_cases <- list( + build_multivalue( + "regime-wealth", + data.frame( + regime = c(0, 0, 1, 1, 2, 2, 1, 2), + wealth = c(0, 0, 1, 1, 1, 1, 0, 0), + Y = c(0, 0, 1, 1, 1, 1, 0, 1) + ), + "Y", c("regime", "wealth"), 0.8 + ), + build_multivalue( + "three-by-three", + data.frame( + A = c(0, 1, 2, 0, 1, 2, 0, 1, 2), + B = c(0, 0, 0, 1, 1, 1, 2, 2, 2), + Y = c(0, 0, 1, 0, 1, 1, 1, 1, 1) + ), + "Y", c("A", "B"), 0.8 + ) +) + # --------------------------------------------------------------------------- # Emit # --------------------------------------------------------------------------- @@ -374,7 +429,8 @@ fixture <- list( pof = pof_cases, analyses = analyses, intermediate = intermediate_cases, - necessity_screens = necessity_screens + necessity_screens = necessity_screens, + multivalue = multivalue_cases ) dir.create(dirname(out_path), recursive = TRUE, showWarnings = FALSE)