ml_meta_perf.terms

  1"""The term library: the vocabulary the equations are allowed to use.
  2
  3A term is a small, named, serialisable expression over one to three raw features --
  4``f``, ``log(f)``, ``1/f``, ``f^2``, ``f1/f2``, ``(f1+f2)/f3`` and friends. Terms are
  5built once from the training columns, evaluated into a design matrix, and the selected
  6ones are what the published equation prints. Keeping them as data (rather than as
  7closures) is what lets a fitted equation round-trip through JSON.
  8
  9The raw features span wildly different magnitudes -- ``gravity`` runs from 1.6 to 1e16 --
 10so composite terms are built over log-compressed operands wherever the feature is
 11strictly positive. Dividing raw ``gravity`` by anything produces a term whose weight is
 121e-16 and whose meaning is unreadable.
 13
 14**This vocabulary assumes strictly positive, continuous features, and degrades quietly
 15when given anything else.** The assumption is worth stating because nothing here raises
 16when it is violated -- the affected terms are simply never generated, and a caller adding
 17a feature sees a smaller library rather than an error. For a feature that reaches zero,
 18``log``, ``sqrt`` and ``1/f`` are all undefined (`Atom.is_defined_on`) and it can never be
 19a denominator (`denominator_atom`), which already excludes ``nr_norm``, ``nr_bin`` and
 20``nr_outliers`` from three quarters of the grammar. For a **binary** feature only ``f`` and
 21``f^2`` survive, and those are the same column twice -- see `Library`. What is left to a
 220/1 column is addition and multiplication, and the two are not equally safe: inside a sum
 23it contributes a level to every row, while inside a product it zeroes the term entirely on
 24the rows where it is off, which is a per-group slope rather than a relationship. Prefer a
 25continuous descriptor that grades the same distinction.
 26
 27Study chapter: [2. The additive model][study-chapter] -- the rationale, in
 28prose, with the figures.
 29
 30[study-chapter]: https://github.com/mariolpantunes/ml-meta-perf/blob/main/assets/docs/02-additive-model.md
 31"""
 32
 33from __future__ import annotations
 34
 35import itertools
 36from dataclasses import dataclass
 37from typing import Literal
 38
 39import numpy as np
 40
 41Transform = Literal["id", "log", "sqrt", "inv", "sq"]
 42Operation = Literal["atom", "ratio", "product", "sum_ratio", "ratio_of_sums"]
 43
 44DENOMINATOR_FLOOR = 1e-9
 45MAX_ABS_ZSCORE = 8.0
 46MIN_RELATIVE_SPREAD = 1e-6
 47MAX_DENOMINATOR_RANGE = 20.0
 48
 49#: Two terms count as the same column when their centred correlation reaches this. At
 50#: ``1 - 1e-9`` the pair is indistinguishable to the solver anyway, so the threshold
 51#: separates "algebraically identical" from "merely similar" rather than imposing a taste.
 52COLLINEARITY_TOLERANCE = 1.0 - 1e-9
 53
 54#: Every elementary transform a unary term may use.
 55TRANSFORMS: tuple[Transform, ...] = ("id", "log", "sqrt", "inv", "sq")
 56
 57_TRANSFORM_FORMAT: dict[Transform, str] = {
 58    "id": "{0}",
 59    "log": "log({0})",
 60    "sqrt": "sqrt({0})",
 61    "inv": "1/{0}",
 62    "sq": "{0}^2",
 63}
 64
 65
 66@dataclass(frozen=True)
 67class Atom:
 68    """One raw feature under one elementary transform."""
 69
 70    feature: str
 71    transform: Transform = "id"
 72
 73    @property
 74    def name(self) -> str:
 75        return _TRANSFORM_FORMAT[self.transform].format(self.feature)
 76
 77    def evaluate(self, columns: dict[str, np.ndarray]) -> np.ndarray:
 78        values = columns[self.feature]
 79        match self.transform:
 80            case "id":
 81                return values
 82            case "log":
 83                return np.log(values)
 84            case "sqrt":
 85                return np.sqrt(values)
 86            case "inv":
 87                return 1.0 / values
 88            case "sq":
 89                return values * values
 90
 91    def is_defined_on(self, columns: dict[str, np.ndarray]) -> bool:
 92        """Whether the transform is real-valued and finite over these columns."""
 93        values = columns.get(self.feature)
 94        if values is None:
 95            return False
 96        if self.transform in ("log", "sqrt", "inv"):
 97            return bool(np.all(values > 0.0))
 98        return True
 99
100
101@dataclass(frozen=True)
102class Term:
103    """A named expression over one to three operands.
104
105    An operand is an `Atom` or, when the term is nested, another `Term`. Nesting lets one
106    term carry a deep expression -- ``([log(a)] / [b]) * [log(c)]`` -- which matters
107    because the equation's length is counted in *terms*, not in operations. A short
108    equation over deep terms can be easier to read than a long one over shallow terms,
109    and ``depth`` is what keeps that trade under control.
110    """
111
112    operation: Operation
113    operands: tuple[Atom | Term, ...]
114
115    @property
116    def name(self) -> str:
117        parts = [operand.name for operand in self.operands]
118        match self.operation:
119            case "atom":
120                return parts[0]
121            case "ratio":
122                return f"[{parts[0]}] / [{parts[1]}]"
123            case "product":
124                return f"[{parts[0]}] * [{parts[1]}]"
125            case "sum_ratio":
126                return f"([{parts[0]}] + [{parts[1]}]) / [{parts[2]}]"
127            case "ratio_of_sums":
128                return f"([{parts[0]}] + [{parts[1]}]) / ([{parts[2]}] + [{parts[3]}])"
129
130    @property
131    def features(self) -> tuple[str, ...]:
132        collected: list[str] = []
133        for operand in self.operands:
134            if isinstance(operand, Term):
135                collected.extend(operand.features)
136            else:
137                collected.append(operand.feature)
138        return tuple(collected)
139
140    @property
141    def depth(self) -> int:
142        """1 for a term over atoms, one more for each level of nesting."""
143        nested = [operand.depth for operand in self.operands if isinstance(operand, Term)]
144        return 1 + max(nested, default=0)
145
146    def evaluate(self, columns: dict[str, np.ndarray]) -> np.ndarray:
147        values = [operand.evaluate(columns) for operand in self.operands]
148        match self.operation:
149            case "atom":
150                return values[0]
151            case "ratio":
152                return values[0] / _guard(values[1])
153            case "product":
154                return values[0] * values[1]
155            case "sum_ratio":
156                return (values[0] + values[1]) / _guard(values[2])
157            case "ratio_of_sums":
158                return (values[0] + values[1]) / _guard(values[2] + values[3])
159
160    def to_dict(self) -> dict[str, object]:
161        return {
162            "operation": self.operation,
163            "operands": [
164                operand.to_dict()
165                if isinstance(operand, Term)
166                else {"feature": operand.feature, "transform": operand.transform}
167                for operand in self.operands
168            ],
169        }
170
171    @classmethod
172    def from_dict(cls, payload: dict[str, object]) -> Term:
173        operands = payload["operands"]
174        if not isinstance(operands, list):
175            raise ValueError("term payload has no operand list")
176        restored: list[Atom | Term] = []
177        for item in operands:  # pyright: ignore[reportUnknownVariableType]
178            if not isinstance(item, dict):
179                raise ValueError("term operand is not an object")
180            if "operation" in item:
181                restored.append(cls.from_dict(item))  # pyright: ignore[reportUnknownArgumentType]
182            else:
183                restored.append(Atom(str(item["feature"]), item["transform"]))  # pyright: ignore[reportArgumentType]
184        return cls(payload["operation"], tuple(restored))  # pyright: ignore[reportArgumentType]
185
186
187def simplify(term: Term) -> Term:
188    """Apply algebraic identities so a term prints in its shortest equivalent form.
189
190    Nesting produces expressions that are correct but redundant: ``([a] / [b]) * [b]``
191    computes ``a`` and says so in six symbols instead of one. Since the whole claim of
192    this project is that the equation can be read, a term that is longer than it needs to
193    be is a defect, not a cosmetic issue.
194
195    Only exact identities are applied -- cancellation of a matching factor against a
196    matching divisor, and division of a term by itself. Nothing is dropped on numerical
197    grounds here; that is ``search.prune``'s job, and keeping the two separate means this
198    function never changes what a term computes.
199    """
200    operands = tuple(simplify(operand) if isinstance(operand, Term) else operand for operand in term.operands)
201    term = Term(term.operation, operands)
202
203    if term.operation == "product":
204        left, right = term.operands
205        # (x / y) * y  ->  x
206        for numerator_side, other in ((left, right), (right, left)):
207            if isinstance(numerator_side, Term) and numerator_side.operation == "ratio":
208                inner_numerator, inner_divisor = numerator_side.operands
209                if inner_divisor == other:
210                    return _as_term(inner_numerator)
211    elif term.operation == "ratio":
212        numerator, divisor = term.operands
213        # x / x  ->  1, which is collinear with the intercept and carries no information
214        if numerator == divisor:
215            return term
216        # (x * y) / y  ->  x
217        if isinstance(numerator, Term) and numerator.operation == "product":
218            first, second = numerator.operands
219            if second == divisor:
220                return _as_term(first)
221            if first == divisor:
222                return _as_term(second)
223    return term
224
225
226def is_trivial(term: Term) -> bool:
227    """Whether a term reduces to a constant, and so duplicates the intercept."""
228    return term.operation == "ratio" and term.operands[0] == term.operands[1]
229
230
231def _as_term(operand: Atom | Term) -> Term:
232    return operand if isinstance(operand, Term) else Term("atom", (operand,))
233
234
235def _guard(denominator: np.ndarray) -> np.ndarray:
236    """Keep a denominator away from zero without changing its sign.
237
238    A backstop only. Terms whose denominator can actually approach zero are never
239    generated (see ``denominator_atom``), so in a library built by ``build_library``
240    this floor is not what any published term relies on.
241    """
242    return np.where(denominator >= 0.0, 1.0, -1.0) * (np.abs(denominator) + DENOMINATOR_FLOOR)
243
244
245def composition_atom(feature: str, columns: dict[str, np.ndarray]) -> Atom:
246    """The operand used when a feature appears as a numerator or inside a sum.
247
248    Strictly positive features enter log-compressed; the rest enter raw. This is the
249    single decision that keeps composite terms on a comparable scale.
250    """
251    values = columns[feature]
252    return Atom(feature, "log" if bool(np.all(values > 0.0)) else "id")
253
254
255def denominator_atom(feature: str, columns: dict[str, np.ndarray]) -> Atom | None:
256    """The operand to divide by, or ``None`` when this feature cannot safely be one.
257
258    Division is the one operation in the vocabulary that can manufacture a term no
259    linear solver can use: as a denominator approaches zero the term's magnitude runs
260    away, one row comes to hold all of its variance, and the fitted weight becomes an
261    artefact of that row rather than a relationship. The library therefore decides
262    eligibility per feature, before any term is built, rather than generating the term
263    and screening it afterwards.
264
265    Two candidate forms are tried in order of preference and the first that qualifies is
266    used:
267
268    * ``log(f)`` -- compresses scale, but is zero at ``f == 1`` and negative below it,
269      so it is only eligible when the feature stays clear of 1;
270    * ``f`` itself -- eligible when the feature stays clear of 0.
271
272    Eligibility is a bound on **dynamic range**, ``max|d| / min|d| <=
273    MAX_DENOMINATOR_RANGE``, not merely ``min|d| > 0``. That is the criterion that
274    matches the failure: a divisor spanning three orders of magnitude makes the ratio
275    span three orders of magnitude too, so one row holds nearly all of the term's
276    variance and its fitted weight describes that row. Requiring only a non-zero
277    minimum would admit exactly those terms.
278
279    Features that qualify under neither form are never used as denominators. In this
280    meta-dataset that rules out ``nr_norm``, ``nr_bin`` and ``nr_outliers`` (which reach
281    0), the log of anything reaching 1, and wide-range counts such as ``nr_inst``
282    (165 to 7.1M) in raw form -- though ``log(nr_inst)`` qualifies comfortably, which is
283    the compression the preference order exists to find.
284    """
285    for candidate in (Atom(feature, "log"), Atom(feature, "id")):
286        if not candidate.is_defined_on(columns):
287            continue
288        with np.errstate(all="ignore"):
289            evaluated = np.abs(candidate.evaluate(columns))
290        if not np.all(np.isfinite(evaluated)):
291            continue
292        smallest = float(np.min(evaluated))
293        if smallest <= DENOMINATOR_FLOOR:
294            continue
295        if float(np.max(evaluated)) / smallest <= MAX_DENOMINATOR_RANGE:
296            return candidate
297    return None
298
299
300def unary_terms(features: tuple[str, ...], columns: dict[str, np.ndarray]) -> list[Term]:
301    """``f``, ``log(f)``, ``sqrt(f)``, ``1/f`` and ``f^2`` for each feature that admits them."""
302    terms: list[Term] = []
303    for feature in features:
304        for transform in TRANSFORMS:
305            atom = Atom(feature, transform)
306            if atom.is_defined_on(columns):
307                terms.append(Term("atom", (atom,)))
308    return terms
309
310
311def pairwise_terms(
312    left: tuple[str, ...],
313    right: tuple[str, ...],
314    columns: dict[str, np.ndarray],
315    *,
316    both_directions: bool,
317) -> list[Term]:
318    """Ratios and products over two feature groups.
319
320    ``both_directions`` controls whether ``a/b`` and ``b/a`` are both emitted. Within a
321    single group the pairs are unordered and one direction suffices; across the dataset
322    and model groups both directions carry different meaning, so both are kept.
323    """
324    pairs = itertools.product(left, right) if both_directions else itertools.combinations(left, 2)
325    terms: list[Term] = []
326    seen: set[str] = set()
327    for first, second in pairs:
328        if first == second:
329            continue
330        a, b = composition_atom(first, columns), composition_atom(second, columns)
331        # Products are always safe; ratios only exist when the divisor is eligible.
332        candidates = [Term("product", (a, b))]
333        denominator_b = denominator_atom(second, columns)
334        if denominator_b is not None:
335            candidates.append(Term("ratio", (a, denominator_b)))
336        if both_directions:
337            denominator_a = denominator_atom(first, columns)
338            if denominator_a is not None:
339                candidates.append(Term("ratio", (b, denominator_a)))
340        for term in candidates:
341            if term.name not in seen:
342                seen.add(term.name)
343                terms.append(term)
344    return terms
345
346
347def sum_ratio_terms(features: tuple[str, ...], columns: dict[str, np.ndarray]) -> list[Term]:
348    """``(f1 + f2) / f3`` over distinct features, with ``f1`` and ``f2`` unordered.
349
350    Pass **every** feature, not one group. An earlier version was called with the dataset
351    features alone, which silently made the highest-arity operation the only one unable to
352    mix a dataset feature with a model one -- precisely the combination that carries E2's
353    entire lift over E1.
354    """
355    terms: list[Term] = []
356    for first, second in itertools.combinations(features, 2):
357        for third in features:
358            if third in (first, second):
359                continue
360            divisor = denominator_atom(third, columns)
361            if divisor is None:
362                continue
363            operands = (
364                composition_atom(first, columns),
365                composition_atom(second, columns),
366                divisor,
367            )
368            terms.append(Term("sum_ratio", operands))
369    return terms
370
371
372def ratio_of_sums_terms(features: tuple[str, ...], columns: dict[str, np.ndarray]) -> list[Term]:
373    """``(f1 + f2) / (f3 + f4)`` -- the four-feature operation.
374
375    Included so that "three features per term is enough" can be **measured** rather than
376    assumed. A sum is a safe divisor far more often than a single feature is, since adding
377    two log-compressed operands moves the result away from zero, so this reaches
378    combinations ``sum_ratio`` cannot.
379    """
380    terms: list[Term] = []
381    for numerator in itertools.combinations(features, 2):
382        rest = [feature for feature in features if feature not in numerator]
383        for divisor in itertools.combinations(rest, 2):
384            operands = (
385                composition_atom(numerator[0], columns),
386                composition_atom(numerator[1], columns),
387                composition_atom(divisor[0], columns),
388                composition_atom(divisor[1], columns),
389            )
390            with np.errstate(all="ignore"):
391                bottom = operands[2].evaluate(columns) + operands[3].evaluate(columns)
392            if not _is_safe_sum(bottom):
393                continue
394            terms.append(Term("ratio_of_sums", operands))
395    return terms
396
397
398def _is_safe_sum(values: np.ndarray) -> bool:
399    """Whether a summed divisor stays clear of zero with a bounded dynamic range."""
400    magnitude = np.abs(values)
401    if not np.all(np.isfinite(magnitude)):
402        return False
403    smallest = float(np.min(magnitude))
404    if smallest <= DENOMINATOR_FLOOR:
405        return False
406    return float(np.max(magnitude)) / smallest <= MAX_DENOMINATOR_RANGE
407
408
409def is_admissible(values: np.ndarray, max_abs_zscore: float = MAX_ABS_ZSCORE) -> bool:
410    """Whether a term is numerically usable as a column of a linear system.
411
412    This is a backstop, not the main line of defence. The terms that used to fail here
413    were unstable by construction -- ratios whose denominator could reach zero -- and
414    those are no longer generated at all (see ``denominator_atom``). What remains is a
415    cheap check that a caller assembling a ``Library`` by hand, or a transform applied
416    to an unexpected feature, cannot smuggle in a column that breaks the solve.
417
418    Three rejections:
419
420    *Not finite, or constant.* A constant column is collinear with the intercept.
421
422    *Effectively constant.* A term whose spread is negligible against its own magnitude
423    passes an absolute variance check but destroys the published equation: fitting
424    happens on standardised terms, so folding the standardisation back into raw units
425    divides the weight by that spread. One near-constant term produced a weight of
426    -1.5e9 against an intercept of +1.5e9 -- algebraically fine, arithmetically fine,
427    and completely unreadable, which defeats the purpose of the exercise.
428
429    *Driven by one row.* A term is rejected when a single row sits more than
430    ``max_abs_zscore`` standard deviations from the mean, because its fitted weight
431    would describe that row rather than a relationship.
432
433    The cap is sample-size dependent and has to be chosen with that in mind: a single
434    outlier among ``n`` rows can reach a z-score of at most about ``sqrt(n)``, so a cap
435    of 8 constrains the 484-row fit but would be inert on the 20-row aggregated fit,
436    which is why the two configurations do not share a value.
437    """
438    if not np.all(np.isfinite(values)):
439        return False
440    spread = float(np.std(values))
441    if spread <= 1e-12:
442        return False
443    if spread < MIN_RELATIVE_SPREAD * (1.0 + abs(float(np.mean(values)))):
444        return False
445    return float(np.max(np.abs(values - np.mean(values)))) / spread <= max_abs_zscore
446
447
448def _unit(values: np.ndarray) -> np.ndarray:
449    """The centred column scaled to unit length -- the form collinearity is judged in.
450
451    Centring first is what makes the test match the fit: fitting happens on standardised
452    terms, so two columns differing by an additive or multiplicative constant are one
453    column as far as the solve is concerned, however different their raw values look.
454    """
455    centred = values - np.mean(values)
456    norm = float(np.linalg.norm(centred))
457    return centred / norm if norm > 0.0 else centred
458
459
460class Library:
461    """A set of terms together with the design matrix they produce.
462
463    Terms are dropped for three reasons, in order: a repeated *name*, failing
464    `is_admissible`, and -- the subtlest -- being the same column as a term already kept.
465
466    That last check cannot be done on names. ``inst_to_attr`` is ``nr_inst / nr_attr`` in
467    this meta-dataset, so ``[log(inst_to_attr)] + [log(nr_attr)]`` **is**
468    ``log(nr_inst)``, exactly, and the grammar generates both. Three such pairs exist in
469    the published 206-term library. A binary feature produces them too and more bluntly:
470    ``log``, ``sqrt`` and ``1/f`` are all undefined at zero, leaving only ``f`` and
471    ``f^2``, which for a 0/1 column are the same numbers under two names.
472
473    Keeping both members of a pair is wasteful rather than dangerous, and the distinction
474    is worth being precise about. The beam search already refuses a candidate whose
475    correlation with a selected term exceeds ``ml_meta_perf.search.COLLINEARITY_LIMIT`` (0.95), so
476    a duplicate pair cannot both be selected on that path and no singular system arises
477    there. What the duplicates cost is candidate-pool slots, search time, and a place in
478    the reported term rankings, where they appear as two independent findings. The check
479    here also covers the paths ``_blocked`` does not: a ``Library`` assembled by hand, and
480    any consumer reading `terms` or `matrix` directly.
481
482    The first term of a pair wins. Generation order runs simple to complex, so the survivor
483    is the shorter form: ``[log(nr_inst)] / [log(Processing Units Number)]`` is kept and
484    ``([log(inst_to_attr)] + [log(nr_attr)]) / [log(Processing Units Number)]`` is dropped.
485    """
486
487    def __init__(
488        self,
489        terms: list[Term],
490        columns: dict[str, np.ndarray],
491        *,
492        max_abs_zscore: float = MAX_ABS_ZSCORE,
493    ) -> None:
494        kept: list[Term] = []
495        vectors: list[np.ndarray] = []
496        seen: set[str] = set()
497        # The survivors' unit vectors as one growing array rather than a list, so the
498        # duplicate test is a single matrix-vector product against everything kept so far
499        # instead of a Python loop of dot products over it. Identical arithmetic and
500        # identical order -- the first term of a pair still wins -- but the loop form ran
501        # 590k times per study and `search.guided_screen` already does it this way.
502        accepted: np.ndarray | None = None
503        for term in terms:
504            if term.name in seen:
505                continue
506            with np.errstate(all="ignore"):
507                values = term.evaluate(columns)
508            if not is_admissible(values, max_abs_zscore):
509                continue
510            unit = _unit(values)
511            if accepted is None:
512                accepted = np.empty((len(terms), values.shape[0]))
513            elif float(np.abs(accepted[: len(kept)] @ unit).max()) > COLLINEARITY_TOLERANCE:
514                continue
515            accepted[len(kept)] = unit
516            seen.add(term.name)
517            kept.append(term)
518            vectors.append(values)
519        if not kept:
520            raise ValueError("term library is empty after filtering")
521        self.terms: list[Term] = kept
522        self.matrix: np.ndarray = np.column_stack(vectors)
523
524    @property
525    def names(self) -> list[str]:
526        return [term.name for term in self.terms]
527
528    @property
529    def feature_groups(self) -> np.ndarray:
530        """Which terms describe the *same combination of raw features*, as group ids.
531
532        One id per term, ``-1`` for a term over a single feature. Two terms share an id
533        when their feature multisets reduce to the same set, however differently they
534        arrange it: ``[log(a)] / [log(b)]`` and ``[log(b)] / [log(a)]`` are one group, and
535        so are ``[a] * [log(b)]`` and ``[a] / [log(b)]``.
536
537        `ml_meta_perf.search.Selector` refuses to place two terms of one group in the same
538        equation. That is a **readability** rule, not a numerical one, and the distinction
539        matters because the numerical guard already passes: the two mirrored pairs this
540        removed from the previous 16-term equation correlated at 0.891 and 0.786, both
541        under ``search.COLLINEARITY_LIMIT``, in a design conditioned at 7.8. Nothing was
542        ill-posed. What was wrong is that the equation spent two of its sixteen slots
543        writing one relationship both ways up -- ``log(PUN)/log(nr_class)`` beside
544        ``log(nr_class)/log(PUN)``, *both* carrying negative weight -- and the term table
545        then reported them as two independent findings, each reading "lowers MCC". Jointly
546        they encode a curvature in one ratio; separately neither sentence is true.
547
548        An equation this study cannot reason about term by term has failed its purpose,
549        so the constraint is part of the protocol rather than an option. It costs nothing
550        measurable: paired over the twenty leave-one-dataset-out (LODO) folds it moves mean
551        absolute error by +0.0018 at twenty terms, worse on ten of twenty datasets, sign
552        test p = 1.000 and a bootstrap interval of [-0.0039, +0.0084] that spans zero.
553
554        A term over one feature is left ungrouped, so ``1/f`` may still sit beside
555        ``f^2``. Those are two points of one curve over one feature, which reads as a
556        shape; the constraint is about a *relationship* being stated twice.
557        """
558        identifiers: dict[frozenset[str], int] = {}
559        groups = np.full(len(self.terms), -1, dtype=np.int64)
560        for index, term in enumerate(self.terms):
561            features = frozenset(term.features)
562            if len(features) < 2:
563                continue
564            groups[index] = identifiers.setdefault(features, len(identifiers))
565        return groups
566
567    def __len__(self) -> int:
568        return len(self.terms)
569
570    def design(self, columns: dict[str, np.ndarray], indices: list[int]) -> np.ndarray:
571        """Re-evaluate a subset of terms on fresh columns, in the given order."""
572        with np.errstate(all="ignore"):
573            return np.column_stack([self.terms[index].evaluate(columns) for index in indices])
574
575
576def build_library(
577    dataset_features: tuple[str, ...],
578    model_features: tuple[str, ...],
579    columns: dict[str, np.ndarray],
580    *,
581    include_sum_ratio: bool = True,
582    max_arity: int = 3,
583    max_abs_zscore: float = MAX_ABS_ZSCORE,
584) -> Library:
585    """Assemble the full candidate vocabulary.
586
587    With ``model_features`` empty this yields the E1 vocabulary: unary terms over the
588    dataset features plus their internal ratios, products and sums. With model features
589    supplied it additionally yields the cross terms, which is where E2 gets its lift --
590    dataset features alone cannot express anything that varies within a dataset.
591
592    ``max_arity`` caps how many distinct raw features one term may combine: 2 for atoms,
593    ratios and products only, 3 to add ``(f1+f2)/f3``, 4 to add ``(f1+f2)/(f3+f4)``. It is
594    a parameter rather than a constant because "three is enough" is a claim that has to be
595    measured, and measuring it means being able to build the alternatives.
596    """
597    # Sorted, so the library is a function of the feature *sets* and not of the order the
598    # caller happened to pass them in. It was not, and the consequence was not cosmetic:
599    # `pairwise_terms` names a product after whichever operand it sees first, so `A * B` and
600    # `B * A` -- the same column, multiplication being commutative -- entered under two names
601    # depending on declaration order. `Library` then de-duplicated by *correlation with terms
602    # already kept*, so which of a near-collinear pair survived also depended on order, and
603    # across five orderings of one six-feature set the library came out at 270, 270, 272, 270
604    # and 271 terms with 24 of 257 pool slots differing. That is the whole of the
605    # LODO ordering band recorded in the branch notes;
606    # it was never the beam's tie-breaking.
607    dataset_features = tuple(sorted(dataset_features))
608    model_features = tuple(sorted(model_features))
609    features = dataset_features + model_features
610    terms = unary_terms(features, columns)
611    terms += pairwise_terms(dataset_features, dataset_features, columns, both_directions=False)
612    if model_features:
613        terms += pairwise_terms(model_features, model_features, columns, both_directions=False)
614        terms += pairwise_terms(dataset_features, model_features, columns, both_directions=True)
615    if include_sum_ratio and max_arity >= 3:
616        # Every feature, not just the dataset ones: the highest-arity operation must be
617        # able to mix the two groups like every other operation can.
618        terms += sum_ratio_terms(features, columns)
619    if max_arity >= 4:
620        terms += ratio_of_sums_terms(features, columns)
621    return Library(terms, columns, max_abs_zscore=max_abs_zscore)