Statistics

The stats stage: assembling the analysis table, testing group differences, correlating features with outcomes, and writing the plain-English report. Principal component analysis (with varimax rotation) lives here too, as does the scorer that applies one or more saved models to a new dataset, and the figures the reports draw (word clouds, and the charts training reports use).

taters.stats.assemble

Build the analysis table: every chosen feature table, one row per text, joined with the metadata columns the statistics need.

The feature extractors each write their own CSV keyed by text_id (or by a multi-column key for media runs). Statistics want ONE wide table -- features side by side, with the grouping and outcome columns from the user's spreadsheet -- and nothing in Taters joined columns before this module: the gathers stack rows, they never widen.

The join is INNER, with accounting. A left join would seed missing values into every downstream analysis for rows some extractor skipped (a zero-token document missing from a term matrix, say); keeping only rows every table has is the honest denominator, and honesty about the cost lives in the manifest: rows in, rows surviving each join, rows each filter removed. Those counts are also written into the report, because "we analyzed 99,881 of your 100,000 rows, and here is where the rest went" is the first thing a reviewer asks.

Filters are the "select * where" of the stage: [column, op, value] triples a row must ALL satisfy to stay. A blank (missing) cell fails every filter -- a row whose word count is unknown does not sneak past word_count >= 25.

assemble_analysis_table

assemble_analysis_table(
    *,
    feature_csvs,
    filter_csvs=(),
    metadata_csv=None,
    metadata_cols=(),
    key_cols=("text_id",),
    text_cols=(),
    split_col=None,
    filters=None,
    bookkeeping="aside",
    bookkeeping_cols=(),
    out_dir="stats_results",
    out_csv=None,
    keep_table=True,
    overwrite_existing=False,
    on_progress=None,
    verbose=True,
    encoding="utf-8-sig"
)

Join feature tables and metadata into one wide analysis table.

Parameters:

Name	Type	Description	Default
feature_csvs	Sequence[PathLike]	The feature tables to include, each with exactly one row per key. Non-numeric columns other than the key are dropped (with a note): the analyses consume numbers, and a stray text column riding along in a feature file is weight, not information.	required
text_cols	Sequence[str]	Columns to carry through even though they are not numbers. Needed because a saved model fitted with a categorical control knows a predictor called gender=male and the table being scored still holds gender as the words the researcher typed -- dropped as "non-numeric", the control could never be rebuilt. Named explicitly, never guessed: a stray text column riding along in a feature file really is weight rather than information.	()
filter_csvs	Sequence[PathLike]	Tables joined in so their columns can be filtered on, but never offered to the analyses as features. A word count computed because someone asked to drop texts under 25 words is a gate, not a predictor -- and it turned up as one in a ridge regression, which is the kind of result that looks fine and is not (a real report). Their columns land in the analysis table and in the sidecar's filters list; only feature_csvs columns reach sets, which is what every analysis reads to decide what a feature is.	()
metadata_csv	Optional[PathLike]	Optional table holding the grouping/outcome columns (typically the run's gathered/metadata.csv). When omitted, the table is features only -- fine for a PCA, useless for an ANOVA.	None
metadata_cols	Sequence[str]	Which metadata columns to carry. Empty means all of them (minus any column literally named text).	()
key_cols	Sequence[str]	The join key. ("text_id",) for text runs; media runs can key on several columns.	('text_id',)
split_col	Optional[str]	A second dimension the feature tables are keyed on but the metadata is not -- "source_col", when a spreadsheet's text columns were measured one at a time. Each text then has one feature row per column, while the participant has one row of metadata, so the features join to each other on the key plus this column and the metadata broadcasts across them. The analyses split on it in turn, so one person's several answers are never treated as several independent observations.	None
filters	Optional[Sequence[Sequence]]	Row filters, [column, op, value] triples that a row must ALL satisfy to stay. Operators: == != < <= > >= in not_in. A blank cell fails every filter. Each filter's removed-row count is recorded.	None
bookkeeping	Literal['aside', 'features']	What to do with the count columns the steps write beside their measures -- a matrix's token_count, a dictionary's WC, a readability step's raw sentence and syllable counts -- which each analyzer declares in the record beside its table. "aside" (the default) joins them into the table, where they can be filtered on and looked at, but keeps them out of the feature sets every analysis reads; "features" analyses them like any other measure. A ridge once learned from token_count because nothing told it apart.	'aside'
bookkeeping_cols	Sequence[str]	Column names to keep aside as well, for a table with no record to declare them (one made by hand, or by an older build).	()
out_dir	PathLike	Where the table goes; out_csv defaults to <out_dir>/analysis_table.csv. Two sidecars land next to it: analysis_table_sets.json (which columns came from which feature table -- the handle for per-feature-set analyses) and assemble_manifest.json (the row accounting).	'stats_results'
out_csv	PathLike	Where the table goes; out_csv defaults to <out_dir>/analysis_table.csv. Two sidecars land next to it: analysis_table_sets.json (which columns came from which feature table -- the handle for per-feature-set analyses) and assemble_manifest.json (the row accounting).	'stats_results'
keep_table	bool	Whether the merged table survives the run. It is written either way -- every analysis reads it from disk -- but with False the report step deletes it at the end. Keep it when you want the dataset the statistics were computed on, which is usually: it is the file to hand to a colleague, or to open in R. Turn it off when it is enormous and derivable, since a wide feature set makes it the largest thing the run produces.	True
overwrite_existing	bool	When False (default) and the table exists, it is returned untouched.	False

Returns:

Type	Description
Path	The analysis table.

Source code in src\taters\stats\assemble.py

def assemble_analysis_table(
    *,
    feature_csvs: Sequence[PathLike],
    filter_csvs: Sequence[PathLike] = (),
    metadata_csv: Optional[PathLike] = None,
    metadata_cols: Sequence[str] = (),
    key_cols: Sequence[str] = ("text_id",),
    text_cols: Sequence[str] = (),
    split_col: Optional[str] = None,
    filters: Optional[Sequence[Sequence]] = None,
    bookkeeping: Literal["aside", "features"] = "aside",
    bookkeeping_cols: Sequence[str] = (),
    out_dir: PathLike = "stats_results",
    out_csv: Optional[PathLike] = None,
    keep_table: bool = True,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[..., None]] = None,
    verbose: bool = True,
    encoding: str = "utf-8-sig",
) -> Path:
    """
    Join feature tables and metadata into one wide analysis table.

    Parameters
    ----------
    feature_csvs
        The feature tables to include, each with exactly one row per key.
        Non-numeric columns other than the key are dropped (with a note):
        the analyses consume numbers, and a stray ``text`` column riding
        along in a feature file is weight, not information.
    text_cols
        Columns to carry through even though they are not numbers. Needed
        because a saved model fitted with a categorical control knows a
        predictor called ``gender=male`` and the table being scored still
        holds ``gender`` as the words the researcher typed -- dropped as
        "non-numeric", the control could never be rebuilt. Named
        explicitly, never guessed: a stray ``text`` column riding along in
        a feature file really is weight rather than information.
    filter_csvs
        Tables joined in so their columns can be *filtered on*, but never
        offered to the analyses as features. A word count computed because
        someone asked to drop texts under 25 words is a gate, not a
        predictor -- and it turned up as one in a ridge regression, which is
        the kind of result that looks fine and is not (a real report). Their
        columns land in the analysis table and in the sidecar's ``filters``
        list; only ``feature_csvs`` columns reach ``sets``, which is what
        every analysis reads to decide what a feature is.
    metadata_csv
        Optional table holding the grouping/outcome columns (typically the
        run's ``gathered/metadata.csv``). When omitted, the table is features
        only -- fine for a PCA, useless for an ANOVA.
    metadata_cols
        Which metadata columns to carry. Empty means all of them (minus any
        column literally named ``text``).
    key_cols
        The join key. ``("text_id",)`` for text runs; media runs can key on
        several columns.
    split_col
        A second dimension the *feature* tables are keyed on but the
        metadata is not -- ``"source_col"``, when a spreadsheet's text
        columns were measured one at a time. Each text then has one feature
        row per column, while the participant has one row of metadata, so
        the features join to each other on the key plus this column and the
        metadata broadcasts across them. The analyses split on it in turn,
        so one person's several answers are never treated as several
        independent observations.
    filters
        Row filters, ``[column, op, value]`` triples that a row must ALL
        satisfy to stay. Operators: ``== != < <= > >= in not_in``. A blank
        cell fails every filter. Each filter's removed-row count is recorded.
    bookkeeping
        What to do with the count columns the steps write beside their
        measures -- a matrix's ``token_count``, a dictionary's ``WC``, a
        readability step's raw sentence and syllable counts -- which each
        analyzer declares in the record beside its table. ``"aside"`` (the
        default) joins them into the table, where they can be filtered on
        and looked at, but keeps them out of the feature sets every analysis
        reads; ``"features"`` analyses them like any other measure. A ridge
        once learned from ``token_count`` because nothing told it apart.
    bookkeeping_cols
        Column names to keep aside as well, for a table with no record to
        declare them (one made by hand, or by an older build).
    out_dir, out_csv
        Where the table goes; ``out_csv`` defaults to
        ``<out_dir>/analysis_table.csv``. Two sidecars land next to it:
        ``analysis_table_sets.json`` (which columns came from which feature
        table -- the handle for per-feature-set analyses) and
        ``assemble_manifest.json`` (the row accounting).
    keep_table
        Whether the merged table survives the run. It is written either way
        -- every analysis reads it from disk -- but with False the report
        step deletes it at the end. Keep it when you want the dataset the
        statistics were computed on, which is usually: it is the file to
        hand to a colleague, or to open in R. Turn it off when it is
        enormous and derivable, since a wide feature set makes it the
        largest thing the run produces.
    overwrite_existing
        When False (default) and the table exists, it is returned untouched.

    Returns
    -------
    Path
        The analysis table.
    """
    import json


    feature_paths = [Path(p) for p in (feature_csvs or [])]
    filter_paths = [Path(p) for p in (filter_csvs or [])]
    if not feature_paths:
        raise ValueError(
            "feature_csvs is empty: nothing to assemble"
            + (" -- every table given was a filter-only one, and an analysis "
               "needs something to analyze" if filter_paths else ""))
    key_cols = [str(k) for k in key_cols]
    if not key_cols:
        raise ValueError("key_cols is empty: there is nothing to join on")
    split_col = str(split_col) if split_col else None
    # the feature tables carry the extra (split) dimension; the metadata doesn't.
    feature_key = key_cols + ([split_col] if split_col else [])

    out_dir = Path(out_dir)
    out_path = Path(out_csv) if out_csv else out_dir / "analysis_table.csv"
    out_path.parent.mkdir(parents=True, exist_ok=True)
    sets_path = out_path.with_name(out_path.stem + "_sets.json")
    manifest_path = out_path.parent / "assemble_manifest.json"

    if reusable(out_path, feature_csvs, filter_csvs, metadata_csv,
                overwrite_existing=overwrite_existing, verbose=verbose,
                what="the analysis table"):
        if verbose:
            print(f"Analysis table already exists; returning existing file: "
                  f"{out_path}")
        return out_path

    if bookkeeping not in ("aside", "features"):
        raise ValueError(
            f"bookkeeping must be 'aside' or 'features', got {bookkeeping!r}")
    manifest: dict = {"key": feature_key, "joins": [], "filters": [],
                      "filter_only": [], "bookkeeping": {},
                      "bookkeeping_mode": bookkeeping,
                      "dropped_columns": {}, "renamed_columns": {},
                      # we write this down rather than just acting on it,
                      # so that a folder with no analysis_table.csv in it
                      # can say whether that was a choice or a failure.
                      "kept": bool(keep_table)}

    # ------------------------------------------------------------------ read
    announce(on_progress, "reading the feature tables")
    frames: List[tuple] = []      # (stem, DataFrame w/ key + kept cols, is_feature)
    seen_stems: dict = {}
    # per table, the numeric columns that aren't actually measures: whatever
    # the analyzer declared in the record beside the table (see
    # `records_settings(bookkeeping=...)`), plus anything the caller names.
    # we join these in but keep them out of the feature sets unless asked --
    # we once had a ridge on a document-term matrix happily learning from
    # `token_count`. oops.
    aside: dict = {}
    named = {str(c) for c in bookkeeping_cols}
    for path, is_feature in ([(p, True) for p in feature_paths]
                             + [(p, False) for p in filter_paths]):
        if not path.is_file():
            raise FileNotFoundError(f"feature table not found: {path}")
        stem = path.stem
        if stem in seen_stems:
            raise ValueError(
                f"two feature tables share the name {stem!r} "
                f"({seen_stems[stem]} and {path}); rename one -- the name "
                f"labels its feature set in every result table.")
        seen_stems[stem] = path
        df = read_str_csv(path, encoding=encoding)
        missing = [k for k in feature_key if k not in df.columns]
        if missing:
            raise ValueError(
                f"{path.name} is missing key column(s) {missing}; it cannot "
                f"be joined into the analysis table.")
        _refuse_duplicate_keys(df, feature_key, path.name)
        kept, dropped = [], []
        wanted = {str(c) for c in text_cols}
        for col in df.columns:
            if col in feature_key:
                continue
            if col in wanted or looks_numeric(df[col].tolist()):
                kept.append(col)
            else:
                dropped.append(col)
        if dropped:
            manifest["dropped_columns"][stem] = dropped
            if verbose:
                print(f"[assemble] {path.name}: dropping non-numeric "
                      f"column(s) {dropped}")
        if not kept:
            raise ValueError(
                f"{path.name} has no numeric feature columns beyond the key; "
                f"there is nothing in it to analyze.")
        frames.append((stem, df[feature_key + kept], is_feature))
        if not is_feature:
            manifest["filter_only"].append(stem)
        elif bookkeeping == "aside":
            declared = set(_declared_bookkeeping(path)) | named
            found = [c for c in kept if c in declared]
            if found:
                aside[stem] = found
                manifest["bookkeeping"][stem] = found

    # ----------------------------------------------------------- metadata
    meta = None
    meta_kept: List[str] = []
    if metadata_csv is not None:
        announce(on_progress, "reading the metadata")
        meta = read_str_csv(Path(metadata_csv), encoding=encoding)
        missing = [k for k in key_cols if k not in meta.columns]
        if missing:
            raise ValueError(
                f"{Path(metadata_csv).name} is missing key column(s) "
                f"{missing}; it cannot anchor the analysis table.")
        _refuse_duplicate_keys(meta, key_cols, Path(metadata_csv).name)
        if metadata_cols:
            absent = [c for c in metadata_cols if c not in meta.columns]
            if absent:
                spare = [c for c in meta.columns if c not in key_cols]
                raise ValueError(
                    f"metadata_cols not in {Path(metadata_csv).name}: "
                    f"{name_a_few(absent)}. It has {len(spare)}: "
                    f"{name_a_few(spare)}.")
            meta_kept = [str(c) for c in metadata_cols]
        else:
            meta_kept = [c for c in meta.columns
                         if c not in key_cols and c != "text"]
        meta = meta[key_cols + meta_kept]

    # ---------------------------------------------------------- collisions
    # here, we deal with name clashes. any feature column used by two feature
    # tables (or clashing with a key/metadata name) gets renamed
    # "<stem>.<col>" in EVERY feature table that has it -- all the colliders,
    # not just the later ones, so that the outcome doesn't depend on file
    # order. keys and metadata keep their names.
    reserved = set(feature_key) | set(meta_kept)
    counts: dict = {}
    for _, df, _is_feature in frames:
        for col in df.columns:
            if col not in feature_key:
                counts[col] = counts.get(col, 0) + 1
    colliding = {c for c, n in counts.items() if n > 1} | \
                {c for c in counts if c in reserved}
    sets: dict = {}
    filter_columns: dict = {}
    bookkeeping_columns: dict = {}
    # we also keep track of where each table came from. without this, a model
    # knows the *names* of its predictors but has no way of getting back to
    # the record of how they were measured -- and that's the whole point of
    # recording it.
    sources: dict = {}
    renamed_frames = []
    for stem, df, is_feature in frames:
        renames = {c: f"{stem}.{c}" for c in df.columns
                   if c not in feature_key and c in colliding}
        if renames:
            df = df.rename(columns=renames)
            manifest["renamed_columns"].update(renames)
        columns_here = [c for c in df.columns if c not in feature_key]
        # only a feature table's columns become a *set*. every analysis reads
        # the sidecar's `sets` to decide what counts as a feature, so if a
        # filter-only table landed there we'd end up predicting from it. a
        # table's bookkeeping columns (renamed along with the rest if they
        # collided) go beside the sets, not in them.
        kept_aside = {renames.get(c, c) for c in aside.get(stem, ())}
        if is_feature and kept_aside:
            bookkeeping_columns[stem] = [c for c in columns_here
                                         if c in kept_aside]
            columns_here = [c for c in columns_here if c not in kept_aside]
        (sets if is_feature else filter_columns)[stem] = columns_here
        if is_feature and stem in seen_stems:
            sources[stem] = str(Path(seen_stems[stem]).resolve())
        renamed_frames.append((stem, df))

    # --------------------------------------------------------------- join
    announce(on_progress, "joining the tables")
    if meta is not None:
        table = meta
        start_label = Path(metadata_csv).name
    else:
        stem, table = renamed_frames[0]
        start_label = stem
        renamed_frames = renamed_frames[1:]
    manifest["rows_start"] = int(len(table))
    manifest["start"] = start_label
    for stem, df in renamed_frames:
        before = int(len(table))
        # the metadata is one row per text; a feature table under a split is
        # several. joining on what they share broadcasts the metadata across
        # all of a text's rows -- that's what "this participant's answers"
        # means, so it's what we want.
        on = [c for c in feature_key if c in table.columns and c in df.columns]
        # rows of this table with no partner so far get lost too, and the
        # before/after count can't see them: a text that's in every table but
        # the first would vanish from the join without showing up in any
        # number. what it's really missing from is an *earlier* table, so
        # that's how we count it.
        unmatched = int((df[on].drop_duplicates()
                         .merge(table[on].drop_duplicates(), on=on,
                                how="left", indicator=True)["_merge"]
                         == "left_only").sum())
        table = table.merge(df, how="inner", on=on)
        manifest["joins"].append({"table": stem, "rows_before": before,
                                  "rows_after": int(len(table)),
                                  "rows_unmatched": unmatched})
        if on_progress is not None:
            on_progress(len(manifest["joins"]), len(frames),
                        f"joined {stem}")

    if len(table) == 0:
        raise ValueError(
            "the join produced no rows: no key value appears in every input. "
            "Check that the metadata and feature tables describe the same "
            "texts and share the same key columns "
            f"({key_cols}). Row counts: {manifest['joins']}")

    # ------------------------------------------------------------- filters
    checked = _validate_filters(filters, table.columns, table)
    for column, op, value in checked:
        before = int(len(table))
        mask = table[column].map(lambda cell: _matches(cell, op, value))
        table = table[mask]
        manifest["filters"].append(
            {"filter": [column, op, value], "removed": before - int(len(table))})
    if checked and len(table) == 0:
        raise ValueError(
            "the filters removed every row. Removed per filter: "
            + "; ".join(f"{f['filter']} removed {f['removed']}"
                        for f in manifest["filters"]))
    manifest["rows_final"] = int(len(table))

    # -------------------------------------------------------------- write
    announce(on_progress, "writing the analysis table")
    with atomic_write(out_path, mode="w", encoding=encoding,
                      newline="") as fh:
        table.to_csv(fh, index=False)
    with atomic_write(sets_path, mode="w", encoding="utf-8") as fh:
        json.dump({"key": feature_key, "split": split_col or "",
                   "metadata": meta_kept, "sets": sets,
                   "sources": sources,
                   "filters": filter_columns,
                   "bookkeeping": bookkeeping_columns}, fh, indent=1)
    manifest["taters"] = taters_version()
    with atomic_write(manifest_path, mode="w", encoding="utf-8") as fh:
        json.dump(manifest, fh, indent=1)

    write_section(out_path.parent, "assemble",
                  _section_md(manifest, sets))
    if verbose:
        print(f"[assemble] {manifest['rows_final']:,} rows × "
              f"{len(table.columns):,} columns -> {out_path}")
    return out_path

taters.stats.group_differences

Group differences: one-way ANOVA per feature, with post-hoc tests that agree with the omnibus test's assumptions.

The default is the classic decomposition -- pooled-variance F with Tukey-Kramer post-hocs (Tukey's HSD generalized to unequal group sizes). welch=True switches BOTH tests: Welch's F (variances not pooled, Satterthwaite degrees of freedom) with Games-Howell post-hocs. The pairing is deliberate: Tukey's q statistic divides by a pooled error term, which is incoherent if you chose Welch precisely because you distrust pooling; Games-Howell is the standard heteroscedastic companion, built from per-pair Welch degrees of freedom and the same studentized-range distribution.

Multiple comparisons are handled twice, at two levels: post-hoc p-values are already family-adjusted within a feature (that is what the studentized range does), and the omnibus p-values are corrected across features within each feature set by the method you choose (p_adjust; Benjamini-Hochberg by default) -- testing 160 cohesion measures at alpha=.05 without correction would hand you eight "findings" by chance alone.

Missing data: listwise per feature (a row with a blank feature value or a blank group cell sits out that feature's test); the per-group n_* columns make the cost visible. NA is never zero.

analyze_group_differences

analyze_group_differences(
    *,
    table_csv,
    group_col,
    feature_sets=None,
    control_cols=(),
    categorical_controls=(),
    pca="off",
    pca_components=0,
    pca_retain="parallel",
    pca_rotation=True,
    pca_max_missing=MAX_MISSING,
    split_col=None,
    welch=False,
    posthoc="auto",
    p_adjust="fdr_bh",
    alpha=0.05,
    out_dir=None,
    overwrite_existing=False,
    on_progress=None,
    verbose=True,
    encoding="utf-8-sig",
    rounding=4
)

Test every feature for differences between the groups in one column.

Parameters:

Name	Type	Description	Default
table_csv	PathLike	The assembled analysis table (see the assemble step).	required
group_col	str	The column naming each row's group. Rows with a blank group sit out.	required
feature_sets		Which features to test: None for all of them as one family, "per_table" to repeat per source feature table (using the sidecar the assemble step wrote), or {name: [columns]}.	None
control_cols	Sequence[str]	Columns to hold constant, which turns each comparison into an analysis of covariance: the groups are compared on what is left once these have had their say, and the table reports each group's adjusted mean -- its predicted value at the average age, the average word count -- rather than its raw one. Continuous controls enter as themselves; categorical ones become indicators against a reference level the report names. eta2 then holds partial eta-squared, the share of what the controls left over, and method says ancova so the two are never confused.	()
categorical_controls	Sequence[str]	Which of control_cols to treat as categories despite looking numeric -- a site coded 1/2/3.	()
pca	str or list of str	Analyze components instead of the raw measures. "off" uses the features as they are; "all" reduces every feature set; a list of set names reduces those and leaves the rest alone -- which is the case worth having, since a hundred dictionary categories are worth reducing and eight readability indices are not. Set per analysis, deliberately: raw variables read better in a correlation table, where each row is a measure you can name, while a ridge over four hundred collinear measures is what components are for. The loadings land beside these results and named after them, because a component means nothing without the table saying what loads on it.	"off"
pca_components	int	How many components to keep, when reducing. 0 decides by the Kaiser criterion -- a starting point, not an answer.	0
pca_retain	('parallel', 'kaiser')	How the component count is chosen when pca_components is 0: parallel analysis keeps a component while its eigenvalue beats what random data of the same size produce at that rank; the Kaiser rule keeps every eigenvalue above 1, which on a wide table is most of them.	"parallel"
pca_rotation	bool	Rotate the components (varimax) so each loads on a small cluster of features and is therefore nameable.	True
pca_max_missing	float	When reducing, set aside any feature missing for more than this fraction of the rows rather than letting it delete them -- the same rule, and the same default, the prediction steps apply to predictors. Of the rows that remain, the components are fitted on those with every kept feature present and a row still missing one is left unscored rather than guessed at; both counts are reported.	0.5
split_col	Optional[str]	Run the whole analysis once per value of this column, labeling each result with it. Set to "source_col" when a spreadsheet's text columns were measured separately: a participant then has one row per column, and analyzing those together would count one person's several answers as several independent observations -- inflating the degrees of freedom, and every p-value with them.	None
welch	bool	False (default): classic pooled ANOVA. True: Welch's F, for when the groups' variances should not be pooled.	False
posthoc	('auto', 'tukey', 'games_howell', 'bonferroni', 'none')	Which pairwise comparison to run, and therefore how the pairwise p-values are corrected: "auto" (default) -- Tukey-Kramer under a pooled ANOVA, Games-Howell under Welch. The pairing matters: Tukey's statistic divides by a pooled error term, which is incoherent once you have chosen Welch precisely because you distrust pooling. "tukey" / "games_howell" -- force one of them. "bonferroni" -- pairwise t-tests, Bonferroni-corrected. "none" -- pairwise t-tests, uncorrected. Six pairs at .05 give a ~26% chance of at least one false positive, so this is a deliberate choice, not a shortcut.	"auto"
p_adjust	('none', 'fdr_bh', 'fdr_by', 'holm', 'bonferroni')	How the omnibus p-values are adjusted across features, which is a separate family from the pairwise one. "fdr_by" is valid under any dependence between features, and stricter for it; "holm" and "bonferroni" control the family-wise error rate instead. With "none" no p_adj column is written -- an unadjusted number under an "adjusted" heading is worse than no column.	"none"
alpha	float	The threshold used for the pairwise_sig summary and the post-hoc confidence intervals. It does not gate what is written -- every test is in the tables.	0.05
out_dir	Optional[PathLike]	Where the two CSVs go; defaults to the analysis table's folder.	None
overwrite_existing	bool	When False (default) and the results already exist, they are returned untouched instead of recomputed.	False
rounding	int	Decimal places in the output tables. Values too small for it (a p of 3e-12) switch to significant digits rather than collapsing to 0.	4

Returns:

Type	Description
Path	group_differences.csv (one row per feature; the pairwise detail is in group_differences_pairwise.csv beside it).

Notes

A p-value printed as 0 means the tail underflowed double precision (roughly p < 1e-300), not that the probability is zero. The columns stay numeric so R, Excel and pandas read them as numbers; report such a value the way every statistics package does, as "p < .001".

Source code in src\taters\stats\group_differences.py

def analyze_group_differences(
    *,
    table_csv: PathLike,
    group_col: str,
    feature_sets=None,
    control_cols: Sequence[str] = (),
    categorical_controls: Sequence[str] = (),
    pca: object = "off",
    pca_components: int = 0,
    pca_retain: Literal["parallel", "kaiser"] = "parallel",
    pca_rotation: bool = True,
    pca_max_missing: float = MAX_MISSING,
    split_col: Optional[str] = None,
    welch: bool = False,
    posthoc: str = "auto",
    p_adjust: str = "fdr_bh",
    alpha: float = 0.05,
    out_dir: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[..., None]] = None,
    verbose: bool = True,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
) -> Path:
    """
    Test every feature for differences between the groups in one column.

    Parameters
    ----------
    table_csv
        The assembled analysis table (see the assemble step).
    group_col
        The column naming each row's group. Rows with a blank group sit out.
    feature_sets
        Which features to test: ``None`` for all of them as one family,
        ``"per_table"`` to repeat per source feature table (using the
        sidecar the assemble step wrote), or ``{name: [columns]}``.
    control_cols
        Columns to hold constant, which turns each comparison into an
        **analysis of covariance**: the groups are compared on what is left
        once these have had their say, and the table reports each group's
        adjusted mean -- its predicted value at the average age, the average
        word count -- rather than its raw one. Continuous controls enter as
        themselves; categorical ones become indicators against a reference
        level the report names. ``eta2`` then holds *partial* eta-squared,
        the share of what the controls left over, and ``method`` says
        ``ancova`` so the two are never confused.
    categorical_controls
        Which of ``control_cols`` to treat as categories despite looking
        numeric -- a site coded 1/2/3.
    pca : str or list of str, default="off"
        Analyze components instead of the raw measures. ``"off"`` uses the
        features as they are; ``"all"`` reduces every feature set; a list of
        set names reduces those and leaves the rest alone -- which is the
        case worth having, since a hundred dictionary categories are worth
        reducing and eight readability indices are not.

        Set per analysis, deliberately: raw variables read better in a
        correlation table, where each row is a measure you can name, while a
        ridge over four hundred collinear measures is what components are
        for. The loadings land beside these results and named after them,
        because a component means nothing without the table saying what
        loads on it.
    pca_components : int, default=0
        How many components to keep, when reducing. ``0`` decides by the
        Kaiser criterion -- a starting point, not an answer.
    pca_retain : {"parallel", "kaiser"}, default="parallel"
        How the component count is chosen when ``pca_components`` is 0:
        parallel analysis keeps a component while its eigenvalue beats what
        random data of the same size produce at that rank; the Kaiser rule
        keeps every eigenvalue above 1, which on a wide table is most of them.
    pca_rotation : bool, default=True
        Rotate the components (varimax) so each loads on a small cluster of
        features and is therefore nameable.
    pca_max_missing : float, default=0.5
        When reducing, set aside any feature missing for more than this
        fraction of the rows rather than letting it delete them -- the same
        rule, and the same default, the prediction steps apply to predictors.
        Of the rows that remain, the components are fitted on those with
        every kept feature present and a row still missing one is left
        unscored rather than guessed at; both counts are reported.
    split_col
        Run the whole analysis once per value of this column, labeling each
        result with it. Set to ``"source_col"`` when a spreadsheet's text
        columns were measured separately: a participant then has one row per
        column, and analyzing those together would count one person's
        several answers as several independent observations -- inflating the
        degrees of freedom, and every p-value with them.
    welch
        False (default): classic pooled ANOVA. True: Welch's F, for when the
        groups' variances should not be pooled.
    posthoc : {"auto", "tukey", "games_howell", "bonferroni", "none"}, default="auto"
        Which pairwise comparison to run, and therefore how the pairwise
        p-values are corrected:

        - ``"auto"`` (default) -- Tukey-Kramer under a pooled ANOVA,
          Games-Howell under Welch. The pairing matters: Tukey's statistic
          divides by a pooled error term, which is incoherent once you have
          chosen Welch precisely because you distrust pooling.
        - ``"tukey"`` / ``"games_howell"`` -- force one of them.
        - ``"bonferroni"`` -- pairwise t-tests, Bonferroni-corrected.
        - ``"none"`` -- pairwise t-tests, uncorrected. Six pairs at .05 give
          a ~26% chance of at least one false positive, so this is a
          deliberate choice, not a shortcut.
    p_adjust : {"none", "fdr_bh", "fdr_by", "holm", "bonferroni"}, default="fdr_bh"
        How the *omnibus* p-values are adjusted across features, which is a
        separate family from the pairwise one. ``"fdr_by"`` is valid under
        any dependence between features, and stricter for it; ``"holm"`` and
        ``"bonferroni"`` control the family-wise error rate instead. With
        ``"none"`` no ``p_adj`` column is written -- an unadjusted number
        under an "adjusted" heading is worse than no column.
    alpha
        The threshold used for the ``pairwise_sig`` summary and the post-hoc
        confidence intervals. It does not gate what is written -- every test
        is in the tables.
    out_dir
        Where the two CSVs go; defaults to the analysis table's folder.
    overwrite_existing
        When False (default) and the results already exist, they are
        returned untouched instead of recomputed.
    rounding
        Decimal places in the output tables. Values too small for it (a p of
        3e-12) switch to significant digits rather than collapsing to 0.

    Returns
    -------
    Path
        ``group_differences.csv`` (one row per feature; the pairwise detail
        is in ``group_differences_pairwise.csv`` beside it).

    Notes
    -----
    A p-value printed as ``0`` means the tail underflowed double precision
    (roughly p < 1e-300), not that the probability is zero. The columns stay
    numeric so R, Excel and pandas read them as numbers; report such a value
    the way every statistics package does, as "p < .001".
    """
    import csv

    import numpy as np

    table_csv = Path(table_csv)
    folder = Path(out_dir) if out_dir else table_csv.parent
    folder.mkdir(parents=True, exist_ok=True)
    main_path = folder / "group_differences.csv"
    pair_path = folder / "group_differences_pairwise.csv"
    if reusable(main_path, table_csv, overwrite_existing=overwrite_existing,
                verbose=verbose, what="the group differences"):
        if verbose:
            print(f"Group differences already exist; returning existing "
                  f"file: {main_path}")
        return main_path

    announce(on_progress, "reading the analysis table")
    from ._common import read_str_csv

    table = read_str_csv(table_csv, encoding=encoding)
    if group_col not in table.columns:
        raise ValueError(
            f"group column {group_col!r} is not in the analysis table, "
            f"which has {len(table.columns)}: "
            f"{name_a_few(list(table.columns))}.")

    labels_all = np.array([str(v).strip() for v in table[group_col].tolist()])
    levels = sorted({v for v in labels_all if v})
    if len(levels) < 2:
        raise ValueError(
            f"column {group_col!r} holds "
            f"{'one group' if levels else 'no group labels'} "
            f"({levels or 'all blank'}); comparing groups needs at least two.")

    # a control isn't a feature. comparing groups on age while adjusting for
    # age would just hand us a row of zeros dressed up as a finding.
    feature_cols = default_feature_cols(
        table, table_csv,
        exclude=(group_col,) + tuple(str(c) for c in control_cols))
    sets = resolve_feature_sets(feature_sets, table_csv=table_csv,
                                feature_cols=feature_cols)

    from . import _controls

    control_matrix, control_names, control_notes, control_spec = \
        _controls.build(table, control_cols, categorical_controls)
    has_controls = bool(control_names)
    if has_controls and welch:
        raise ValueError(
            "welch=True cannot be combined with control columns: Welch's F "
            "exists precisely to avoid pooling the error term, and an "
            "analysis of covariance pools it by construction. Drop the "
            "controls, or drop welch.")

    method = "ancova" if has_controls else ("welch" if welch else "anova")
    # unset reads as the documented default, for the same reason as p_adjust
    # below: a null in a preset means "I didn't choose", and not choosing
    # ought to get you the coherent pairing (Tukey with a pooled F,
    # Games-Howell with Welch).
    posthoc = str(posthoc or "auto").strip().lower()
    if posthoc not in ("auto", "tukey", "games_howell", "bonferroni", "none"):
        raise ValueError(
            f"posthoc must be 'auto', 'tukey', 'games_howell', 'bonferroni' "
            f"or 'none', got {posthoc!r}")
    if posthoc == "auto":
        posthoc = "games_howell" if welch else "tukey"
    if has_controls and posthoc == "games_howell":
        raise ValueError(
            "posthoc='games_howell' cannot be combined with control "
            "columns, for the same reason as welch: it is built for "
            "unpooled variances and an analysis of covariance pools them.")
    posthoc_name = {"tukey": "tukey_hsd", "games_howell": "games_howell",
                    "bonferroni": "t_bonferroni", "none": "t_uncorrected"}[posthoc]
    # a preset that says `p_adjust:` with nothing after it means "I didn't
    # set this", so we fall back to the default -- it does NOT mean "correct
    # nothing". only the explicit string "none" turns correction off.
    p_adjust = "fdr_bh" if p_adjust in (None, "") else str(p_adjust)
    if p_adjust not in P_ADJUST_METHODS:
        raise ValueError(
            f"unknown p-value adjustment {p_adjust!r}; choose one of "
            f"{sorted(P_ADJUST_METHODS)}")
    adjusting = p_adjust != "none"
    desc_header = [f"{stat}_{level}" for level in levels
                   for stat in ("n", "mean", "sd")]
    if has_controls:
        # the adjusted means go next to the raw ones rather than replacing
        # them. the difference between the two columns is what the controls
        # did, and hiding the raw mean would hide that too.
        desc_header += [f"adj_mean_{level}" for level in levels]
    lead = [split_col] if split_col else []
    main_header = (lead + ["feature_set", "feature"] + desc_header +
                   ["F", "df1", "df2", "p"] + (["p_adj"] if adjusting else []) +
                   ["eta2", "method", "pairwise_sig", "note"])
    # we name this column for what's actually in it: under posthoc="none" the
    # pairwise p is an uncorrected t-test p, and calling it `p_adj` (like we
    # used to) said otherwise.
    pair_header = lead + ["feature_set", "feature", "group_1", "group_2",
                          "n_1", "n_2", "mean_diff", "ci_low", "ci_high",
                          "d", "p" if posthoc == "none" else "p_adj", "method"]

    main_rows, pair_rows = [], []
    total = sum(len(cols) for cols in sets.values())
    done = 0
    n_tested = 0
    top_effects = []           # (eta2, set, feature, sig_summary)

    # we parse the numeric columns once and mask them per subset rather than
    # re-parsing for each one. parsing is the expensive part, and the subsets
    # are all just views of the same rows anyway.
    columns = {c: numeric_column(table[c].tolist(), column=c)
               for cols in sets.values() for c in cols}
    # swap in components for raw measures, if this analysis asked for them --
    # we do this per analysis and per feature set, since the right answer
    # differs between them.
    sets, columns, pca_notes, _reductions = reduce_sets(
        sets, columns, pca=pca, n_components=pca_components,
        rotation=pca_rotation, out_stem=folder / "group_differences",
        encoding=encoding, rounding=rounding, verbose=verbose,
        max_missing=pca_max_missing, retain=pca_retain)

    subsets = list(split_rows(table, split_col))
    total *= len(subsets)

    for subset, subset_mask in subsets:
        prefix = [subset] if split_col else []
        for set_name, cols in sets.items():
            omnibus_ps = []
            set_rows = []
            for feature in cols:
                done += 1
                if on_progress is not None:
                    on_progress(done, total, f"testing {feature}")
                values_all = columns[feature]
                valid = ~np.isnan(values_all) & (labels_all != "")
                if has_controls:
                    valid = valid & ~np.isnan(control_matrix).any(axis=1)
                if subset_mask is not None:
                    valid = valid & subset_mask
                values, labels = values_all[valid], labels_all[valid]

                desc = _descriptives(values, labels, levels)
                desc_cells = []
                for level in levels:
                    n, mean, sd = desc[level]
                    desc_cells += [str(n), fmt(mean, rounding), fmt(sd, rounding)]

                blank_adjusted = [""] * len(levels) if has_controls else []
                usable = [lv for lv in levels if desc[lv][0] >= 2]
                note = ""
                excluded = [lv for lv in levels if lv not in usable]
                if excluded:
                    note = f"excluded (n<2): {', '.join(excluded)}"

                if len(usable) < 2:
                    omnibus_ps.append(float("nan"))
                    why = (note + "; " if note else "") + \
                        "fewer than 2 groups with data"
                    set_rows.append(prefix + [set_name, feature] + desc_cells
                                    + blank_adjusted
                                    + [""] * (6 if adjusting else 5)
                                    + [method, "", why])
                    continue

                stats = [(desc[lv][0], desc[lv][1], desc[lv][2] ** 2)
                         for lv in usable]
                if np.all(values[np.isin(labels, usable)] ==
                          values[np.isin(labels, usable)][0]):
                    omnibus_ps.append(float("nan"))
                    row = (prefix + [set_name, feature] + desc_cells
                           + blank_adjusted
                           + [""] * (6 if adjusting else 5)
                           + [method, "",
                              (note + "; " if note else "") + "constant feature"])
                    set_rows.append(row)
                    continue

                fit = None
                if has_controls:
                    keep = np.isin(labels, usable)
                    fit = _ancova(values[keep], labels[keep], usable,
                                  control_matrix[valid][keep])
                    if fit is None:
                        omnibus_ps.append(float("nan"))
                        set_rows.append(
                            prefix + [set_name, feature] + desc_cells
                            + blank_adjusted + [""] * (6 if adjusting else 5)
                            + [method, "",
                               (note + "; " if note else "")
                               + "too few rows for the controls"])
                        continue
                    f_value, df1, df2 = fit["f"], fit["df1"], fit["df2"]
                    p, eta2 = fit["p"], fit["eta2"]
                    ms_within, df_within = None, None
                    # now we fill in the adjusted-mean cells where the blanks were.
                    adjusted_cells = [
                        fmt(fit["adjusted"].get(lv, float("nan")), rounding)
                        for lv in levels]
                    desc_cells = desc_cells + adjusted_cells
                elif welch:
                    f_value, df1, df2, p = _welch_anova(stats)
                    ms_within, df_within = None, None
                    eta2 = _eta_squared(stats)
                else:
                    f_value, df1, df2, p, ms_within = _classic_anova(stats)
                    df_within = df2
                    eta2 = _eta_squared(stats)
                omnibus_ps.append(p)
                n_tested += 1

                k = len(usable)
                sig_parts = []
                all_pairs = _pairs(usable)
                for a, b in all_pairs:
                    if fit is not None:
                        diff, lo, hi, p_adj = _ancova_pair(
                            fit, usable, a, b, n_pairs=len(all_pairs),
                            posthoc=posthoc, alpha=alpha)
                    elif posthoc == "games_howell":
                        diff, lo, hi, p_adj = _games_howell_pair(
                            desc, a, b, k=k, alpha=alpha)
                    elif posthoc == "tukey":
                        # Tukey needs a pooled error term, and Welch's omnibus
                        # never computed one -- so we compute it here, so that
                        # the choice is there whichever F was asked for.
                        pooled = ms_within if ms_within is not None else \
                            _pooled_ms(stats)
                        pooled_df = df_within if df_within is not None else \
                            sum(n for n, _, _ in stats) - len(stats)
                        diff, lo, hi, p_adj = _tukey_pair(
                            desc, a, b, k=k, ms_within=pooled,
                            df_within=pooled_df, alpha=alpha)
                    else:
                        diff, lo, hi, p_adj = _pairwise_t(
                            desc, a, b, welch=welch,
                            ms_within=(ms_within if ms_within is not None
                                       else _pooled_ms(stats)),
                            df_within=(df_within if df_within is not None
                                       else sum(n for n, _, _ in stats) - len(stats)),
                            n_pairs=len(all_pairs),
                            bonferroni=(posthoc == "bonferroni"), alpha=alpha)
                    d = _cohens_d(desc, a, b)
                    pair_rows.append(prefix + [
                        set_name, feature, a, b, str(desc[a][0]),
                        str(desc[b][0]), fmt(diff, rounding), fmt(lo, rounding),
                        fmt(hi, rounding), fmt(d, rounding), fmt(p_adj, rounding),
                        posthoc_name])
                    if p_adj < alpha:
                        sig_parts.append(f"{a}>{b}" if diff > 0 else f"{b}>{a}")

                set_rows.append(
                    prefix + [set_name, feature] + desc_cells
                    + [fmt(f_value, rounding), fmt(df1, rounding),
                       fmt(df2, rounding), fmt(p, rounding)]
                    + ([""] if adjusting else [])   # p_adj; we fill this in below
                    + [fmt(eta2, rounding), method, "; ".join(sig_parts), note])
                top_effects.append(
                    (eta2, f"{set_name} @ {subset}" if subset else set_name,
                     feature, "; ".join(sig_parts) or "no pair survived"))

            # the adjustment family is the set. each feature set is one batch
            # of hypotheses, and correcting across sets would punish someone
            # for having asked a second question.
            if adjusting:
                adjusted = adjust_pvalues(omnibus_ps, p_adjust)
                adj_col = main_header.index("p_adj")
                for row, value in zip(set_rows, adjusted, strict=True):
                    row[adj_col] = fmt(value, rounding)
            main_rows.extend(set_rows)

    announce(on_progress, "writing the results")
    with atomic_write(main_path, mode="w", encoding=encoding,
                      newline="") as fh:
        w = csv.writer(fh)
        w.writerow(main_header)
        w.writerows(main_rows)
    with atomic_write(pair_path, mode="w", encoding=encoding,
                      newline="") as fh:
        w = csv.writer(fh)
        w.writerow(pair_header)
        w.writerows(pair_rows)

    verdict_col = main_header.index("p_adj" if adjusting else "p")
    survivors = sum(1 for row in main_rows
                    if row[verdict_col] != ""
                    and float(row[verdict_col]) < alpha)
    section = _section_md(
        group_col=group_col, levels=levels, method=method,
        posthoc=posthoc_name, n_tested=n_tested, survivors=survivors,
        alpha=alpha, top_effects=top_effects, sets=sets,
        p_adjust=p_adjust, subsets=[s for s, _ in subsets if s],
        split_col=split_col, control_names=control_names,
        control_notes=control_notes)
    section += components_appendix(pca_notes)
    write_section(folder, "group-differences", section)
    if verbose:
        print(f"[group-differences] {n_tested} feature(s) tested across "
              f"{len(levels)} group(s); {survivors} survive FDR "
              f"-> {main_path}")
    return main_path

taters.stats.correlations

Correlations between features and outcomes, laid out for perusal.

Features as rows, outcomes as column blocks -- for each outcome its r, its p, its FDR-adjusted p, and the pairwise N -- because that is how a researcher actually reads this table: scan down an outcome's column for the big coefficients, then check that the N behind each one is respectable.

Missing data is handled by pairwise deletion: each (feature, outcome) cell uses exactly the rows where both values exist, and the _n column says how many that was. NA is never zero, and a cell whose pairwise-complete subset is too small (n < 3) or constant is left honestly blank -- with its _n still filled, so the blank explains itself.

Spearman coefficients rank each pair's complete subset (average ranks for ties), matching scipy.stats.spearmanr on that subset; a global ranking would let rows missing from one pair distort another's coefficients.

analyze_correlations

analyze_correlations(
    *,
    table_csv,
    outcome_cols,
    feature_sets=None,
    control_cols=(),
    categorical_controls=(),
    pca="off",
    pca_components=0,
    pca_retain="parallel",
    pca_rotation=True,
    pca_max_missing=MAX_MISSING,
    split_col=None,
    method="pearson",
    p_adjust="fdr_bh",
    out_dir=None,
    overwrite_existing=False,
    on_progress=None,
    verbose=True,
    encoding="utf-8-sig",
    rounding=4
)

Correlate every feature with every outcome column.

Parameters:

Name	Type	Description	Default
table_csv	PathLike	The assembled analysis table.	required
outcome_cols	Sequence[str]	The outcome columns (numeric). Each becomes a block of four columns in the output: <outcome>_r, _p, _p_adj, _n.	required
feature_sets		None for all features as one family, "per_table" to repeat per source feature table, or {name: [columns]}.	None
control_cols	Sequence[str]	Columns to hold constant, turning every coefficient into a partial correlation: the relationship that is left once these have had their say. Continuous controls (age, word count) enter as themselves; categorical ones (gender, site) are expanded into indicators against a reference level, which the report names. The degrees of freedom fall by one per control column, and rows missing any control drop out -- so the _n columns can be smaller than they were, which is the honest cost of adjusting.	()
categorical_controls	Sequence[str]	Which of control_cols to treat as categories even though they look numeric -- a site coded 1/2/3, a condition coded 0/1/2. Left to itself, a column whose every value parses as a number is treated as continuous.	()
pca	str or list of str	Analyze components instead of the raw measures. "off" uses the features as they are; "all" reduces every feature set; a list of set names reduces those and leaves the rest alone -- which is the case worth having, since a hundred dictionary categories are worth reducing and eight readability indices are not. Set per analysis, deliberately: raw variables read better in a correlation table, where each row is a measure you can name, while a ridge over four hundred collinear measures is what components are for. The loadings land beside these results and named after them, because a component means nothing without the table saying what loads on it.	"off"
pca_components	int	How many components to keep, when reducing. 0 decides by the Kaiser criterion -- a starting point, not an answer.	0
pca_retain	('parallel', 'kaiser')	How the component count is chosen when pca_components is 0: parallel analysis keeps a component while its eigenvalue beats what random data of the same size produce at that rank; the Kaiser rule keeps every eigenvalue above 1, which on a wide table is most of them.	"parallel"
pca_rotation	bool	Rotate the components (varimax) so each loads on a small cluster of features and is therefore nameable.	True
pca_max_missing	float	When reducing, set aside any feature missing for more than this fraction of the rows rather than letting it delete them -- the same rule, and the same default, the prediction steps apply to predictors. Of the rows that remain, the components are fitted on those with every kept feature present and a row still missing one is left unscored rather than guessed at; both counts are reported.	0.5
split_col	Optional[str]	Run the whole analysis once per value of this column, labeling each result with it. Set to "source_col" when a spreadsheet's text columns were measured separately: a participant then has one row per column, and analyzing those together would count one person's several answers as several independent observations -- inflating the degrees of freedom, and every p-value with them.	None
method	('pearson', 'spearman', 'both')	"both" writes two files and returns the Pearson one. The type is spelled out here because it is what the options screen reads: a closed set becomes a picker, and without one the user is left typing into a box with no way to learn that "spearman" is a word it knows.	"pearson"
p_adjust	('none', 'fdr_bh', 'fdr_by', 'holm', 'bonferroni')	How p-values are adjusted for the number of correlations tested: "fdr_by" is valid under any dependence between features, and stricter for it; "holm" and "bonferroni" control the family-wise error rate instead; "none" corrects nothing. The family is one feature set's whole feature-by-outcome grid, per method. With "none" no _p_adj columns are written -- an unadjusted number under an "adjusted" heading is worse than no column.	"none"
out_dir	Optional[PathLike]	Where the CSV(s) go; defaults to the analysis table's folder.	None
overwrite_existing	bool	When False (default) and the results already exist, they are returned untouched instead of recomputed.	False
rounding	int	Decimal places in the output tables. Values too small for it (a p of 3e-12) switch to significant digits rather than collapsing to 0.	4

Returns:

Type	Description
Path	correlations_<method>.csv (the Pearson file under "both").

Notes

A p-value printed as 0 means the tail underflowed double precision (roughly p < 1e-300), not that the probability is zero. The columns stay numeric so R, Excel and pandas read them as numbers; report such a value the way every statistics package does, as "p < .001".

Source code in src\taters\stats\correlations.py

def analyze_correlations(
    *,
    table_csv: PathLike,
    outcome_cols: Sequence[str],
    feature_sets=None,
    control_cols: Sequence[str] = (),
    categorical_controls: Sequence[str] = (),
    pca: object = "off",
    pca_components: int = 0,
    pca_retain: Literal["parallel", "kaiser"] = "parallel",
    pca_rotation: bool = True,
    pca_max_missing: float = MAX_MISSING,
    split_col: Optional[str] = None,
    method: str = "pearson",
    p_adjust: str = "fdr_bh",
    out_dir: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[..., None]] = None,
    verbose: bool = True,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
) -> Path:
    """
    Correlate every feature with every outcome column.

    Parameters
    ----------
    table_csv
        The assembled analysis table.
    outcome_cols
        The outcome columns (numeric). Each becomes a block of four columns
        in the output: ``<outcome>_r``, ``_p``, ``_p_adj``, ``_n``.
    feature_sets
        ``None`` for all features as one family, ``"per_table"`` to repeat
        per source feature table, or ``{name: [columns]}``.
    control_cols
        Columns to hold constant, turning every coefficient into a **partial**
        correlation: the relationship that is left once these have had their
        say. Continuous controls (age, word count) enter as themselves;
        categorical ones (gender, site) are expanded into indicators against
        a reference level, which the report names. The degrees of freedom
        fall by one per control column, and rows missing any control drop
        out -- so the ``_n`` columns can be smaller than they were, which is
        the honest cost of adjusting.
    categorical_controls
        Which of ``control_cols`` to treat as categories even though they
        look numeric -- a site coded 1/2/3, a condition coded 0/1/2. Left to
        itself, a column whose every value parses as a number is treated as
        continuous.
    pca : str or list of str, default="off"
        Analyze components instead of the raw measures. ``"off"`` uses the
        features as they are; ``"all"`` reduces every feature set; a list of
        set names reduces those and leaves the rest alone -- which is the
        case worth having, since a hundred dictionary categories are worth
        reducing and eight readability indices are not.

        Set per analysis, deliberately: raw variables read better in a
        correlation table, where each row is a measure you can name, while a
        ridge over four hundred collinear measures is what components are
        for. The loadings land beside these results and named after them,
        because a component means nothing without the table saying what
        loads on it.
    pca_components : int, default=0
        How many components to keep, when reducing. ``0`` decides by the
        Kaiser criterion -- a starting point, not an answer.
    pca_retain : {"parallel", "kaiser"}, default="parallel"
        How the component count is chosen when ``pca_components`` is 0:
        parallel analysis keeps a component while its eigenvalue beats what
        random data of the same size produce at that rank; the Kaiser rule
        keeps every eigenvalue above 1, which on a wide table is most of them.
    pca_rotation : bool, default=True
        Rotate the components (varimax) so each loads on a small cluster of
        features and is therefore nameable.
    pca_max_missing : float, default=0.5
        When reducing, set aside any feature missing for more than this
        fraction of the rows rather than letting it delete them -- the same
        rule, and the same default, the prediction steps apply to predictors.
        Of the rows that remain, the components are fitted on those with
        every kept feature present and a row still missing one is left
        unscored rather than guessed at; both counts are reported.
    split_col
        Run the whole analysis once per value of this column, labeling each
        result with it. Set to ``"source_col"`` when a spreadsheet's text
        columns were measured separately: a participant then has one row per
        column, and analyzing those together would count one person's
        several answers as several independent observations -- inflating the
        degrees of freedom, and every p-value with them.
    method : {"pearson", "spearman", "both"}, default="pearson"
        ``"both"`` writes two files and returns the Pearson one. The type is
        spelled out here because it is what the options screen reads: a
        closed set becomes a picker, and without one the user is left typing
        into a box with no way to learn that "spearman" is a word it knows.
    p_adjust : {"none", "fdr_bh", "fdr_by", "holm", "bonferroni"}, default="fdr_bh"
        How p-values are adjusted for the number of correlations tested:
        ``"fdr_by"`` is valid under any dependence between features, and
        stricter for it; ``"holm"`` and ``"bonferroni"`` control the
        family-wise error rate instead; ``"none"`` corrects nothing. The family is one feature set's
        whole feature-by-outcome grid, per method. With ``"none"`` no
        ``_p_adj`` columns are written -- an unadjusted number under an
        "adjusted" heading is worse than no column.
    out_dir
        Where the CSV(s) go; defaults to the analysis table's folder.
    overwrite_existing
        When False (default) and the results already exist, they are
        returned untouched instead of recomputed.
    rounding
        Decimal places in the output tables. Values too small for it (a p of
        3e-12) switch to significant digits rather than collapsing to 0.

    Returns
    -------
    Path
        ``correlations_<method>.csv`` (the Pearson file under ``"both"``).

    Notes
    -----
    A p-value printed as ``0`` means the tail underflowed double precision
    (roughly p < 1e-300), not that the probability is zero. The columns stay
    numeric so R, Excel and pandas read them as numbers; report such a value
    the way every statistics package does, as "p < .001".
    """
    import csv

    method = str(method).strip().lower()
    if method not in ("pearson", "spearman", "both"):
        raise ValueError(
            f"method must be 'pearson', 'spearman' or 'both', got {method!r}")
    outcome_cols = [str(c) for c in outcome_cols]
    if not outcome_cols:
        raise ValueError("outcome_cols is empty: nothing to correlate with")
    # a preset that says `p_adjust:` with nothing after it means "I didn't
    # set this", so we fall back to the default -- it does NOT mean "correct
    # nothing". only the explicit string "none" turns correction off.
    p_adjust = "fdr_bh" if p_adjust in (None, "") else str(p_adjust)
    if p_adjust not in P_ADJUST_METHODS:
        raise ValueError(
            f"unknown p-value adjustment {p_adjust!r}; choose one of "
            f"{sorted(P_ADJUST_METHODS)}")
    adjusting = p_adjust != "none"

    table_csv = Path(table_csv)
    folder = Path(out_dir) if out_dir else table_csv.parent
    folder.mkdir(parents=True, exist_ok=True)
    methods = ("pearson", "spearman") if method == "both" else (method,)
    primary = folder / f"correlations_{methods[0]}.csv"
    if reusable(primary, table_csv, overwrite_existing=overwrite_existing,
                verbose=verbose, what="the correlations"):
        if verbose:
            print(f"Correlations already exist; returning existing file: "
                  f"{primary}")
        return primary

    announce(on_progress, "reading the analysis table")
    from ._common import read_str_csv

    table = read_str_csv(table_csv, encoding=encoding)
    missing = [c for c in outcome_cols if c not in table.columns]
    if missing:
        raise ValueError(
            f"outcome column(s) not in the analysis table: "
            f"{name_a_few(missing)}. It has {len(table.columns)}: "
            f"{name_a_few(list(table.columns))}.")

    # a control isn't a feature. if we left it in, `age` would get correlated
    # with the outcome while age was being partialled out of it -- we'd get a
    # row of zeros dressed up as a finding.
    feature_cols = default_feature_cols(
        table, table_csv,
        exclude=tuple(outcome_cols) + tuple(str(c) for c in control_cols))
    sets = resolve_feature_sets(feature_sets, table_csv=table_csv,
                                feature_cols=feature_cols)

    from . import _controls

    control_matrix, control_names, control_notes, control_spec = \
        _controls.build(table, control_cols, categorical_controls)

    outcomes = {c: numeric_column(table[c].tolist(), column=c)
                for c in outcome_cols}
    features = {c: numeric_column(table[c].tolist(), column=c)
                for cols in sets.values() for c in cols}
    # swap in components for raw measures, if this analysis asked for them --
    # we do this per analysis and per feature set, since the right answer
    # differs between them.
    sets, features, pca_notes, _reductions = reduce_sets(
        sets, features, pca=pca, n_components=pca_components,
        rotation=pca_rotation, out_stem=folder / "correlations",
        encoding=encoding, rounding=rounding, verbose=verbose,
        max_missing=pca_max_missing, retain=pca_retain)


    lead = [split_col] if split_col else []
    subsets = list(split_rows(table, split_col))
    header = lead + ["feature_set", "feature"]
    for outcome in outcome_cols:
        header += [f"{outcome}_r", f"{outcome}_p"]
        if adjusting:
            header.append(f"{outcome}_p_adj")
        header.append(f"{outcome}_n")

    total = (sum(len(cols) for cols in sets.values()) * len(methods)
             * len(subsets))
    done = 0
    strongest = {}                 # method -> [(abs r, r, set, feature, outcome, n)]
    written = []
    for m in methods:
        rows = []
        for subset, mask in subsets:
          prefix = [subset] if split_col else []
          for set_name, cols in sets.items():
            cells = {}             # (feature, outcome) -> (r, p, n)
            for feature in cols:
                done += 1
                if on_progress is not None:
                    on_progress(done, total, f"correlating {feature}")
                for outcome in outcome_cols:
                    x = features[feature]
                    y = outcomes[outcome]
                    c = control_matrix
                    if mask is not None:
                        x, y = x[mask], y[mask]
                        c = c[mask] if c.shape[1] else c
                    cells[(feature, outcome)] = _pair_r(
                        x, y, spearman=(m == "spearman"), controls=c)
            # one correction family per (set, method). every feature x
            # outcome test goes in it, since those are the tests we're
            # scanning across.
            keys = [(f, o) for f in cols for o in outcome_cols]
            adjusted = dict(zip(
                keys, adjust_pvalues([cells[k][1] for k in keys], p_adjust),
                strict=True))
            for feature in cols:
                row = prefix + [set_name, feature]
                for outcome in outcome_cols:
                    r, p, n = cells[(feature, outcome)]
                    q = adjusted[(feature, outcome)]
                    row += [fmt(r, rounding), fmt(p, rounding)]
                    if adjusting:
                        row.append(fmt(q, rounding))
                    row.append(str(n))
                    if r == r:     # i.e. not NaN
                        strongest.setdefault(m, []).append(
                            (abs(r), r,
                             f"{set_name} @ {subset}" if subset else set_name,
                             feature, outcome, n, q))
                rows.append(row)
        path = folder / f"correlations_{m}.csv"
        with atomic_write(path, mode="w", encoding=encoding, newline="") as fh:
            w = csv.writer(fh)
            w.writerow(header)
            w.writerows(rows)
        written.append(path)

    section = _section_md(
        methods=methods, outcome_cols=outcome_cols, sets=sets,
        strongest=strongest, files=[p.name for p in written],
        p_adjust=p_adjust, subsets=[s for s, _ in subsets if s],
        split_col=split_col, control_names=control_names,
        control_notes=control_notes)
    section += components_appendix(pca_notes)
    write_section(folder, "correlations", section)
    if verbose:
        print(f"[correlations] {len(features)} feature(s) × "
              f"{len(outcome_cols)} outcome(s), {', '.join(methods)} "
              f"-> {primary}")
    return primary

taters.stats.ridge

Ridge regression with cross-validation: predict an outcome from the features, honestly, and keep the model.

Language features are many and collinear -- three readability indices are three views of sentence length -- which is exactly the situation ordinary least squares handles worst and ridge handles well: the penalty trades a little bias for a large drop in variance, and it never has to invert a singular matrix. What ridge does not do is choose its own penalty, so this searches a grid of alphas by k-fold cross-validation and reports the out-of-fold performance at the chosen one. In-sample R-squared is reported too, next to it, because the gap between them is the whole story of whether a model learned anything or memorized.

The arithmetic is plain numpy on purpose. One economy SVD per fold serves the entire alpha grid -- coefficients for any alpha are V diag(s/(s^2+a)) U' y -- so a 33-point grid costs one decomposition rather than 33 solves, and the numbers are exactly reproducible from the stored model rather than depending on a solver's iteration count.

A fit is an instrument, not a result: the model file carries the predictor names, the training centering and scaling, the chosen alpha and the grid it came from, and the coefficients, so the same model can score a new dataset later -- standardized against the training sample, which is what makes the scores comparable. That is the same fit-once/apply-many discipline the MEM topic model and PCA use here.

apply_ridge_csv

apply_ridge_csv(
    *,
    model_json,
    input_csv,
    out_csv=None,
    overwrite_existing=False,
    on_progress=None,
    verbose=True,
    encoding="utf-8-sig",
    rounding=4,
    id_cols=None
)

Score a new table with a saved ridge model.

Predictors are matched by name, in any column order; every other column passes through as an identifier. A predictor the new table lacks is refused by name rather than imputed -- a model scored on a feature it never sees is not the model that was validated.

Parameters:

Name	Type	Description	Default
model_json	PathLike	A model written by :func:fit_ridge_csv (or a folder holding exactly one).	required
input_csv	PathLike	The table to score: an analysis table, or any feature table carrying the model's predictors.	required
out_csv	Optional[PathLike]	Defaults to <input stem>_ridge_predictions.csv beside the input.	None
overwrite_existing	bool	When False (default) and the predictions already exist, they are returned untouched instead of recomputed.	False
rounding	int	Decimal places in the predicted values.	4

Returns:

Type	Description
Path	The predictions table: the input's identifier columns plus one pred_<outcome> column per outcome the model holds.

Source code in src\taters\stats\ridge.py

def apply_ridge_csv(
    *,
    model_json: PathLike,
    input_csv: PathLike,
    out_csv: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[..., None]] = None,
    verbose: bool = True,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
    id_cols: Optional[Sequence[str]] = None,
) -> Path:
    """
    Score a new table with a saved ridge model.

    Predictors are matched **by name**, in any column order; every other
    column passes through as an identifier. A predictor the new table lacks
    is refused by name rather than imputed -- a model scored on a feature it
    never sees is not the model that was validated.

    Parameters
    ----------
    model_json
        A model written by :func:`fit_ridge_csv` (or a folder holding exactly
        one).
    input_csv
        The table to score: an analysis table, or any feature table carrying
        the model's predictors.
    out_csv
        Defaults to ``<input stem>_ridge_predictions.csv`` beside the input.
    overwrite_existing
        When False (default) and the predictions already exist, they are
        returned untouched instead of recomputed.
    rounding
        Decimal places in the predicted values.

    Returns
    -------
    Path
        The predictions table: the input's identifier columns plus one
        ``pred_<outcome>`` column per outcome the model holds.
    """
    import numpy as np

    got = prepare_apply(
        model_json=model_json, load=_load_model, input_csv=input_csv,
        out_csv=out_csv, default_suffix="_ridge_predictions.csv",
        overwrite_existing=overwrite_existing, on_progress=on_progress,
        verbose=verbose, encoding=encoding, id_cols=id_cols)
    if isinstance(got, Path):
        return got
    model, design = got.model, got.design
    columns = {}
    skipped = {}
    for outcome, block in model["outcomes"].items():
        preds = predict_rows(design, block["kept"],
                             np.asarray(block["mu"]),
                             np.asarray(block["sigma"]),
                             np.asarray(block["coef"]),
                             float(block["intercept"]), rounding)
        label = output_label(model, outcome)
        columns[f"pred_{label}"] = preds
        skipped[label] = sum(1 for p in preds if p is None)

    return write_predictions(got, columns, skipped, encoding=encoding,
                             on_progress=on_progress, verbose=verbose,
                             tag="ridge")

default_alphas

default_alphas()

The default penalty grid, as a plain list (JSON-storable).

Source code in src\taters\stats\ridge.py

def default_alphas() -> List[float]:
    """The default penalty grid, as a plain list (JSON-storable)."""
    import numpy as np

    return [float(a) for a in np.logspace(_DEFAULT_ALPHA_LOG_LOW,
                                          _DEFAULT_ALPHA_LOG_HIGH,
                                          _DEFAULT_ALPHA_COUNT)]

fit_ridge_csv

fit_ridge_csv(
    *,
    table_csv,
    outcome_cols,
    feature_sets=None,
    control_cols=(),
    categorical_controls=(),
    control_combos="none_and_all",
    set_combos="subsets",
    pca="off",
    pca_components=0,
    pca_retain="parallel",
    pca_rotation=True,
    split_col=None,
    alphas=None,
    n_folds=5,
    stratify=True,
    zscore=True,
    max_missing=MAX_MISSING,
    seed=0,
    out_dir=None,
    out_models_dir=None,
    overwrite_existing=False,
    on_progress=None,
    verbose=True,
    encoding="utf-8-sig",
    rounding=4
)

Fit a cross-validated ridge model per outcome, and keep it.

Parameters:

Name	Type	Description	Default
table_csv	PathLike	The assembled analysis table.	required
outcome_cols	Sequence[str]	The numeric column(s) to predict. Each gets its own model, its own chosen penalty and its own row of performance.	required
feature_sets		None to use every feature as one predictor set, "per_table" to fit a separate model per source feature table -- which is how you answer "which feature set actually predicts this?", since the resulting rows sit side by side in one metrics table -- or an explicit {name: [columns]}.	None
control_cols	Sequence[str]	Columns to fit alongside -- and instead of -- the language, so the language's own contribution can be read off. With controls the run fits three models per outcome on one common sample: the controls alone, the language alone, and both; the combined model's row carries delta_r2_over_controls, which is what "language adds this much over age and gender" means. Continuous controls enter as themselves, categorical ones as indicators against a reference level. The shared sample is the part that matters: fitting each model on whatever rows it happened to have would make the difference between them part sample and part language, with no way to tell which.	()
categorical_controls	Sequence[str]	Which of control_cols to treat as categories despite looking numeric.	()
set_combos	('none', 'each_and_all', 'subsets')	When the feature tables are analyzed together (feature_sets not given), how to combine them: subsets fits every combination of the tables -- each alone, every pair, and so on up to all of them -- so the report can say what each table adds to the others; each_and_all fits each alone and all together; none fits only all together. Every combination is fitted for up to five tables; above that, each and all. A combination is named by its members (dictionary+readability) in the tables and the models.	"none"
control_combos	('none_and_all', 'each', 'subsets')	How many control sets to try. The default pair answers the usual question. "each" adds one model per single control, which is how you see which of them is doing the work; "subsets" runs all 2^k of them, which is thorough and grows the table accordingly.	"none_and_all"
pca	str or list of str	Analyze components instead of the raw measures. "off" uses the features as they are; "all" reduces every feature set; a list of set names reduces those and leaves the rest alone -- which is the case worth having, since a hundred dictionary categories are worth reducing and eight readability indices are not. Set per analysis, deliberately: raw variables read better in a correlation table, where each row is a measure you can name, while a ridge over four hundred collinear measures is what components are for. The loadings land beside these results and named after them, because a component means nothing without the table saying what loads on it.	"off"
pca_components	int	How many components to keep, when reducing. 0 decides by the Kaiser criterion -- a starting point, not an answer.	0
pca_retain	('parallel', 'kaiser')	How the component count is chosen when pca_components is 0: parallel analysis keeps a component while its eigenvalue beats what random data of the same size produce at that rank; the Kaiser rule keeps every eigenvalue above 1, which on a wide table is most of them.	"parallel"
pca_rotation	bool	Rotate the components (varimax) so each loads on a small cluster of features and is therefore nameable.	True
split_col	Optional[str]	Fit separately for each value of this column, labeling each result with it. Set to "source_col" when a spreadsheet's text columns were measured separately: a participant then has one row per column, and cross-validating over them together leaks -- the same person appears in the training and the held-out fold, and the reported performance is not out-of-sample at all.	None
alphas	Optional[Sequence[float]]	Penalties to search. Default: 33 points from 1e-3 to 1e5. Every value must be greater than zero; alpha = 0 is ordinary least squares, which is exactly what ridge exists to avoid on collinear features.	None
n_folds	int	Cross-validation folds. The reported performance is out-of-fold: each row is predicted by a model that never saw it.	5
stratify	bool	Make the folds alike: each gets the same spread of the outcome and sizes within one of each other (rows are ordered by the outcome and dealt round the folds). Off, the folds are a plain seeded shuffle, which can put most of the high scores in one fold.	True
zscore	bool	Standardize predictors on the training statistics. Coefficients are then per standard deviation and comparable to each other; with False they stay in the features' own units and are not.	True
max_missing	float	Set aside any predictor missing for more than this fraction of the rows, rather than letting it delete them. Some cohesion measures compare paragraphs two apart and so are blank for any text with fewer than three; left in, 60 such columns once deleted 885 of 938 rows and the model was fitted on 52. Set it to 1.0 to keep every column and accept the row loss; what was dropped is always reported.	MAX_MISSING
seed	int	Fixes the fold shuffle. Stored in the model, because a cross-validated number nobody can reproduce is a number nobody can check.	0
out_dir	Optional[PathLike]	Where the tables and the model files go; both default beside the analysis table (models under models/).	None
out_models_dir	Optional[PathLike]	Where the tables and the model files go; both default beside the analysis table (models under models/).	None
overwrite_existing	bool	When False (default) and the metrics table exists, it is returned untouched.	False
rounding	int	Decimal places in the output tables.	4

Returns:

Type	Description
Path	ridge_cv_metrics.csv -- one row per (feature set, outcome), which is also the table to sort when comparing feature sets.

Notes

Rows missing any predictor or the outcome sit out that model's fit entirely (listwise), and n_used in the metrics table says how many were left. A predictor that never varies in the training rows is dropped by name rather than silently contributing nothing.

Source code in src\taters\stats\ridge.py

def fit_ridge_csv(
    *,
    table_csv: PathLike,
    outcome_cols: Sequence[str],
    feature_sets=None,
    control_cols: Sequence[str] = (),
    categorical_controls: Sequence[str] = (),
    control_combos: str = "none_and_all",
    set_combos: Literal["none", "each_and_all", "subsets"] = "subsets",
    pca: object = "off",
    pca_components: int = 0,
    pca_retain: Literal["parallel", "kaiser"] = "parallel",
    pca_rotation: bool = True,
    split_col: Optional[str] = None,
    alphas: Optional[Sequence[float]] = None,
    n_folds: int = 5,
    stratify: bool = True,
    zscore: bool = True,
    max_missing: float = MAX_MISSING,
    seed: int = 0,
    out_dir: Optional[PathLike] = None,
    out_models_dir: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[..., None]] = None,
    verbose: bool = True,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
) -> Path:
    """
    Fit a cross-validated ridge model per outcome, and keep it.

    Parameters
    ----------
    table_csv
        The assembled analysis table.
    outcome_cols
        The numeric column(s) to predict. Each gets its own model, its own
        chosen penalty and its own row of performance.
    feature_sets
        ``None`` to use every feature as one predictor set, ``"per_table"``
        to fit a separate model per source feature table -- which is how you
        answer "which feature set actually predicts this?", since the
        resulting rows sit side by side in one metrics table -- or an
        explicit ``{name: [columns]}``.
    control_cols
        Columns to fit alongside -- and instead of -- the language, so the
        language's own contribution can be read off. With controls the run
        fits three models per outcome on **one common sample**: the controls
        alone, the language alone, and both; the combined model's row
        carries ``delta_r2_over_controls``, which is what "language adds
        this much over age and gender" means. Continuous controls enter as
        themselves, categorical ones as indicators against a reference
        level.

        The shared sample is the part that matters: fitting each model on
        whatever rows it happened to have would make the difference between
        them part sample and part language, with no way to tell which.
    categorical_controls
        Which of ``control_cols`` to treat as categories despite looking
        numeric.
    set_combos : {"none", "each_and_all", "subsets"}, default="subsets"
        When the feature tables are analyzed together (``feature_sets`` not
        given), how to combine them: ``subsets`` fits every combination of
        the tables -- each alone, every pair, and so on up to all of them --
        so the report can say what each table adds to the others;
        ``each_and_all`` fits each alone and all together; ``none`` fits
        only all together. Every combination is fitted for up to five
        tables; above that, each and all. A combination is named by its
        members (``dictionary+readability``) in the tables and the models.
    control_combos : {"none_and_all", "each", "subsets"}, default="none_and_all"
        How many control sets to try. The default pair answers the usual
        question. ``"each"`` adds one model per single control, which is how
        you see which of them is doing the work; ``"subsets"`` runs all
        2^k of them, which is thorough and grows the table accordingly.
    pca : str or list of str, default="off"
        Analyze components instead of the raw measures. ``"off"`` uses the
        features as they are; ``"all"`` reduces every feature set; a list of
        set names reduces those and leaves the rest alone -- which is the
        case worth having, since a hundred dictionary categories are worth
        reducing and eight readability indices are not.

        Set per analysis, deliberately: raw variables read better in a
        correlation table, where each row is a measure you can name, while a
        ridge over four hundred collinear measures is what components are
        for. The loadings land beside these results and named after them,
        because a component means nothing without the table saying what
        loads on it.
    pca_components : int, default=0
        How many components to keep, when reducing. ``0`` decides by the
        Kaiser criterion -- a starting point, not an answer.
    pca_retain : {"parallel", "kaiser"}, default="parallel"
        How the component count is chosen when ``pca_components`` is 0:
        parallel analysis keeps a component while its eigenvalue beats what
        random data of the same size produce at that rank; the Kaiser rule
        keeps every eigenvalue above 1, which on a wide table is most of them.
    pca_rotation : bool, default=True
        Rotate the components (varimax) so each loads on a small cluster of
        features and is therefore nameable.
    split_col
        Fit separately for each value of this column, labeling each result
        with it. Set to ``"source_col"`` when a spreadsheet's text columns
        were measured separately: a participant then has one row per column,
        and cross-validating over them together leaks -- the same person
        appears in the training and the held-out fold, and the reported
        performance is not out-of-sample at all.
    alphas
        Penalties to search. Default: 33 points from 1e-3 to 1e5. Every value
        must be greater than zero; alpha = 0 is ordinary least squares, which
        is exactly what ridge exists to avoid on collinear features.
    n_folds
        Cross-validation folds. The reported performance is out-of-fold:
        each row is predicted by a model that never saw it.
    stratify
        Make the folds alike: each gets the same spread of the outcome and
        sizes within one of each other (rows are ordered by the outcome and
        dealt round the folds). Off, the folds are a plain seeded shuffle,
        which can put most of the high scores in one fold.
    zscore
        Standardize predictors on the training statistics. Coefficients are
        then per standard deviation and comparable to each other; with False
        they stay in the features' own units and are not.
    max_missing
        Set aside any predictor missing for more than this fraction of the
        rows, rather than letting it delete them. Some cohesion measures
        compare paragraphs two apart and so are blank for any text with
        fewer than three; left in, 60 such columns once deleted 885 of 938
        rows and the model was fitted on 52. Set it to 1.0 to keep every
        column and accept the row loss; what was dropped is always reported.
    seed
        Fixes the fold shuffle. Stored in the model, because a
        cross-validated number nobody can reproduce is a number nobody can
        check.
    out_dir, out_models_dir
        Where the tables and the model files go; both default beside the
        analysis table (models under ``models/``).
    overwrite_existing
        When False (default) and the metrics table exists, it is returned
        untouched.
    rounding
        Decimal places in the output tables.

    Returns
    -------
    Path
        ``ridge_cv_metrics.csv`` -- one row per (feature set, outcome), which
        is also the table to sort when comparing feature sets.

    Notes
    -----
    Rows missing any predictor or the outcome sit out that model's fit
    entirely (listwise), and ``n_used`` in the metrics table says how many
    were left. A predictor that never varies in the training rows is dropped
    by name rather than silently contributing nothing.
    """
    import numpy as np

    table_csv, folder, models_dir = output_folders(table_csv, out_dir,
                                                   out_models_dir)
    metrics_path = folder / "ridge_cv_metrics.csv"
    coef_path = folder / "ridge_coefficients.csv"
    path_path = folder / "ridge_alpha_path.csv"
    if reusable(metrics_path, table_csv, overwrite_existing=overwrite_existing,
                verbose=verbose, what="the ridge results"):
        if verbose:
            print(f"Ridge results already exist; returning existing file: "
                  f"{metrics_path}")
        return metrics_path

    outcome_cols = [str(c) for c in outcome_cols]
    inputs = prepare_fit(
        table_csv=table_csv, folder=folder, stem="ridge",
        outcome_cols=outcome_cols, feature_sets=feature_sets,
        control_cols=control_cols, categorical_controls=categorical_controls,
        control_combos=control_combos, pca=pca, pca_components=pca_components,
        pca_rotation=pca_rotation, pca_retain=pca_retain, split_col=split_col,
        set_combos=set_combos, alphas=alphas,
        n_folds=n_folds, max_missing=max_missing, encoding=encoding,
        rounding=rounding, verbose=verbose, on_progress=on_progress,
        zero_alpha_reason=("alpha = 0 is ordinary least squares, which is "
                           "what ridge exists to avoid here"))
    n_folds = int(n_folds)
    table, grid, sets, columns = (inputs.table, inputs.grid, inputs.sets,
                                  inputs.columns)
    pca_notes = inputs.pca_notes
    control_matrix, control_names, control_notes, control_spec = (
        inputs.control_matrix, inputs.control_names, inputs.control_notes,
        inputs.control_spec)
    control_sets, subsets, lead, total = (inputs.control_sets, inputs.subsets,
                                          inputs.lead, inputs.total)
    for c in outcome_cols:
        _refuse_categorical_outcome(c, table[c].tolist())
    outcomes = {c: numeric_column(table[c].tolist(), column=c)
                for c in outcome_cols}
    # (outcome, control-set, with-language) -> pooled R2, so that we can say
    # how much the combined model added over the controls alone.
    model_scores: dict = {}

    metrics_rows, coef_rows, path_rows, fold_rows = [], [], [], []
    best: dict = {}
    gaps: list = []
    done = 0

    for subset, subset_mask in subsets:
      prefix = [subset] if split_col else []
      for set_name, all_cols in sets.items():
        cols, sparse, design_all, design = set_design(
            set_name, all_cols, columns, subset_mask, max_missing,
            verbose=verbose, tag="ridge")
        gaps += gap_lines(set_name, cols, design)
        model_doc = model_header(
            kind="taters-ridge-model", fmt=RIDGE_MODEL_FORMAT,
            set_name=set_name, cols=cols, table_csv=table_csv,
            reduction=inputs.reductions.get(set_name),
            cv={"n_folds": n_folds, "seed": seed, "stratify": bool(stratify),
                "selection": "min_mean_out_of_fold_rmse_prefer_larger_alpha",
                "alphas": grid},
            zscore=zscore)
        predictions = {}
        # what each outcome's own model file needs beyond the shared header:
        # its coefficient block, and the predictors and controls it was fitted
        # with (they can differ between outcomes when a constant column gets
        # dropped for one and not another)
        per_outcome: dict = {}

        for outcome in outcome_cols:
            y_all = (outcomes[outcome] if subset_mask is None
                     else outcomes[outcome][subset_mask])
            ctrl = (control_matrix if subset_mask is None
                    else control_matrix[subset_mask])

            # ONE complete-case sample for every model of this outcome. if we
            # did listwise deletion over different column sets, the
            # controls-only model would get more rows than the combined one,
            # and the difference between their scores (the whole point of
            # running both) would be part sample and part language. no good.
            complete = ~np.isnan(design).any(axis=1) & ~np.isnan(y_all)
            if ctrl.shape[1]:
                complete = complete & ~np.isnan(ctrl).any(axis=1)
            y = y_all[complete]
            n_used = int(complete.sum())
            # how many rows COULD have been used, i.e. the outcome is present.
            # we print this next to n_used because the gap between the two is
            # a question anyone reading these results needs to be able to ask.
            n_available = int((~np.isnan(y_all)).sum())
            floor = max(_MIN_ROWS_PER_FOLD * n_folds, _MIN_ROWS)
            if n_used < floor:
                raise ValueError(
                    f"{outcome} on feature set {set_name!r}"
                    + (f" ({split_col} {subset!r})" if subset else "")
                    + f": only {n_used} row(s) have every predictor"
                    + (", every control" if ctrl.shape[1] else "")
                    + f" and the outcome, and {n_folds}-fold "
                    f"cross-validation needs at least {floor}. Use fewer "
                    f"features, fill the gaps, or fewer folds.")

            keep_best = None
            for control_set in control_sets:
                idx = [control_names.index(c) for c in control_names
                       if c.split("=")[0] in control_set]
                for with_language in (False, True):
                    if not with_language and not idx:
                        # "neither the controls nor the language" isn't a
                        # model, it's just the mean. skip it.
                        continue
                    done += 1
                    if on_progress is not None:
                        on_progress(done, total,
                                    f"fitting {outcome} on {set_name}")
                    blocks, names = [], []
                    if idx:
                        blocks.append(ctrl[complete][:, idx])
                        names += [control_names[i] for i in idx]
                    if with_language:
                        blocks.append(design[complete])
                        names += list(cols)
                    x = np.hstack(blocks)
                    label = ("controls+language" if idx and with_language
                             else "controls" if idx else "language")

                    fit, dropped = _cv_ridge(
                        x, y, names=names, grid=grid, n_folds=n_folds,
                        seed=seed, zscore=zscore, stratify=stratify)
                    if fit is None:
                        raise ValueError(
                            f"{outcome} on feature set {set_name!r}: every "
                            f"predictor is constant across the usable rows, "
                            f"so there is nothing to predict from.")
                    if dropped and verbose:
                        print(f"[ridge] {set_name}/{outcome} ({label}): "
                              f"dropping constant predictor(s) {dropped}")

                    cv = fit["cv"]
                    combo = "|".join(control_set)
                    # language's own contribution: same rows, same folds,
                    # same penalty grid -- we hold everything equal except
                    # whether the words were in the model.
                    delta = ""
                    if idx and with_language:
                        base = model_scores.get((outcome, combo, False))
                        if base is not None and base == base:
                            delta = fmt(cv["r2"] - base, rounding)
                    model_scores[(outcome, combo, with_language)] = cv["r2"]

                    metrics_rows.append(prefix + [
                        set_name, str(len(inputs.set_members.get(set_name, [set_name]))),
                        str(len(all_cols)), outcome, label, combo,
                        str(n_used), str(n_available),
                        str(len(sparse)), str(len(fit["kept"])),
                        str(len(dropped)), str(n_folds), str(bool(zscore)),
                        fmt(fit["alpha"], 8),
                        fmt(cv["r2"], rounding), delta,
                        fmt(cv["r2_folds"], rounding),
                        fmt(cv["r2_folds_se"], rounding),
                        fmt(cv["r"], rounding), fmt(cv["r_p"], rounding),
                        fmt(cv["r_folds"], rounding),
                        fmt(cv["rho"], rounding), fmt(cv["rho_p"], rounding),
                        fmt(cv["mse"], rounding), fmt(cv["rmse"], rounding),
                        fmt(cv["mae"], rounding),
                        fmt(cv["baseline_mae"], rounding),
                        fmt(fit["train_r2"], rounding),
                        fmt(fit["intercept"], rounding)])
                    for a, rmse in zip(grid, fit["rmse_by_alpha"],
                                       strict=True):
                        path_rows.append(prefix + [set_name, outcome, label,
                                                   fmt(a, 8),
                                                   fmt(float(rmse), rounding)])
                    for f, m in enumerate(fit["per_fold"]):
                        fold_rows.append(prefix + [
                            set_name, outcome, label, str(f + 1),
                            str(int((fit["folds"] == f).sum())),
                            fmt(m["r2"], rounding), fmt(m["r"], rounding),
                            fmt(m["rho"], rounding), fmt(m["mae"], rounding)])
                    coef = np.asarray(fit["coef"])
                    for i in np.argsort(-np.abs(coef)):
                        coef_rows.append(
                            (tuple(prefix) + (set_name, label, combo,
                                              names[fit["kept"][i]]),
                             outcome, float(coef[i])))

                    if cv["r2"] == cv["r2"] and cv["r2"] > best.get(
                            outcome, (-1e9,))[0]:
                        best[outcome] = (
                            cv["r2"],
                            f"{set_name} @ {subset}" if subset else set_name,
                            fit["alpha"])
                    # the model we keep (and apply later on) is the richest
                    # one: everything the run had to predict with.
                    if with_language:
                        keep_best = (fit, names, idx)

            fit, names, control_idx = keep_best
            # the recipe for the control columns, not just their names. the
            # next table this model scores will have `gender` as whatever
            # words the researcher typed, and which level is the reference
            # lives only here.
            per_outcome[outcome] = {
                "controls": [control_spec[i] for i in control_idx],
                "predictors": names,
                "block": {
                    "n_rows": n_used,
                    "alpha": fit["alpha"],
                    "kept": fit["kept"],
                    "mu": fit["mu"],
                    "sigma": fit["sigma"],
                    "coef": fit["coef"],
                    "intercept": fit["intercept"],
                    "cv": fit["cv"],
                },
            }
            # we score through the same helper that apply uses, over every
            # row of the table -- so that a fit and an apply on the same data
            # agree to the byte. that's what keeps the two from drifting
            # apart.
            full = np.hstack([ctrl[:, [control_names.index(c) for c in names
                                       if c in control_names]],
                              design]) if control_names else design
            predictions[outcome] = predict_rows(
                full, fit["kept"], np.asarray(fit["mu"]),
                np.asarray(fit["sigma"]), np.asarray(fit["coef"]),
                fit["intercept"], rounding)


        model_doc["subset"] = subset
        models_dir.mkdir(parents=True, exist_ok=True)
        stem = slug(set_name) + (f"__{slug(subset)}" if subset else "")
        # one model file per outcome. we used to bundle every outcome of a
        # feature set into one file, and somebody who'd predicted five
        # personality traits got one model back and asked where the other
        # four were. a model is "the thing that predicts X", so each X gets
        # its own file, its own name in the library, and its own import
        for outcome, parts in per_outcome.items():
            doc = dict(model_doc)
            doc["outcomes"] = {outcome: parts["block"]}
            doc["controls"] = parts["controls"]
            doc["predictors"] = parts["predictors"]
            model_path = models_dir / f"ridge__{stem}__{slug(outcome)}.json"
            with atomic_write(model_path, mode="w", encoding="utf-8") as fh:
                json.dump(doc, fh, indent=1)
        _write_predictions(folder / f"ridge_predictions__{stem}.csv",
                           table if subset_mask is None else table[subset_mask],
                           outcome_cols, predictions, encoding)

    announce(on_progress, "writing the results")
    metrics_header = (
           lead + ["feature_set",
                   # how many tables the set is made of, and how many
                   # feature columns they brought along (before the sparse
                   # ones were set aside; `n_predictors` is what the model
                   # actually used, controls included) -- so that a row's
                   # score can be read against how much went into it.
                   "n_feature_sets", "n_features", "outcome",
                   # which model this row is: the controls alone, the
                   # language alone, or both -- and which controls, exactly.
                   "model", "controls",
                   "n_used", "n_available", "n_dropped_sparse",
                   "n_predictors",
                   "n_dropped_constant", "n_folds", "zscore", "alpha",
                   # variance explained, pooled and per-fold (the mean and
                   # the standard error of that mean), plus what the
                   # language added over the controls alone.
                   "cv_r2", "delta_r2_over_controls",
                   "cv_r2_folds", "cv_r2_folds_se",
                   # predicted-vs-observed association, three ways.
                   "cv_r", "cv_r_p", "cv_r_folds", "cv_rho", "cv_rho_p",
                   # error, and the floor we should read it against.
                   "cv_mse", "cv_rmse", "cv_mae", "baseline_mae",
                   "train_r2", "intercept"])
    _write_csv(metrics_path, metrics_header, metrics_rows, encoding)
    _write_coefficients(coef_path, coef_rows, outcome_cols, lead, rounding,
                        encoding)
    _write_csv(path_path,
           lead + ["feature_set", "outcome", "model", "alpha", "cv_rmse"],
           path_rows, encoding)
    _write_csv(folder / "ridge_folds.csv",
           lead + ["feature_set", "outcome", "model", "fold", "n", "r2", "r",
                   "rho", "mae"], fold_rows, encoding)

    section = _section_md(
        sets=sets, outcome_cols=outcome_cols, best=best, n_folds=n_folds,
        stratify=stratify,
        zscore=zscore, metrics_rows=metrics_rows,
        metrics_header=metrics_header, control_names=control_names,
        control_notes=control_notes, model_scores=model_scores,
        gap_notes=gaps_section(gaps, max_missing=max_missing),
        set_notes=inputs.set_notes)
    section += components_appendix(pca_notes)
    write_section(folder, "ridge", section)
    if verbose:
        print(f"[ridge] {len(metrics_rows)} model(s) -> {metrics_path}")
    return metrics_path

predict_rows

predict_rows(
    cells, kept, mu, sigma, coef, intercept, rounding
)

Predictions for a block of rows -- the one scorer fit and apply share.

cells is (n_rows, n_predictors) in the model's predictor order. A row with a missing value among the predictors it needs cannot be scored: the missing value carries through the arithmetic as NaN -- deliberately, not incidentally -- and the row comes back as None for the caller to write as a blank cell. A prediction built from a filled-in guess would look exactly like a real one in the output.

Only the predictors the model actually kept are read, so a row missing a predictor that was dropped as constant still scores.

Source code in src\taters\stats\ridge.py

def predict_rows(cells, kept, mu, sigma, coef, intercept, rounding: int):
    """
    Predictions for a block of rows -- the one scorer fit and apply share.

    ``cells`` is (n_rows, n_predictors) in the model's predictor order. A row
    with a missing value among the predictors it needs cannot be scored: the
    missing value carries through the arithmetic as NaN -- deliberately, not
    incidentally -- and the row comes back as ``None`` for the caller to
    write as a blank cell. A prediction built from a filled-in guess would
    look exactly like a real one in the output.

    Only the predictors the model actually kept are read, so a row missing a
    predictor that was dropped as constant still scores.
    """
    import numpy as np

    x = np.asarray(cells, dtype=float)[:, kept]
    predicted = ((x - mu) / sigma) @ coef + intercept
    return [None if np.isnan(v) else round(float(v), rounding)
            for v in predicted]

taters.stats.report

The plain-English report: one markdown file summarizing a run's statistics.

Each analysis writes its own section fragment into stats_results/_sections/ as it runs (see :func:taters.stats._common.write_section); this module just concatenates them in their fixed order under a dated header. It computes no statistics of its own -- every number in the report was produced, and is owned, by the module that understands it. That division is what lets a future analysis join the report by writing one fragment, with no edits here.

write_stats_report

write_stats_report(
    *,
    stats_dir="stats_results",
    title=None,
    out_md=None,
    keep_table=True,
    overwrite_existing=False,
    verbose=True
)

Concatenate the analyses' section fragments into report.md.

Parameters:

Name	Type	Description	Default
stats_dir	PathLike	The stats output folder holding _sections/.	'stats_results'
title	Optional[str]	Heading for the report; defaults to "Statistical results".	None
out_md	Optional[PathLike]	Where to write; defaults to <stats_dir>/report.md.	None
keep_table	bool	Whether to leave the merged analysis table in place. Deleting it is this step's job rather than the assemble step's for the obvious reason: every analysis reads that table, so the only moment it is safe to remove is after the last one has run, and this is the step that runs last. The sidecar describing it goes with it -- a map of a file that is not there helps nobody -- while the row accounting in assemble_manifest.json stays, because how many rows were analyzed is part of the result.	True
overwrite_existing	bool	Unlike the analyses, the report is cheap and derived, so the default skip is mostly about symmetry: pass True to rebuild after re-running any analysis.	False

Returns:

Type	Description
Path	The report file.

Source code in src\taters\stats\report.py

def write_stats_report(
    *,
    stats_dir: PathLike = "stats_results",
    title: Optional[str] = None,
    out_md: Optional[PathLike] = None,
    keep_table: bool = True,
    overwrite_existing: bool = False,
    verbose: bool = True,
) -> Path:
    """
    Concatenate the analyses' section fragments into ``report.md``.

    Parameters
    ----------
    stats_dir
        The stats output folder holding ``_sections/``.
    title
        Heading for the report; defaults to "Statistical results".
    out_md
        Where to write; defaults to ``<stats_dir>/report.md``.
    keep_table
        Whether to leave the merged analysis table in place. Deleting it is
        this step's job rather than the assemble step's for the obvious
        reason: every analysis reads that table, so the only moment it is
        safe to remove is after the last one has run, and this is the step
        that runs last. The sidecar describing it goes with it -- a map of a
        file that is not there helps nobody -- while the row accounting in
        ``assemble_manifest.json`` stays, because how many rows were
        analyzed is part of the result.
    overwrite_existing
        Unlike the analyses, the report is cheap and derived, so the default
        skip is mostly about symmetry: pass True to rebuild after re-running
        any analysis.

    Returns
    -------
    Path
        The report file.
    """
    import datetime

    stats_dir = Path(stats_dir)
    out_path = Path(out_md) if out_md else stats_dir / "report.md"
    sections_dir = stats_dir / SECTIONS_DIR
    from ._common import reusable

    if reusable(out_path, *(sorted(sections_dir.glob("*.md"))
                            if sections_dir.is_dir() else []),
                overwrite_existing=overwrite_existing, verbose=verbose,
                what="the report"):
        if verbose:
            print(f"Report already exists; returning existing file: {out_path}")
        return out_path
    fragments = sorted(sections_dir.glob("*.md")) if sections_dir.is_dir() else []
    if not fragments:
        raise FileNotFoundError(
            f"no report sections found under {sections_dir}. Run at least "
            f"one analysis first -- each one writes its own section.")

    version = taters_version()
    header = [
        f"# {title or 'Statistical results'}",
        "",
        f"*Written {datetime.date.today().isoformat()}"
        + (f" by taters {version}" if version else " by taters")
        + f"; tables live beside this file in `{stats_dir.name}/`.*",
        "",
    ]
    body = []
    for fragment in fragments:
        body.append(fragment.read_text(encoding="utf-8").rstrip())
        body.append("")

    out_path.parent.mkdir(parents=True, exist_ok=True)
    with atomic_write(out_path, mode="w", encoding="utf-8") as fh:
        fh.write("\n".join(header + body).rstrip() + "\n")
    if not keep_table:
        for name in ("analysis_table.csv", "analysis_table_sets.json"):
            target = stats_dir / name
            if target.is_file():
                target.unlink()
                if verbose:
                    print(f"[report] removed {name} (keep_table is off)")
    if verbose:
        print(f"[report] {len(fragments)} section(s) -> {out_path}")
    return out_path

taters.stats.pca

Streaming, exact PCA with varimax rotation, over CSVs of any length.

The engine behind the MEM topic model, and a standalone tool for reducing any feature table (dictionary scores, readability measures, embeddings) before a downstream model. Two design facts:

• Exact, one pass, bounded memory. Feature tables are long (rows) but narrow (columns), so the fit streams the rows once, accumulating the column sums and the p-by-p cross-product in chunks (BLAS does the work: ~50x faster than row-at-a-time), then eigendecomposes the correlation matrix in memory. No randomized SVD, no seed, no dask: the decomposition is verified against scikit-learn's full SVD to machine precision, and the finished rotated solution against R's psych::principal to ~2 decimal places on a real 1,000-term MEM (see varimax -- Kaiser normalization plus iterate-to-convergence, deliberately correcting the dask script this module replaces, whose rotation was unnormalized and capped at 20 sweeps). Randomized methods earn their keep when p is tens of thousands; if that day comes, it is a solver option behind this same interface.

• A fit is a reusable instrument. fit_pca_csv writes a model file with the feature names, the training means and deviations, and the projection; apply_pca_csv matches features by name in any new table, standardizes with the TRAINING statistics, and projects. Applying a model to its own training table reproduces the training scores exactly.

Eigenvalues are those of the correlation matrix (they average 1.0), so both retention rules mean what the textbooks say: parallel analysis (the default for n_components=0) keeps a component while its eigenvalue beats what random data of the same size produce at that rank, and the Kaiser criterion keeps every eigenvalue >= 1. (dask-ml and scikit-learn scale variances by n-1 rather than n; that uniform sqrt((n-1)/n) on loadings is a bookkeeping convention, not a disagreement.)

apply_in_memory

apply_in_memory(x, axes)

Project rows onto axes from :func:fit_in_memory.

Standardized with the axes' own means and deviations, not ones recomputed from the rows being projected -- which is what lets a saved reduction score a new table on the original sample's scale.

Source code in src\taters\stats\pca.py

def apply_in_memory(x, axes: dict):
    """Project rows onto axes from :func:`fit_in_memory`.

    Standardized with the axes' **own** means and deviations, not ones
    recomputed from the rows being projected -- which is what lets a saved
    reduction score a new table on the original sample's scale.
    """
    import numpy as np

    a = np.asarray(x, dtype=np.float64)[:, axes["kept"]]
    mu = np.asarray(axes["mu"], dtype=np.float64)
    sigma = np.asarray(axes["sigma"], dtype=np.float64)
    return ((a - mu) / sigma) @ np.asarray(axes["projection"],
                                           dtype=np.float64)

apply_pca_csv

apply_pca_csv(
    model_json,
    input_csv,
    *,
    out_scores_csv=None,
    overwrite_existing=False,
    on_progress=None,
    encoding="utf-8-sig",
    rounding=4
)

Score a new feature table on a saved PCA model.

The model's features are matched by name in the new table's header (any column order; extra columns pass through as identifiers), then each row is standardized with the training means and deviations and projected. A feature the model needs that the table lacks is refused -- filling it with anything would silently move every score.

Parameters:

Name	Type	Description	Default
model_json	PathLike	A *_pca_model.json written by :func:fit_pca_csv.	required
input_csv	PathLike	Any CSV whose header contains every feature the model was fit on.	required
out_scores_csv	optional	Defaults to <input>_pca_scores.csv beside the input.	None
overwrite_existing	bool	If False and the output already exists, skip and return it.	False
encoding	str	Encoding for reading and writing CSV files.	"utf-8-sig"
rounding	int	Decimal places for the component scores.	4

Returns:

Type	Description
Path	out_scores_csv: the input's non-feature columns, then the model's components.

Source code in src\taters\stats\pca.py

def apply_pca_csv(
    model_json: PathLike,
    input_csv: PathLike,
    *,
    out_scores_csv: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[[int, int], None]] = None,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
) -> Path:
    """
    Score a new feature table on a saved PCA model.

    The model's features are matched **by name** in the new table's header
    (any column order; extra columns pass through as identifiers), then each
    row is standardized with the *training* means and deviations and
    projected. A feature the model needs that the table lacks is refused --
    filling it with anything would silently move every score.

    Parameters
    ----------
    model_json
        A ``*_pca_model.json`` written by :func:`fit_pca_csv`.
    input_csv
        Any CSV whose header contains every feature the model was fit on.
    out_scores_csv : optional
        Defaults to ``<input>_pca_scores.csv`` beside the input.
    overwrite_existing : bool, default=False
        If ``False`` and the output already exists, skip and return it.
    encoding : str, default="utf-8-sig"
        Encoding for reading and writing CSV files.
    rounding : int, default=4
        Decimal places for the component scores.

    Returns
    -------
    Path
        ``out_scores_csv``: the input's non-feature columns, then
        the model's components.
    """
    import numpy as np

    model = load_pca_model(model_json)
    features: List[str] = list(model["features"])
    fit = model["model"]
    kept = np.asarray(fit["kept"], dtype=int)
    mu = np.asarray(fit["mu"], dtype=np.float64)
    sigma = np.asarray(fit["sigma"], dtype=np.float64)
    projection = np.asarray(fit["projection"], dtype=np.float64)
    comp_names = list(fit["components"])

    input_csv = Path(input_csv)
    if not input_csv.exists():
        raise FileNotFoundError(f"input_csv not found: {input_csv}")
    scores_path = Path(out_scores_csv) if out_scores_csv else \
        input_csv.with_name(f"{input_csv.stem}_pca_scores.csv")
    scores_path.parent.mkdir(parents=True, exist_ok=True)
    from ._common import reusable

    if reusable(scores_path, input_csv, model_json,
                overwrite_existing=overwrite_existing, verbose=True,
                what="the PCA scores"):
        print("PCA scores output file already exists; returning existing file.")
        return scores_path

    header = _read_header(input_csv, encoding)
    positions = {name: i for i, name in enumerate(header)}
    missing = [f for f in features if f not in positions]
    if missing:
        raise ValueError(
            f"{input_csv} is missing {len(missing)} feature(s) the model was "
            f"fit on ({', '.join(missing[:6])}{'…' if len(missing) > 6 else ''}). "
            "The new table must carry every original feature, by name."
        )
    feature_pos = [positions[f] for f in features]
    id_pos = [i for i, name in enumerate(header) if name not in set(features)]

    ticker = Ticker(on_progress, count_rows(input_csv, on_progress=on_progress))
    with input_csv.open("r", newline="", encoding=encoding) as f, \
            atomic_write(scores_path, newline="", encoding=encoding) as out:
        reader = csv.reader(f)
        next(reader)
        writer = csv.writer(out)
        writer.writerow([*(header[i] for i in id_pos), *comp_names])
        for row in reader:
            ticker.tick(message="scoring components")
            cells = [row[i] for i in feature_pos]
            writer.writerow([*(row[i] for i in id_pos),
                             *project_row(cells, kept, mu, sigma, projection,
                                          rounding)])
    return scores_path

component_names

component_names(n)

Component_1 … Component_n -- one name everywhere. The file API wrote PC_1 and the in-analysis reduction Component_1 for the same thing, so a reader met two names for one idea in one results folder. Component_1 … Component_n -- what the analyses call them, so a loadings table and a results table can be read side by side.

Source code in src\taters\stats\pca.py

def component_names(n: int) -> list:
    """``Component_1 … Component_n`` -- one name everywhere. The file API
    wrote ``PC_1`` and the in-analysis reduction ``Component_1`` for the same
    thing, so a reader met two names for one idea in one results folder.
    ``Component_1 … Component_n`` -- what the analyses call them, so a
    loadings table and a results table can be read side by side."""
    return [f"Component_{i + 1}" for i in range(n)]

describe_retention

describe_retention(retention)

One sentence on how the component count was chosen, for a report.

Source code in src\taters\stats\pca.py

def describe_retention(retention: dict) -> str:
    """One sentence on how the component count was chosen, for a report."""
    rule = retention.get("rule")
    if rule == "parallel":
        eig = retention.get("unrotated_eigenvalues") or []
        thr = retention.get("thresholds") or []
        k = max(1, len(eig) - 1) if len(eig) > 1 else len(eig)
        tail = ""
        if len(eig) > k:
            tail = (f"; the next eigenvalue, {eig[k]:.2f}, fell short of its "
                    f"chance level of {thr[k]:.2f}")
        return (f"chosen by parallel analysis: a component is kept while its "
                f"eigenvalue beats the {retention.get('percentile', 95):g}th "
                f"percentile of {retention.get('draws', PARALLEL_DRAWS)} random "
                f"data sets of the same size{tail}")
    if rule == "kaiser":
        return "chosen by the Kaiser rule: every component with an eigenvalue above 1"
    return "the number asked for"

fit_axes

fit_axes(
    n,
    sums,
    cross,
    *,
    n_components,
    rotation,
    retain="parallel",
    on_progress=None
)

From streamed moments to finished axes.

Returns (kept, mu, sigma, loadings, projection, eigenvalues, pct, retention): kept indexes the non-constant columns; loadings is the rotated feature-by-component matrix (interpretation); projection scores a row -- standardized kept cells @ projection -- and both share one component order (descending rotated eigenvalue) and one sign convention (each component's strongest feature loads positive), so reruns cannot come back mirror-flipped or shuffled. retention says how the count was chosen: the rule, and for parallel analysis the unrotated eigenvalues beside the chance thresholds up to the first one that fell short, so the decision can be read and argued with.

n_components > 0 is honored as asked; 0 chooses by retain, one of :data:RETAIN_RULES.

Source code in src\taters\stats\pca.py

def fit_axes(n: int, sums, cross, *, n_components: int, rotation: bool,
             retain: str = "parallel", on_progress=None):
    """
    From streamed moments to finished axes.

    Returns ``(kept, mu, sigma, loadings, projection, eigenvalues, pct,
    retention)``: ``kept`` indexes the non-constant columns; ``loadings`` is
    the rotated feature-by-component matrix (interpretation); ``projection``
    scores a row -- standardized kept cells @ projection -- and both share
    one component order (descending rotated eigenvalue) and one sign
    convention (each component's strongest feature loads positive), so
    reruns cannot come back mirror-flipped or shuffled. ``retention`` says
    how the count was chosen: the rule, and for parallel analysis the
    unrotated eigenvalues beside the chance thresholds up to the first one
    that fell short, so the decision can be read and argued with.

    ``n_components > 0`` is honored as asked; ``0`` chooses by ``retain``,
    one of :data:`RETAIN_RULES`.
    """
    import numpy as np

    mu_all = sums / n
    var_all = np.maximum(cross.diagonal() / n - mu_all ** 2, 0.0)
    kept = np.flatnonzero(var_all > 0)
    if kept.size == 0:
        raise ValueError("Every feature column is constant; there is nothing "
                         "for a PCA to explain.")
    mu = mu_all[kept]
    sigma = np.sqrt(var_all[kept])

    cov = cross[np.ix_(kept, kept)] / n - np.outer(mu, mu)
    corr = cov / np.outer(sigma, sigma)
    corr = np.clip((corr + corr.T) / 2, -1.0, 1.0)      # force it to be symmetric

    eigvals, eigvecs = np.linalg.eigh(corr)             # these come out ascending, so flip
    eigvals = np.maximum(eigvals[::-1], 0.0)
    eigvecs = eigvecs[:, ::-1]

    retention: dict = {"rule": "asked", "n_components": int(n_components)}
    if n_components > 0:
        k = min(int(n_components), kept.size)
    elif retain == "parallel":
        thresholds = parallel_thresholds(n, kept.size, on_progress=on_progress)
        # we keep going while the eigenvalue beats chance at its rank, and
        # stop at the first one that doesn't (that's Horn's procedure). a
        # later one that happens to beat its threshold isn't a component
        # worth naming.
        k = 0
        while k < kept.size and eigvals[k] > thresholds[k]:
            k += 1
        k = max(1, k)
        shown = min(k + 1, kept.size)
        retention = {"rule": "parallel", "draws": PARALLEL_DRAWS,
                     "percentile": PARALLEL_PERCENTILE,
                     "unrotated_eigenvalues": [float(v) for v in eigvals[:shown]],
                     "thresholds": [float(v) for v in thresholds[:shown]]}
    elif retain == "kaiser":
        # Kaiser's rule: keep the components that explain more than one
        # feature's worth of variance.
        k = max(1, int(np.sum(eigvals >= 1.0)))
        retention = {"rule": "kaiser",
                     "unrotated_eigenvalues": [float(v) for v in eigvals[:k + 1]]}
    else:
        raise ValueError(f"retain must be one of {RETAIN_RULES}, got {retain!r}")

    unrotated = eigvecs[:, :k] * np.sqrt(eigvals[:k])   # the loadings, before rotation
    if rotation and k > 1:
        loadings, R = varimax(unrotated)
    else:
        loadings, R = unrotated, np.eye(k)

    comp_eig = (loadings ** 2).sum(axis=0)
    order = np.argsort(comp_eig)[::-1]
    loadings = loadings[:, order]
    comp_eig = comp_eig[order]
    flips = np.where(loadings[np.abs(loadings).argmax(axis=0),
                              np.arange(k)] < 0, -1.0, 1.0)
    loadings = loadings * flips

    # the scores have to go through the same rotation, order, and sign flips.
    projection = (eigvecs[:, :k] @ R)[:, order] * flips
    pct = comp_eig / kept.size * 100
    return kept, mu, sigma, loadings, projection, comp_eig, pct, retention

fit_in_memory

fit_in_memory(
    x, *, n_components=0, rotation=True, retain="parallel"
)

Fit the same axes as :func:fit_pca_csv, from a matrix already in memory.

The analyses hold their design matrix in memory already -- the whole filtered table, since a reduction is fitted once per analysis before the analysis runs (see the planning decisions) -- so they cannot go through the CSV-streaming path, which reads a file. Same moments, same :func:fit_axes, same sign and order conventions, so a reduction fitted here and one fitted from a file agree.

Returns the dict :func:apply_in_memory consumes, plus the loadings and eigenvalues a reader wants to inspect.

Source code in src\taters\stats\pca.py

def fit_in_memory(x, *, n_components: int = 0, rotation: bool = True,
                  retain: str = "parallel"):
    """
    Fit the same axes as :func:`fit_pca_csv`, from a matrix already in memory.

    The analyses hold their design matrix in memory already -- the whole
    filtered table, since a reduction is fitted once per analysis before the
    analysis runs (see the planning decisions) -- so they cannot go through
    the CSV-streaming path, which reads a file. Same moments, same
    :func:`fit_axes`, same sign and order conventions, so a reduction fitted
    here and one fitted from a file agree.

    Returns the dict :func:`apply_in_memory` consumes, plus the loadings and
    eigenvalues a reader wants to inspect.
    """
    import numpy as np

    a = np.asarray(x, dtype=np.float64)
    n = int(a.shape[0])
    if n < 3:
        raise ValueError(
            f"a PCA needs at least three rows to find a direction in, got {n}.")
    sums = a.sum(axis=0)
    cross = a.T @ a
    kept, mu, sigma, loadings, projection, eigenvalues, pct, retention = fit_axes(
        n, sums, cross, n_components=n_components, rotation=rotation,
        retain=retain)
    return {
        "retention": retention,
        "kept": [int(i) for i in kept],
        "mu": [float(v) for v in mu],
        "sigma": [float(v) for v in sigma],
        "projection": [[float(v) for v in row] for row in projection],
        "loadings": [[float(v) for v in row] for row in loadings],
        "eigenvalues": [float(v) for v in eigenvalues],
        "pct_variance": [float(v) for v in pct],
        "n_components": int(projection.shape[1]),
    }

fit_pca_csv

fit_pca_csv(
    input_csv,
    *,
    start_col=2,
    n_components=0,
    retain="parallel",
    rotation=True,
    out_scores_csv=None,
    out_model_json=None,
    out_loadings_csv=None,
    out_eigenvalues_csv=None,
    overwrite_existing=False,
    on_progress=None,
    encoding="utf-8-sig",
    rounding=4
)

Fit a varimax-rotated PCA to the numeric columns of any feature CSV.

Parameters:

Name	Type	Description	Default
input_csv	PathLike	Any wide feature table: identifier column(s) first, numeric feature columns after. Every Taters feature file (dictionary scores, readability, a document-term matrix) has this shape.	required
start_col	int	1-based position of the first feature column, exactly as in the house dask script this replaces. Everything before it is carried into the scores file unchanged as identifiers. The default fits tables shaped text_id, <features...>.	2
n_components	int	0 picks automatically by retain; an explicit number is honored up to the count of non-constant features.	0
retain	('parallel', 'kaiser')	How 0 chooses. parallel is parallel analysis: a component is kept while its eigenvalue beats the 95th percentile of what fifty random data sets of the same size produce at the same rank, so chance structure is not kept. kaiser keeps every eigenvalue above 1, which on a wide table keeps most of them. The decision is written into the model file so it can be read and argued with.	"parallel"
rotation	bool	Varimax-rotate. Off gives the raw principal axes.	True
out_scores_csv	Optional[PathLike]	Defaults: <input>_pca_scores.csv (identifiers + Component_1..Component_k) beside the input, with the model / loadings / eigenvalue files named the same way.	None
out_model_json	Optional[PathLike]	Defaults: <input>_pca_scores.csv (identifiers + Component_1..Component_k) beside the input, with the model / loadings / eigenvalue files named the same way.	None
out_loadings_csv	Optional[PathLike]	Defaults: <input>_pca_scores.csv (identifiers + Component_1..Component_k) beside the input, with the model / loadings / eigenvalue files named the same way.	None
out_eigenvalues_csv	Optional[PathLike]	Defaults: <input>_pca_scores.csv (identifiers + Component_1..Component_k) beside the input, with the model / loadings / eigenvalue files named the same way.	None
overwrite_existing	bool	If False and the scores file already exists, skip and return it.	False
encoding	str	Encoding for reading and writing CSV files.	"utf-8-sig"
rounding	int	Decimal places for scores, loadings, and eigenvalues.	4

Returns:

Type	Description
Path	out_scores_csv.

Source code in src\taters\stats\pca.py

def fit_pca_csv(
    input_csv: PathLike,
    *,
    start_col: int = 2,
    n_components: int = 0,
    retain: Literal["parallel", "kaiser"] = "parallel",
    rotation: bool = True,
    out_scores_csv: Optional[PathLike] = None,
    out_model_json: Optional[PathLike] = None,
    out_loadings_csv: Optional[PathLike] = None,
    out_eigenvalues_csv: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[[int, int], None]] = None,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
) -> Path:
    """
    Fit a varimax-rotated PCA to the numeric columns of any feature CSV.

    Parameters
    ----------
    input_csv
        Any wide feature table: identifier column(s) first, numeric feature
        columns after. Every Taters feature file (dictionary scores,
        readability, a document-term matrix) has this shape.
    start_col : int, default=2
        1-based position of the first *feature* column, exactly as in the
        house dask script this replaces. Everything before it is carried
        into the scores file unchanged as identifiers. The default fits
        tables shaped ``text_id, <features...>``.
    n_components : int, default=0
        ``0`` picks automatically by ``retain``; an explicit number is
        honored up to the count of non-constant features.
    retain : {"parallel", "kaiser"}, default="parallel"
        How ``0`` chooses. ``parallel`` is parallel analysis: a component is
        kept while its eigenvalue beats the 95th percentile of what fifty
        random data sets of the same size produce at the same rank, so
        chance structure is not kept. ``kaiser`` keeps every eigenvalue above
        1, which on a wide table keeps most of them. The decision is written
        into the model file so it can be read and argued with.
    rotation : bool, default=True
        Varimax-rotate. Off gives the raw principal axes.
    out_scores_csv, out_model_json, out_loadings_csv, out_eigenvalues_csv
        Defaults: ``<input>_pca_scores.csv`` (identifiers + ``Component_1..Component_k``)
        beside the input, with the model / loadings / eigenvalue files named
        the same way.
    overwrite_existing : bool, default=False
        If ``False`` and the scores file already exists, skip and return it.
    encoding : str, default="utf-8-sig"
        Encoding for reading and writing CSV files.
    rounding : int, default=4
        Decimal places for scores, loadings, and eigenvalues.

    Returns
    -------
    Path
        ``out_scores_csv``.
    """
    input_csv = Path(input_csv)
    if not input_csv.exists():
        raise FileNotFoundError(f"input_csv not found: {input_csv}")
    if start_col < 1:
        raise ValueError(f"start_col is 1-based; got {start_col}")

    stem = input_csv.stem
    scores_path = Path(out_scores_csv) if out_scores_csv else \
        input_csv.with_name(f"{stem}_pca_scores.csv")
    model_path = Path(out_model_json) if out_model_json else \
        input_csv.with_name(f"{stem}_pca_model.json")
    loadings_path = Path(out_loadings_csv) if out_loadings_csv else \
        input_csv.with_name(f"{stem}_pca_loadings.csv")
    eigen_path = Path(out_eigenvalues_csv) if out_eigenvalues_csv else \
        input_csv.with_name(f"{stem}_pca_eigenvalues.csv")
    scores_path.parent.mkdir(parents=True, exist_ok=True)

    from ._common import reusable

    if reusable(scores_path, input_csv, overwrite_existing=overwrite_existing,
                verbose=True, what="the PCA scores"):
        print("PCA scores output file already exists; returning existing file.")
        return scores_path

    header = _read_header(input_csv, encoding)
    if start_col > len(header):
        raise ValueError(f"start_col={start_col} but {input_csv} has only "
                         f"{len(header)} column(s).")
    id_names = header[:start_col - 1]
    features = header[start_col - 1:]

    warn_if_wide(len(features))
    n, sums, cross = stream_moments(input_csv, encoding=encoding,
                                    skip_cols=start_col - 1,
                                    on_progress=on_progress)
    if n < 3:
        raise ValueError(f"Only {n} row(s) in {input_csv}; PCA needs data.")
    announce(on_progress, "extracting components")
    kept, mu, sigma, loadings, projection, eig, pct, retention = fit_axes(
        n, sums, cross, n_components=n_components, rotation=rotation,
        retain=retain, on_progress=on_progress)
    k = projection.shape[1]
    comp_names = component_names(k)

    kept_list = kept.tolist()
    dropped = [features[j] for j in range(len(features))
               if j not in set(kept_list)]
    if dropped:
        import warnings
        warnings.warn(
            f"{len(dropped)} constant feature column(s) carried no signal "
            f"and were left out: {', '.join(dropped[:8])}"
            f"{'…' if len(dropped) > 8 else ''}")

    with atomic_write(loadings_path, newline="", encoding=encoding) as f:
        writer = csv.writer(f)
        writer.writerow(["feature", *comp_names])
        for i, j in enumerate(kept_list):
            writer.writerow([features[j],
                             *(round(float(v), rounding) for v in loadings[i])])
    with atomic_write(eigen_path, newline="", encoding=encoding) as f:
        writer = csv.writer(f)
        writer.writerow(["component", "eigenvalue", "pct_variance"])
        for name, e, p in zip(comp_names, eig, pct):
            writer.writerow([name, round(float(e), rounding),
                             round(float(p), rounding)])

    from datetime import date

    model = {
        "kind": "taters-pca-model",
        "format": PCA_MODEL_FORMAT,
        "created": date.today().isoformat(),
        "taters": taters_version(),
        "features": features,
        "model": {"n_rows": n, "rotation": bool(rotation), "retention": retention,
                  "components": comp_names, "kept": kept_list,
                  "mu": mu.tolist(), "sigma": sigma.tolist(),
                  "projection": projection.tolist(),
                  "eigenvalues": eig.tolist(),
                  "pct_variance": pct.tolist()},
    }
    with atomic_write(model_path, encoding="utf-8") as f:
        json.dump(model, f, indent=1)

    ticker = Ticker(on_progress, n)
    with input_csv.open("r", newline="", encoding=encoding) as f, \
            atomic_write(scores_path, newline="", encoding=encoding) as out:
        reader = csv.reader(f)
        next(reader)
        writer = csv.writer(out)
        writer.writerow([*id_names, *comp_names])
        for row in reader:
            ticker.tick(message="scoring components")
            writer.writerow([*row[:start_col - 1],
                             *project_row(row[start_col - 1:], kept, mu,
                                          sigma, projection, rounding)])
    return scores_path

parallel_thresholds

parallel_thresholds(
    n,
    p,
    *,
    draws=PARALLEL_DRAWS,
    percentile=PARALLEL_PERCENTILE,
    seed=0,
    chunk_rows=_CHUNK_ROWS,
    on_progress=None
)

The eigenvalues chance alone produces: for each rank, the percentile of the correlation-matrix eigenvalues of draws random standard-normal data sets with n rows and p columns.

Memory-safe by construction: no random data set is ever held whole. Each draw is generated chunk_rows rows at a time and folded into a running p-by-p cross-product, exactly as :func:stream_moments folds a real table, so the cost in memory is one more p-by-p matrix -- the same order as the fit already holds -- whatever n is. The cost in time is draws decompositions; fifty is the usual number and takes seconds at a few hundred features.

Source code in src\taters\stats\pca.py

def parallel_thresholds(n: int, p: int, *, draws: int = PARALLEL_DRAWS,
                        percentile: float = PARALLEL_PERCENTILE, seed: int = 0,
                        chunk_rows: int = _CHUNK_ROWS, on_progress=None):
    """
    The eigenvalues chance alone produces: for each rank, the ``percentile``
    of the correlation-matrix eigenvalues of ``draws`` random standard-normal
    data sets with ``n`` rows and ``p`` columns.

    Memory-safe by construction: no random data set is ever held whole. Each
    draw is generated ``chunk_rows`` rows at a time and folded into a running
    ``p``-by-``p`` cross-product, exactly as :func:`stream_moments` folds a
    real table, so the cost in memory is one more ``p``-by-``p`` matrix --
    the same order as the fit already holds -- whatever ``n`` is. The cost in
    time is ``draws`` decompositions; fifty is the usual number and takes
    seconds at a few hundred features.
    """
    import numpy as np

    rng = np.random.default_rng(seed)
    n, p = int(n), int(p)
    eigs = np.empty((draws, p), dtype=np.float64)
    for d in range(draws):
        if on_progress is not None:
            on_progress(d, draws, "parallel analysis: eigenvalues of random data")
        sums = np.zeros(p)
        cross = np.zeros((p, p))
        done = 0
        while done < n:
            rows = min(chunk_rows, n - done)
            block = rng.standard_normal((rows, p))
            sums += block.sum(axis=0)
            cross += block.T @ block
            done += rows
        mu = sums / n
        cov = cross / n - np.outer(mu, mu)
        sd = np.sqrt(np.maximum(np.diag(cov), 1e-12))
        corr = cov / np.outer(sd, sd)
        eigs[d] = np.linalg.eigvalsh((corr + corr.T) / 2)[::-1]
    return np.percentile(eigs, percentile, axis=0)

project_row

project_row(cells, kept, mu, sigma, projection, rounding)

One row's component scores from raw feature cells -- shared by every fit and every apply, so the two can never drift.

Source code in src\taters\stats\pca.py

def project_row(cells, kept, mu, sigma, projection, rounding: int) -> list:
    """One row's component scores from raw feature cells -- shared by every
    fit and every apply, so the two can never drift."""
    import numpy as np

    x = np.asarray(cells, dtype=np.float64)
    scores = ((x[kept] - mu) / sigma) @ projection
    return [round(float(s), rounding) for s in scores]

stream_moments

stream_moments(
    csv_path,
    *,
    encoding,
    skip_cols,
    on_progress=None,
    message="measuring the table"
)

One pass over the rows: n, per-column sums, and X^T X for the columns after the first skip_cols. Rows are accumulated in chunks so the cross-product is a handful of BLAS calls per chunk rather than an outer product per row -- same numbers, a large constant factor faster.

Source code in src\taters\stats\pca.py

def stream_moments(csv_path: Path, *, encoding: str, skip_cols: int,
                   on_progress=None,
                   message: str = "measuring the table") -> Tuple[int, "object", "object"]:
    """
    One pass over the rows: n, per-column sums, and X^T X for the columns
    after the first ``skip_cols``. Rows are accumulated in chunks so the
    cross-product is a handful of BLAS calls per chunk rather than an outer
    product per row -- same numbers, a large constant factor faster.
    """
    import numpy as np

    sums = cross = None
    n = 0
    chunk: List[List[str]] = []

    def _flush():
        nonlocal sums, cross, n
        if not chunk:
            return
        block = np.asarray(chunk, dtype=np.float64)
        if sums is None:
            sums = np.zeros(block.shape[1], dtype=np.float64)
            cross = np.zeros((block.shape[1], block.shape[1]), dtype=np.float64)
        sums += block.sum(axis=0)
        cross += block.T @ block
        n += block.shape[0]
        chunk.clear()

    ticker = Ticker(on_progress, count_rows(csv_path, on_progress=on_progress))
    with Path(csv_path).open("r", newline="", encoding=encoding) as f:
        reader = csv.reader(f)
        next(reader)                                    # skip the header; caller's job
        for row in reader:
            ticker.tick(message=message)
            chunk.append(row[skip_cols:])
            if len(chunk) >= _CHUNK_ROWS:
                _flush()
        _flush()
    return n, sums, cross

varimax

varimax(
    loadings, gamma=1.0, q=2000, tol=1e-13, normalize=True
)

Orthogonal varimax rotation of an (n_features, n_components) loadings matrix; returns (rotated_loadings, rotation_matrix).

normalize=True is Kaiser normalization -- each feature's loading row is scaled to unit communality while the rotation is chosen, which is what stats::varimax, psych::principal and SPSS all do by default. Without it, high-communality features dominate the criterion and the solution genuinely differs (a real MEM run diverged from its R twin until this matched; the worst theme correlated at 0.39 unnormalized and 0.996 normalized). Normalization only steers the choice of rotation: the returned loadings are still loadings @ R exactly, so a scorer using the same R stays coherent with them.

The iteration budget is a ceiling, not a schedule -- the loop stops on convergence, and the predecessor's cap of 20 sweeps was itself a source of divergence at MEM sizes (1,000 terms x 25 themes converges in ~100).

Source code in src\taters\stats\pca.py

def varimax(loadings, gamma: float = 1.0, q: int = 2000, tol: float = 1e-13,
            normalize: bool = True):
    """
    Orthogonal varimax rotation of an (n_features, n_components) loadings
    matrix; returns (rotated_loadings, rotation_matrix).

    ``normalize=True`` is **Kaiser normalization** -- each feature's loading
    row is scaled to unit communality while the rotation is *chosen*, which
    is what ``stats::varimax``, ``psych::principal`` and SPSS all do by
    default. Without it, high-communality features dominate the criterion
    and the solution genuinely differs (a real MEM run diverged from its R
    twin until this matched; the worst theme correlated at 0.39 unnormalized
    and 0.996 normalized). Normalization only steers the choice of rotation:
    the returned loadings are still ``loadings @ R`` exactly, so a scorer
    using the same ``R`` stays coherent with them.

    The iteration budget is a ceiling, not a schedule -- the loop stops on
    convergence, and the predecessor's cap of 20 sweeps was itself a source
    of divergence at MEM sizes (1,000 terms x 25 themes converges in ~100).
    """
    import numpy as np

    p, k = loadings.shape
    if normalize:
        h = np.sqrt((loadings ** 2).sum(axis=1))
        h[h == 0] = 1.0
        working = loadings / h[:, None]
    else:
        working = loadings
    R = np.eye(k)
    d = 0.0
    converged = False
    for _ in range(q):
        Lambda = working @ R
        u, s, vh = np.linalg.svd(
            working.T @ (Lambda ** 3 - (gamma / p) * Lambda
                         @ np.diag(np.sum(Lambda ** 2, axis=0)))
        )
        R_new = u @ vh
        d_old = d
        d = float(np.sum(s))
        R = R_new
        if d_old != 0 and d / d_old < 1 + tol:
            converged = True
            break
    if not converged:
        # running out of sweeps, in practice, means near-tied components: the
        # criterion is almost flat and no rotation is uniquely "the" varimax
        # solution. the rotation we hand back is still exact and
        # deterministic, but quietly passing it off as converged is exactly
        # the kind of plausible-looking-numbers failure we're trying to
        # avoid, so we say so.
        import warnings
        warnings.warn(
            f"varimax did not converge within {q} sweeps; the rotation "
            "returned is the best found. Near-tied components usually cause "
            "this -- the themes are valid but their exact rotation is not "
            "unique, and a different tool may land elsewhere. Fewer "
            "components often stabilizes it.")
    return loadings @ R, R

warn_if_wide

warn_if_wide(n_features, *, warn_gb=1.0)

Memory truth-telling before the work starts.

The fit is memory-safe against corpus LENGTH by construction -- rows stream through in chunks and are gone. Its one memory cost grows with table WIDTH: the p-by-p cross-product/correlation matrix plus the eigendecomposition's workspace, roughly three p*p float64 buffers. At the default 250-term vocabulary that is under 2 MB; at 5,000 features about 0.6 GB; at 20,000 it would be ~10 GB, which deserves a warning before hours of streaming, not an OOM after.

Source code in src\taters\stats\pca.py

def warn_if_wide(n_features: int, *, warn_gb: float = 1.0) -> None:
    """
    Memory truth-telling before the work starts.

    The fit is memory-safe against corpus LENGTH by construction -- rows
    stream through in chunks and are gone. Its one memory cost grows with
    table WIDTH: the p-by-p cross-product/correlation matrix plus the
    eigendecomposition's workspace, roughly three p*p float64 buffers. At the
    default 250-term vocabulary that is under 2 MB; at 5,000 features about
    0.6 GB; at 20,000 it would be ~10 GB, which deserves a warning before
    hours of streaming, not an OOM after.
    """
    approx_gb = 3 * (n_features ** 2) * 8 / 1024 ** 3
    if approx_gb >= warn_gb:
        import warnings
        warnings.warn(
            f"{n_features} feature columns: the exact fit holds about "
            f"{approx_gb:.1f} GB of p-by-p working memory. Consider fewer "
            "features (e.g. a smaller vocab_top_n).")

taters.stats.describe

Descriptive statistics for every feature table a run writes.

The first thing anyone does with a new measure is look at its distribution: how many texts have it, what its mean and spread are, whether it is skewed into a corner, whether a suspicious maximum is one text or a hundred. Until now that meant opening each CSV in a spreadsheet and doing it by hand, per column, per table. This step does it once for every numeric column of every feature table -- count, missing, mean, standard deviation, quartiles, range, skewness, kurtosis, zeros, distinct values -- and writes one table of descriptives per feature table under stats_descriptives/, whether or not any statistics were asked for.

The numbers are the textbook sample statistics (standard deviation with n-1; the bias-adjusted skewness and excess kurtosis that SPSS, R's psych::describe with type=2 and Excel report), so they can be pasted into a methods table beside anyone else's. Every cell is computed from the values that are present; a blank cell is missing, and is counted as such.

One file per table, written only when the table is newer than it, so a second run of an unchanged pipeline writes nothing -- the resume contract.

describe_features

describe_features(
    feature_csvs,
    out_dir="stats_descriptives",
    *,
    overwrite_existing=False,
    encoding="utf-8-sig",
    rounding=4,
    verbose=True,
    on_progress=None
)

Write a table of descriptive statistics for every feature table.

Parameters:

Name	Type	Description	Default
feature_csvs	sequence of str or Path	The feature tables to describe -- one row per text, identifier columns first, measures after. Every Taters feature file has this shape; the pipeline hands over every table the run wrote.	required
out_dir	str or Path	Where the descriptives go: <out_dir>/<table stem>.csv for each table, and a README.md listing them.	"stats_descriptives"
overwrite_existing	bool	Recompute a table's descriptives even when they are newer than it.	False
encoding	str	The tables' encoding, and the descriptives'.	"utf-8-sig"
rounding	int	Decimal places written.	4
verbose	bool	The usual.	True
on_progress	bool	The usual.	True

Returns:

Type	Description
Path	out_dir.

Notes

Every numeric column is described except the identifiers (text_id, source, speaker). Columns holding words are left out, since a mean of them means nothing. Count columns the steps write beside their measures (token_count, WC) are described like any other -- they are the first thing to look at when a measure looks strange.

Source code in src\taters\stats\describe.py

def describe_features(
    feature_csvs: Sequence[PathLike],
    out_dir: PathLike = "stats_descriptives",
    *,
    overwrite_existing: bool = False,
    encoding: str = "utf-8-sig",
    rounding: int = 4,
    verbose: bool = True,
    on_progress: Optional[Callable[..., None]] = None,
) -> Path:
    """
    Write a table of descriptive statistics for every feature table.

    Parameters
    ----------
    feature_csvs : sequence of str or Path
        The feature tables to describe -- one row per text, identifier
        columns first, measures after. Every Taters feature file has this
        shape; the pipeline hands over every table the run wrote.
    out_dir : str or Path, default "stats_descriptives"
        Where the descriptives go: ``<out_dir>/<table stem>.csv`` for each
        table, and a ``README.md`` listing them.
    overwrite_existing : bool, default False
        Recompute a table's descriptives even when they are newer than it.
    encoding : str, default "utf-8-sig"
        The tables' encoding, and the descriptives'.
    rounding : int, default 4
        Decimal places written.
    verbose, on_progress
        The usual.

    Returns
    -------
    Path
        ``out_dir``.

    Notes
    -----
    Every numeric column is described except the identifiers (``text_id``,
    ``source``, ``speaker``). Columns holding words are left out, since a
    mean of them means nothing. Count columns the steps write beside their
    measures (``token_count``, ``WC``) are described like any other -- they
    are the first thing to look at when a measure looks strange.
    """
    out_dir = Path(out_dir)
    paths = [Path(p) for p in feature_csvs]
    if not paths:
        if verbose:
            print("[describe] no feature tables to describe.")
        return out_dir
    written = reused = 0
    index: List[str] = []
    for i, path in enumerate(paths):
        if on_progress is not None:
            on_progress(i, len(paths), f"describing {path.name}")
        if not path.is_file():
            raise FileNotFoundError(f"feature table not found: {path}")
        out = out_dir / f"{path.stem}.csv"
        if reusable(out, path, overwrite_existing=overwrite_existing,
                    verbose=verbose, what=f"the descriptives for {path.name}"):
            reused += 1
        else:
            rows = _describe_table(path, encoding)
            out_dir.mkdir(parents=True, exist_ok=True)
            with atomic_write(out, mode="w", newline="", encoding=encoding) as fh:
                writer = csv.writer(fh)
                writer.writerow(COLUMNS)
                for row in rows:
                    writer.writerow([row[0]] + [_fmt(v, rounding) for v in row[1:]])
            written += 1
        index.append(out.name)
    announce(on_progress, "writing the descriptives index")
    lines = ["# Descriptive statistics", "",
             "One table per feature table: for every numeric measure, how many "
             "texts have it (`n`) and how many do not (`missing`), the mean and "
             "standard deviation (n-1), the minimum, quartiles and maximum, the "
             "skewness and excess kurtosis (bias-adjusted, as SPSS and R's "
             "psych report them; zero for a normal distribution), how many "
             "values are exactly zero, and how many distinct values there are. "
             "Identifier columns are left out.", ""]
    lines += [f"- `{name}`" for name in index]
    _write_if_changed(out_dir / "README.md", "\n".join(lines))
    if verbose:
        print(f"[describe] {written} table(s) described, {reused} already "
              f"current -> {out_dir}")
    return out_dir

describe_values

describe_values(values)

The descriptives of one column, from its values with blanks as NaN.

sd needs two values, skew three and kurtosis four; below that they are None rather than a division by zero dressed as a number. Skewness is the bias-adjusted sample statistic (G1) and kurtosis the bias-adjusted excess kurtosis (G2), zero for a normal distribution -- the pair a methods table expects.

Source code in src\taters\stats\describe.py

def describe_values(values) -> Dict[str, Optional[float]]:
    """
    The descriptives of one column, from its values with blanks as NaN.

    ``sd`` needs two values, ``skew`` three and ``kurtosis`` four; below
    that they are ``None`` rather than a division by zero dressed as a
    number. Skewness is the bias-adjusted sample statistic (G1) and kurtosis
    the bias-adjusted *excess* kurtosis (G2), zero for a normal distribution
    -- the pair a methods table expects.
    """
    import numpy as np

    x = np.asarray(values, dtype=float)
    present = x[~np.isnan(x)]
    n = int(present.size)
    out: Dict[str, Optional[float]] = {
        "n": n, "missing": int(x.size - n), "mean": None, "sd": None,
        "min": None, "q1": None, "median": None, "q3": None, "max": None,
        "skew": None, "kurtosis": None, "zeros": int(np.sum(present == 0)),
        "unique": int(np.unique(present).size),
    }
    if n == 0:
        return out
    mean = float(present.mean())
    out["mean"] = mean
    q1, median, q3 = (float(v) for v in np.percentile(present, [25, 50, 75]))
    out.update(min=float(present.min()), q1=q1, median=median, q3=q3,
               max=float(present.max()))
    if n >= 2:
        sd = float(present.std(ddof=1))
        out["sd"] = sd
        if sd > 0:
            z = (present - mean) / sd
            if n >= 3:
                out["skew"] = float(np.sum(z ** 3) * n / ((n - 1) * (n - 2)))
            if n >= 4:
                m4 = float(np.sum(z ** 4))
                out["kurtosis"] = float(
                    (n * (n + 1) * m4 / ((n - 1) * (n - 2) * (n - 3)))
                    - 3 * (n - 1) ** 2 / ((n - 2) * (n - 3)))
        else:
            out["skew"] = 0.0 if n >= 3 else None
            out["kurtosis"] = None
    return out

taters.score_model

Score this dataset with a model somebody already fitted.

A saved model is the one artifact in Taters meant to leave the run that made it: fit a ridge on a corpus where you have the outcome, then apply it to a corpus where you do not. That is one job as far as the researcher is concerned, and it was three -- a MEM topic model, a ridge and a classifier each had their own menu entry, their own vocabulary and their own place in the pipeline, and the ridge's place was at the very end, among the statistics, where it is least useful: what comes back is a feature, a column per text, and the next thing anyone wants to do with it is analyze it like any other.

So there is one entry, it lives with feature extraction, and it works out what to do from the model file. What a model needs is a fact about the model, not a question for the user:

• a MEM topic model needs text, and re-derives its own vocabulary from it;
• a ridge or a classifier needs the feature columns it was fitted on, by name, which means the steps that produce them have to have run.

The second case is the interesting one, and the honest thing to do when those columns are absent is to name them and stop -- not to impute, and not to score a model on a feature it has never seen and call the result a prediction.

What this does not yet do

A model fitted on features that were themselves transformed -- sentence embeddings reduced to fifty dimensions by PCA, say -- needs that whole recipe replayed before its weights mean anything, and a model file records its predictor names but not how they were made. So the columns have to already exist under the names the model knows. Getting from "here is a folder of text" to "here are the exact fifty numbers this model wants" is a pipeline the model would have to carry with it, and that is not built yet; until it is, this step tells you precisely which columns it wanted.

score_with_model

score_with_model(
    *,
    model_json,
    feature_csvs=(),
    metadata_csv=None,
    key_cols=("text_id",),
    keep_inputs=False,
    unverified_ok=(),
    allow_unrecorded=False,
    csv_path=None,
    txt_dir=None,
    analysis_csv=None,
    gathered_csv=None,
    out_csv=None,
    overwrite_existing=False,
    on_progress=None,
    verbose=True,
    workers=None,
    device="auto",
    encoding="utf-8-sig",
    rounding=6,
    text_cols=("text",),
    id_cols=None,
    mode="concat",
    group_by=None,
    delimiter=",",
    joiner=" ",
    num_buckets=512,
    max_open_bucket_files=64,
    tmp_root=None,
    recursive=True,
    pattern=None,
    id_from="stem",
    include_source_path=True
)

Score this dataset with a saved model, whatever kind of model it is.

Parameters:

Name	Type	Description	Default
model_json	Union[PathLike, Sequence[PathLike]]	One or more models saved by earlier runs -- topic models, ridges, classifiers, word vectors, fine-tuned predictors -- as a file, a list of files, or a folder holding them (every model in it). Each model says what kind it is and what it needs; nothing here has to be told. Several models are scored one after another, each against the feature tables it was fitted on, into one file per model and one merged table; two models whose names read the same are refused (rename one under Settings → Manage saved models).	required
feature_csvs	Sequence[PathLike]	The feature tables this run produced, joined on key_cols to find the columns a ridge or classifier was fitted on. Ignored by a model that works from text.	()
metadata_csv	Optional[PathLike]	A table keyed like the feature tables that carries the spreadsheet's own columns -- needed only by a model fitted with controls, whose age or gender no feature step produces. In a pipeline the metadata gather writes it; the composer adds that step whenever a chosen model needs controls.	None
key_cols	Sequence[str]	The columns identifying a row, used to join the feature tables.	('text_id',)
keep_inputs	bool	Keep the joined feature table the scoring read. Off by default: the row accounting that explains an unscored row is reported either way, and the table itself is a copy of columns that already exist.	False
unverified_ok	Sequence[str]	Feature tables to score against even though their measuring settings cannot be checked -- named one at a time, never a blanket flag, because a single boolean would let you waive one unverifiable table and silently waive the acoustics table sitting next to it. This waives not knowing. It cannot waive a known mismatch: when both the model and the table state their settings and the settings differ, scoring is refused outright.	()
allow_unrecorded	bool	Score a model that predates settings-recording at all. Such a model cannot be checked in either direction, so this is refused by default: it is exactly the case that reported a prediction seven years wrong without complaint.	False
csv_path	Optional[PathLike]	Where the text comes from, for a model that reads text: a spreadsheet, a folder of .txt files, or an already-gathered analysis-ready table. Exactly one applies; a model that reads features ignores all four.	None
txt_dir	Optional[PathLike]	Where the text comes from, for a model that reads text: a spreadsheet, a folder of .txt files, or an already-gathered analysis-ready table. Exactly one applies; a model that reads features ignores all four.	None
analysis_csv	Optional[PathLike]	Where the text comes from, for a model that reads text: a spreadsheet, a folder of .txt files, or an already-gathered analysis-ready table. Exactly one applies; a model that reads features ignores all four.	None
gathered_csv	Optional[PathLike]	Where the text comes from, for a model that reads text: a spreadsheet, a folder of .txt files, or an already-gathered analysis-ready table. Exactly one applies; a model that reads features ignores all four.	None
out_csv	Optional[PathLike]	Where the scores go. With one model, its scores; default features/model_scores/<model name>.csv. With several, the merged table -- an outer join of every model's scores on key_cols, every score column prefixed <model name>__ so two models predicting the same outcome cannot collide -- default features/model_scores.csv, with each model's own file beside it under <out_csv stem>/<model name>.csv and its own unscored-row accounting; the merged _unscored.csv names the model per row.	None
workers	int	Parallelism for a text model, which has real work to do per row. Ignored by the others, whose arithmetic is one matrix multiply.	None
device	('auto', 'cuda', 'cpu')	Where a text model's tagger runs, if it was built with the stanza engine -- a runtime choice, deliberately not stored in the model.	"auto"
overwrite_existing	bool	When False (default) and the scores exist, they are returned untouched.	False
rounding	int	Decimal places in the scores.	6
text_cols	Sequence[str]	How to read a spreadsheet of text: which columns hold it, which identify a row, and whether several text columns are joined ("concat") or scored one at a time ("separate"). Only used with csv_path.	('text',)
id_cols	Sequence[str]	How to read a spreadsheet of text: which columns hold it, which identify a row, and whether several text columns are joined ("concat") or scored one at a time ("separate"). Only used with csv_path.	('text',)
mode	Sequence[str]	How to read a spreadsheet of text: which columns hold it, which identify a row, and whether several text columns are joined ("concat") or scored one at a time ("separate"). Only used with csv_path.	('text',)
group_by	Sequence[str]	How to read a spreadsheet of text: which columns hold it, which identify a row, and whether several text columns are joined ("concat") or scored one at a time ("separate"). Only used with csv_path.	('text',)
delimiter	Sequence[str]	How to read a spreadsheet of text: which columns hold it, which identify a row, and whether several text columns are joined ("concat") or scored one at a time ("separate"). Only used with csv_path.	('text',)
joiner	Sequence[str]	How to read a spreadsheet of text: which columns hold it, which identify a row, and whether several text columns are joined ("concat") or scored one at a time ("separate"). Only used with csv_path.	('text',)
num_buckets	int	Spill settings for grouping a large spreadsheet, passed straight through to the gatherer.	512
max_open_bucket_files	int	Spill settings for grouping a large spreadsheet, passed straight through to the gatherer.	512
tmp_root	int	Spill settings for grouping a large spreadsheet, passed straight through to the gatherer.	512
recursive	bool	How to read a folder of .txt files. Only used with txt_dir.	True
pattern	bool	How to read a folder of .txt files. Only used with txt_dir.	True
id_from	bool	How to read a folder of .txt files. Only used with txt_dir.	True
include_source_path	bool	How to read a folder of .txt files. Only used with txt_dir.	True

Returns:

Type	Description
Path	The scores: one row per text, with the model's own output columns. A row the model could not score is blank rather than absent, so the table still joins back to everything else.

Raises:

Type	Description
ValueError	When the model needs feature columns this run does not have. The message names them, because the fix is to add the steps that produce them and nothing else in the message helps with that.

Source code in src\taters\score_model.py

def score_with_model(
    *,
    model_json: Union[PathLike, Sequence[PathLike]],

    # ----- what a feature-reading model scores -----
    feature_csvs: Sequence[PathLike] = (),
    metadata_csv: Optional[PathLike] = None,
    key_cols: Sequence[str] = ("text_id",),
    keep_inputs: bool = False,
    unverified_ok: Sequence[str] = (),
    allow_unrecorded: bool = False,

    # ----- what a text-reading model scores. pick one of these three; the
    # ----- arguments belonging to the other two get ignored. this is the
    # ----- same input contract every text step in Taters has, so that the
    # ----- composer can wire this one up exactly like the rest.
    csv_path: Optional[PathLike] = None,
    txt_dir: Optional[PathLike] = None,
    analysis_csv: Optional[PathLike] = None,
    gathered_csv: Optional[PathLike] = None,

    # ----- output -----
    out_csv: Optional[PathLike] = None,
    overwrite_existing: bool = False,
    on_progress: Optional[Callable[..., None]] = None,
    verbose: bool = True,
    workers: Optional[int] = None,
    device: str = "auto",
    encoding: str = "utf-8-sig",
    rounding: int = 6,

    # ====== CSV GATHER OPTIONS (used when csv_path is provided) ======
    text_cols: Sequence[str] = ("text",),
    id_cols: Optional[Sequence[str]] = None,
    mode: str = "concat",
    group_by: Optional[Sequence[str]] = None,
    delimiter: str = ",",
    joiner: str = " ",
    num_buckets: int = 512,
    max_open_bucket_files: int = 64,
    tmp_root: Optional[PathLike] = None,

    # ====== TXT FOLDER GATHER OPTIONS (used when txt_dir is provided) ======
    recursive: bool = True,
    pattern: Optional[str] = None,
    id_from: str = "stem",
    include_source_path: bool = True,
) -> Path:
    """
    Score this dataset with a saved model, whatever kind of model it is.

    Parameters
    ----------
    model_json
        One or more models saved by earlier runs -- topic models, ridges,
        classifiers, word vectors, fine-tuned predictors -- as a file, a
        list of files, or a folder holding them (every model in it). Each
        model says what kind it is and what it needs; nothing here has to
        be told. Several models are scored one after another, each against
        the feature tables *it* was fitted on, into one file per model and
        one merged table; two models whose names read the same are refused
        (rename one under Settings → Manage saved models).
    feature_csvs
        The feature tables this run produced, joined on ``key_cols`` to find
        the columns a ridge or classifier was fitted on. Ignored by a model
        that works from text.
    metadata_csv
        A table keyed like the feature tables that carries the spreadsheet's
        own columns -- needed only by a model fitted with controls, whose
        ``age`` or ``gender`` no feature step produces. In a pipeline the
        metadata gather writes it; the composer adds that step whenever a
        chosen model needs controls.
    key_cols
        The columns identifying a row, used to join the feature tables.
    keep_inputs
        Keep the joined feature table the scoring read. Off by default: the
        row accounting that explains an unscored row is reported either way,
        and the table itself is a copy of columns that already exist.
    unverified_ok
        Feature tables to score against even though their measuring settings
        cannot be checked -- named one at a time, never a blanket flag,
        because a single boolean would let you waive one unverifiable table
        and silently waive the acoustics table sitting next to it.

        This waives *not knowing*. It cannot waive a known mismatch: when
        both the model and the table state their settings and the settings
        differ, scoring is refused outright.
    allow_unrecorded
        Score a model that predates settings-recording at all. Such a model
        cannot be checked in either direction, so this is refused by default:
        it is exactly the case that reported a prediction seven years wrong
        without complaint.
    csv_path, txt_dir, analysis_csv, gathered_csv
        Where the text comes from, for a model that reads text: a
        spreadsheet, a folder of ``.txt`` files, or an already-gathered
        analysis-ready table. Exactly one applies; a model that reads
        features ignores all four.
    out_csv
        Where the scores go. With one model, its scores; default
        ``features/model_scores/<model name>.csv``. With several, the
        **merged** table -- an outer join of every model's scores on
        ``key_cols``, every score column prefixed ``<model name>__`` so
        two models predicting the same outcome cannot collide -- default
        ``features/model_scores.csv``, with each model's own file beside it
        under ``<out_csv stem>/<model name>.csv`` and its own unscored-row
        accounting; the merged ``_unscored.csv`` names the model per row.
    workers : int, optional
        Parallelism for a text model, which has real work to do per row.
        Ignored by the others, whose arithmetic is one matrix multiply.
    device : {"auto", "cuda", "cpu"}
        Where a text model's tagger runs, if it was built with the stanza
        engine -- a runtime choice, deliberately not stored in the model.
    overwrite_existing
        When False (default) and the scores exist, they are returned
        untouched.
    rounding
        Decimal places in the scores.
    text_cols, id_cols, mode, group_by, delimiter, joiner
        How to read a spreadsheet of text: which columns hold it, which
        identify a row, and whether several text columns are joined
        (``"concat"``) or scored one at a time (``"separate"``). Only used
        with ``csv_path``.
    num_buckets, max_open_bucket_files, tmp_root
        Spill settings for grouping a large spreadsheet, passed straight
        through to the gatherer.
    recursive, pattern, id_from, include_source_path
        How to read a folder of ``.txt`` files. Only used with ``txt_dir``.

    Returns
    -------
    pathlib.Path
        The scores: one row per text, with the model's own output columns.
        A row the model could not score is blank rather than absent, so the
        table still joins back to everything else.

    Raises
    ------
    ValueError
        When the model needs feature columns this run does not have. The
        message names them, because the fix is to add the steps that produce
        them and nothing else in the message helps with that.
    """
    from .helpers.progress import announce
    from .stats._common import reusable

    infos = _models(model_json)
    several = len(infos) > 1
    if not several:
        info = infos[0]
        out_path = Path(out_csv) if out_csv else (
            Path("features") / "model_scores" / f"{slug(info.name)}.csv")
        tables = _tables_for(info, feature_csvs)
        # gate first, *then* check whether the scores already exist. if we
        # did it the other way around, a stale scores file sitting next to a
        # feature table that somebody re-measured with different settings
        # would sail right through, and that's exactly how this bites people
        # in real life (re-running a results folder after tweaking something)
        _check_provenance(info, tables, unverified_ok=unverified_ok,
                          allow_unrecorded=allow_unrecorded, verbose=verbose)
        return _score_one(
            info, tables, out_path, metadata_csv=metadata_csv, key_cols=key_cols,
            keep_inputs=keep_inputs, overwrite_existing=overwrite_existing,
            on_progress=on_progress, verbose=verbose, encoding=encoding,
            rounding=rounding, workers=workers, device=device,
            text_input=dict(
                csv_path=csv_path, txt_dir=txt_dir, analysis_csv=analysis_csv,
                gathered_csv=gathered_csv, text_cols=text_cols, id_cols=id_cols,
                mode=mode, group_by=group_by, delimiter=delimiter, joiner=joiner,
                num_buckets=num_buckets, max_open_bucket_files=max_open_bucket_files,
                tmp_root=tmp_root, recursive=recursive, pattern=pattern,
                id_from=id_from, include_source_path=include_source_path))

    # several models. the order here matters: we run every model's gate
    # before we score anything, because nobody wants a five-model run to
    # die on model #4 after chewing through the first three. so we pick each
    # model's tables, check them all, and complain about everything at once
    merged_path = Path(out_csv) if out_csv else Path("features") / "model_scores.csv"
    per_dir = merged_path.parent / merged_path.stem
    chosen_tables = []
    problems = []
    for info in infos:
        try:
            tables = _tables_for(info, feature_csvs)
            _check_provenance(info, tables, unverified_ok=unverified_ok,
                              allow_unrecorded=allow_unrecorded, verbose=verbose)
        except ValueError as e:
            problems.append(str(e))
            tables = []
        chosen_tables.append(tables)
    if problems:
        if len(problems) == 1:
            raise ValueError(problems[0])
        raise ValueError(
            f"{len(problems)} of the {len(infos)} models cannot be scored on this "
            f"run's features:\n- " + "\n- ".join(problems))

    text_input = dict(
        csv_path=csv_path, txt_dir=txt_dir, analysis_csv=analysis_csv,
        gathered_csv=gathered_csv, text_cols=text_cols, id_cols=id_cols,
        mode=mode, group_by=group_by, delimiter=delimiter, joiner=joiner,
        num_buckets=num_buckets, max_open_bucket_files=max_open_bucket_files,
        tmp_root=tmp_root, recursive=recursive, pattern=pattern, id_from=id_from,
        include_source_path=include_source_path)
    per_paths = []
    for info, tables in zip(infos, chosen_tables):
        per_path = per_dir / f"{slug(info.name)}.csv"
        per_paths.append(_score_one(
            info, tables, per_path, metadata_csv=metadata_csv, key_cols=key_cols,
            keep_inputs=keep_inputs, overwrite_existing=overwrite_existing,
            on_progress=on_progress, verbose=verbose, encoding=encoding,
            rounding=rounding, workers=workers, device=device,
            text_input=text_input))

    merged_path.parent.mkdir(parents=True, exist_ok=True)
    if reusable(merged_path, *per_paths, *(i.path for i in infos),
                overwrite_existing=overwrite_existing, verbose=verbose,
                what="the merged model scores"):
        return merged_path
    announce(on_progress, f"merging the scores of {len(infos)} models")
    _merge_scores(infos, per_paths, merged_path, key_cols=key_cols, encoding=encoding)
    _merge_unscored(infos, per_paths, merged_path, encoding=encoding)
    if verbose:
        print(f"[score] {len(infos)} models scored; each in {per_dir}, all in "
              f"{merged_path} with columns prefixed by the model's name")
    return merged_path

taters.figures.wordclouds

Which words go into which cloud, for every kind of result Taters writes.

Three entry points, one per stage of a run:

• :func:stats_wordclouds reads the tables under a stats_results folder -- ridge and classifier coefficients, correlations, pairwise group differences, PCA loadings, and the analysis table itself for per-group frequencies -- and draws one cloud per direction per outcome, plus a report section that shows them.
• :func:theme_wordclouds draws one cloud per theme of a topic model from its loadings table.
• :func:frequency_wordclouds draws the corpus's most frequent words from a frequency list.

Every function reads finished tables and computes no statistic of its own. The choosing -- which rows, which model, which threshold -- is in pure functions over lists of rows, tested without a picture ever being drawn; the drawing is :mod:taters.figures.render. A cloud that would be empty is not drawn and the report says why, because a missing picture reads as a bug and "no feature passed p < .05" reads as a finding.

The folders follow the analyses: <analysis>/<feature set>/ holds that set's clouds, with components/ and themes/ inside it for what its predictors are made of -- a ridge over topic-model themes shows "Theme_5" as a word, and the picture that says what Theme 5 is belongs beside it, not three folders away.

Pictures are redrawn only when the table they came from is newer than they are, the same rule every statistics step follows, and the report section is rewritten only when its text changes -- so a second run of an unchanged pipeline writes nothing, which is the resume contract the whole pipeline keeps.

CloudSpec dataclass

CloudSpec(
    analysis,
    name,
    title,
    legend,
    words=list(),
    note="",
    folder="",
    set_name="",
    stat="",
    kind="cloud",
    inputs=list(),
)

One picture to draw, or one sentence about why there is none.

folder is where it goes under the figures root and name the file name without extension; analysis the heading it lists under in the report and set_name the feature set it belongs to; words the (label, weight) pairs; note replaces the picture when words is empty; stat names the statistic the weights are, for the sentence that says where a theme appeared ("TIPI_Open: β = +0.32"); kind is cloud, component or theme; inputs the tables it was drawn from, which decide whether it needs redrawing.

relative_png property

relative_png

The picture's path under the figures root, POSIX style.

clouds_for_classifier

clouds_for_classifier(rows, allowed)

Two clouds per class per outcome per feature set from classifier_coefficients.csv: what pushes a text toward the class, and what pushes it away.

Source code in src\taters\figures\wordclouds.py

def clouds_for_classifier(rows: Sequence[Row], allowed: Callable[[str, str], bool]
                          ) -> List[CloudSpec]:
    """Two clouds per class per outcome per feature set from
    ``classifier_coefficients.csv``: what pushes a text toward the class,
    and what pushes it away."""
    if not rows:
        return []
    heading, folder = ANALYSES["classifier"]
    split = _split_col(rows)
    out: List[CloudSpec] = []
    keyed = _grouped(rows, split or "", "feature_set", "outcome", "class")
    for (sv, fs, outcome, klass), group in keyed.items():
        if "+" in fs:
            continue        # a combination of tables; its members get the clouds
        model_rows = _full_model(group)
        words = [(_label(r["predictor"], fs), _num(r.get("coef")))
                 for r in model_rows if allowed(fs, r["predictor"])]
        pos, neg = _signed([(t, w) for t, w in words if w])
        where = _where(split, sv)
        out += _pair(
            heading, folder, fs,
            f"{slug(outcome)}__{slug(klass)}{_tag(split, sv)}", pos, neg,
            titles=(f"Classifier: {outcome} = {klass} — features pushing toward "
                    f"this class ({fs}){where}",
                    f"Classifier: {outcome} = {klass} — features pushing away "
                    f"from this class ({fs}){where}"),
            legends=(f"Bigger and darker = larger standardized coefficient. "
                     f"Blue: pushes toward {klass}.",
                     f"Bigger and darker = larger standardized coefficient. "
                     f"Red: pushes away from {klass}."),
            notes=(f"no feature pushed toward {outcome} = {klass} ({fs}).",
                   f"no feature pushed away from {outcome} = {klass} ({fs})."),
            stat=f"{outcome} = {klass}: coefficient")
    return out

clouds_for_components

clouds_for_components(rows, *, analysis, component_words)

One mixed-sign cloud per component from a *_pca_loadings*.csv: the component_words features that load most strongly, blue for a positive loading and red for a negative one. Filed under the analysis and set they belong to, in a components/ folder of their own.

Source code in src\taters\figures\wordclouds.py

def clouds_for_components(rows: Sequence[Row], *, analysis: str,
                          component_words: int) -> List[CloudSpec]:
    """One mixed-sign cloud per component from a ``*_pca_loadings*.csv``:
    the ``component_words`` features that load most strongly, blue for a
    positive loading and red for a negative one. Filed under the analysis
    and set they belong to, in a ``components/`` folder of their own."""
    if not rows:
        return []
    heading, folder = ANALYSES.get(
        analysis, (analysis.replace("_", " ").capitalize(), slug(analysis)))
    cols = list(rows[0].keys())
    components = cols[cols.index("feature") + 1:]
    out: List[CloudSpec] = []
    for (fs,), group in _grouped(rows, "feature_set").items():
        for comp in components:
            words = [(_label(r["feature"], fs), _num(r.get(comp))) for r in group]
            words = [(t, w) for t, w in words if w]
            words.sort(key=lambda tw: -abs(tw[1]))
            chosen = words[:max(0, int(component_words))]
            if not chosen:
                continue        # this set has fewer components than the widest one
            out.append(CloudSpec(
                analysis=heading, folder=f"{folder}/{slug(fs)}/components",
                set_name=fs, kind="component", name=slug(comp),
                title=f"{heading}: {comp} of {fs} — what loads on it",
                legend="Bigger and darker = larger loading. Blue loads "
                       "positively, red negatively.",
                words=chosen))
    return out

clouds_for_correlations

clouds_for_correlations(rows, *, max_p, method='')

Two clouds per outcome per feature set from a correlations table: features correlated positively and negatively, sized by r, only those under max_p (adjusted when the table has adjusted p-values).

Source code in src\taters\figures\wordclouds.py

def clouds_for_correlations(rows: Sequence[Row], *, max_p: float,
                            method: str = "") -> List[CloudSpec]:
    """Two clouds per outcome per feature set from a correlations table:
    features correlated positively and negatively, sized by r, only those
    under ``max_p`` (adjusted when the table has adjusted p-values)."""
    if not rows:
        return []
    heading, folder = ANALYSES["correlations"]
    split = _split_col(rows)
    cols = list(rows[0].keys())
    outcomes = [c[:-2] for c in cols
                if c.endswith("_r") and f"{c[:-2]}_p" in cols]
    adjusted = any(c.endswith("_p_adj") for c in cols)
    which = f"{method} " if method else ""
    adj = " (adjusted)" if adjusted else ""
    out: List[CloudSpec] = []
    for (sv, fs), group in _grouped(rows, split or "", "feature_set").items():
        for outcome in outcomes:
            words = [(_label(r["feature"], fs), _num(r.get(f"{outcome}_r")))
                     for r in group if _p_of(r, outcome, max_p)]
            pos, neg = _signed([(t, w) for t, w in words if w])
            where = _where(split, sv)
            legend = (f"Bigger and darker = larger |r|; only {which}correlations "
                      f"with p < {max_p:g}{adj} are shown.")
            out += _pair(
                heading, folder, fs, f"{slug(outcome)}{_tag(split, sv)}", pos, neg,
                titles=(f"Correlations: {outcome} — features correlated "
                        f"positively ({fs}){where}",
                        f"Correlations: {outcome} — features correlated "
                        f"negatively ({fs}){where}"),
                legends=(legend, legend),
                notes=(f"no feature was positively correlated with {outcome} at "
                       f"p < {max_p:g}{adj} ({fs}).",
                       f"no feature was negatively correlated with {outcome} at "
                       f"p < {max_p:g}{adj} ({fs})."),
                stat=f"{outcome}: r")
    return out

clouds_for_frequencies

clouds_for_frequencies(rows, *, top_words)

The corpus cloud from a frequency list: the top_words most frequent n-grams after :func:skip_redundant, sized by frequency.

Source code in src\taters\figures\wordclouds.py

def clouds_for_frequencies(rows: Sequence[Row], *, top_words: int) -> List[CloudSpec]:
    """The corpus cloud from a frequency list: the ``top_words`` most
    frequent n-grams after :func:`skip_redundant`, sized by frequency."""
    if not rows:
        return []
    phrases = []
    for r in rows:
        f = _num(r.get("frequency"))
        text = (r.get("ngram") or "").strip()
        if not text or not f:
            continue
        phrases.append((text, f))
    chosen = skip_redundant(phrases)[:max(0, int(top_words))]
    if rows and rows[0].get("pos") is not None:
        tags = {r.get("ngram", ""): r.get("pos", "") for r in rows}
        chosen = [(f"{t} ({tags[t]})" if tags.get(t) else t, f) for t, f in chosen]
    return [CloudSpec(
        analysis="Most frequent terms",
        name="top_words",
        title=f"The {len(chosen)} most frequent terms in the corpus",
        legend="Bigger and darker = more frequent. Phrases that only repeat "
               "a more frequent word are left out.",
        words=chosen,
        note="" if chosen else "the frequency list is empty.")]

clouds_for_frequencies_by_group

clouds_for_frequencies_by_group(
    table_rows, sets, *, group_col, max_words
)

One cloud per group level per document-term set: the terms used most in that group's texts, sized by their summed count in the analysis table.

The analysis table already joins the matrix and the grouping column, so this is a column sum per level and nothing is re-tokenized.

Source code in src\taters\figures\wordclouds.py

def clouds_for_frequencies_by_group(table_rows: Sequence[Row],
                                    sets: Dict[str, List[str]], *,
                                    group_col: str, max_words: int
                                    ) -> List[CloudSpec]:
    """
    One cloud per group level per document-term set: the terms used most in
    that group's texts, sized by their summed count in the analysis table.

    The analysis table already joins the matrix and the grouping column, so
    this is a column sum per level and nothing is re-tokenized.
    """
    if not table_rows or not group_col or group_col not in table_rows[0]:
        return []
    heading, folder = _FREQ_BY_GROUP
    out: List[CloudSpec] = []
    for fs, cols in sets.items():
        if not fs.startswith("doc_term_matrix"):
            continue
        totals: Dict[str, Dict[str, float]] = {}
        for r in table_rows:
            level = str(r.get(group_col) or "").strip()
            if not level:
                continue
            bucket = totals.setdefault(level, {})
            for c in cols:
                v = _num(r.get(c))
                if v:
                    bucket[c] = bucket.get(c, 0.0) + v
        for level in sorted(totals):
            words = sorted(((_label(c, fs), v) for c, v in totals[level].items()),
                           key=lambda tw: -tw[1])[:max(0, int(max_words))]
            if not words:
                continue
            out.append(CloudSpec(
                analysis=heading, folder=f"{folder}/{slug(fs)}", set_name=fs,
                name=f"{slug(group_col)}-{slug(level)}",
                title=f"Most frequent terms — {group_col} = {level} ({fs})",
                legend=f"Bigger and darker = used more in the texts of "
                       f"{group_col} = {level}.",
                words=words))
    return out

clouds_for_group_differences

clouds_for_group_differences(rows, *, max_p)

Two clouds per pair of groups per feature set from group_differences_pairwise.csv: what is higher in the first group and what is higher in the second, sized by Cohen's d.

Source code in src\taters\figures\wordclouds.py

def clouds_for_group_differences(rows: Sequence[Row], *, max_p: float
                                 ) -> List[CloudSpec]:
    """Two clouds per pair of groups per feature set from
    ``group_differences_pairwise.csv``: what is higher in the first group
    and what is higher in the second, sized by Cohen's d."""
    if not rows:
        return []
    heading, folder = ANALYSES["group_differences"]
    split = _split_col(rows)
    cols = list(rows[0].keys())
    pcol = "p_adj" if "p_adj" in cols else "p"
    shown_p = pcol.replace("_", " ")
    out: List[CloudSpec] = []
    keyed = _grouped(rows, split or "", "feature_set", "group_1", "group_2")
    for (sv, fs, g1, g2), group in keyed.items():
        words = []
        for r in group:
            p = _num(r.get(pcol))
            d = _num(r.get("d"))
            if p is not None and p < max_p and d:
                words.append((_label(r["feature"], fs), d))
        pos, neg = _signed(words)
        base = f"{slug(g1)}_vs_{slug(g2)}{_tag(split, sv)}"
        where = _where(split, sv)
        for chosen, higher, other, color in ((pos, g1, g2, "Blue"),
                                              (neg, g2, g1, "Red")):
            out.append(CloudSpec(
                analysis=heading, folder=f"{folder}/{slug(fs)}", set_name=fs,
                name=f"{base}__higher_in_{slug(higher)}",
                title=f"Group differences: {g1} vs {g2} — higher in {higher} "
                      f"({fs}){where}",
                legend=f"Bigger and darker = larger Cohen's d; only pairs "
                       f"with {shown_p} < {max_p:g} are shown. "
                       f"{color}: higher in {higher}.",
                words=chosen,
                note="" if chosen else
                f"no feature was higher in {higher} than in {other} at "
                f"{shown_p} < {max_p:g} ({fs}).",
                stat=f"{g1} vs {g2}: d"))
    return out

clouds_for_neighbors

clouds_for_neighbors(rows, *, top_words)

One cloud per probe from a word-vector model's neighbors table (probe, rank, word, similarity), sized by similarity. A probe the model does not know has a row with blanks, and gets a sentence.

Source code in src\taters\figures\wordclouds.py

def clouds_for_neighbors(rows: Sequence[Row], *, top_words: int) -> List[CloudSpec]:
    """One cloud per probe from a word-vector model's neighbors table
    (``probe, rank, word, similarity``), sized by similarity. A probe the
    model does not know has a row with blanks, and gets a sentence."""
    if not rows:
        return []
    by_probe: Dict[str, List[Word]] = {}
    for r in rows:
        probe = (r.get("probe") or "").strip()
        if not probe:
            continue
        by_probe.setdefault(probe, [])
        sim = _num(r.get("similarity"))
        word = (r.get("word") or "").strip()
        if word and sim is not None:
            by_probe[probe].append((word, sim))
    out: List[CloudSpec] = []
    for probe, words in by_probe.items():
        words.sort(key=lambda tw: -tw[1])
        chosen = [(w, s) for w, s in words if s > 0][:max(0, int(top_words))]
        out.append(CloudSpec(
            analysis="Nearest neighbors",
            name=slug(probe),
            title=f"Words closest to “{probe}” in the model",
            legend="Bigger and darker = more similar (cosine).",
            words=chosen,
            note="" if chosen else f"“{probe}” is not in the model's vocabulary."))
    return out

clouds_for_ridge

clouds_for_ridge(rows, allowed)

Two clouds per outcome per feature set from ridge_coefficients.csv: the predictors of a higher score and of a lower one, sized by their standardized coefficient in the full model.

Source code in src\taters\figures\wordclouds.py

def clouds_for_ridge(rows: Sequence[Row], allowed: Callable[[str, str], bool]
                     ) -> List[CloudSpec]:
    """
    Two clouds per outcome per feature set from ``ridge_coefficients.csv``:
    the predictors of a higher score and of a lower one, sized by their
    standardized coefficient in the full model.
    """
    if not rows:
        return []
    heading, folder = ANALYSES["ridge"]
    split = _split_col(rows)
    cols = list(rows[0].keys())
    outcomes = cols[cols.index("predictor") + 1:]
    out: List[CloudSpec] = []
    for (sv, fs), group in _grouped(rows, split or "", "feature_set").items():
        if "+" in fs:
            continue        # a combination of tables; its members get the clouds
        model_rows = _full_model(group)
        for outcome in outcomes:
            words = [(_label(r["predictor"], fs), _num(r.get(outcome)))
                     for r in model_rows if allowed(fs, r["predictor"])]
            pos, neg = _signed([(t, w) for t, w in words if w])
            where = _where(split, sv)
            out += _pair(
                heading, folder, fs, f"{slug(outcome)}{_tag(split, sv)}", pos, neg,
                titles=(f"Ridge: {outcome} — features predicting higher scores "
                        f"({fs}){where}",
                        f"Ridge: {outcome} — features predicting lower scores "
                        f"({fs}){where}"),
                legends=(f"Bigger and darker = larger standardized coefficient. "
                         f"Blue: predicts a higher {outcome}.",
                         f"Bigger and darker = larger standardized coefficient. "
                         f"Red: predicts a lower {outcome}."),
                notes=(f"no feature had a positive coefficient for {outcome} ({fs}).",
                       f"no feature had a negative coefficient for {outcome} ({fs})."),
                stat=f"{outcome}: β")
    return out

clouds_for_theme_predictors

clouds_for_theme_predictors(
    loadings,
    mentions,
    *,
    analysis,
    folder,
    set_name,
    top_words
)

What each theme that appears in a set's clouds is made of.

A ridge over topic-model themes draws "Theme_5" as a word, which says nothing until the theme's own words are in view. For every theme named in mentions (theme -> where it appeared, as (statistic, value) pairs such as ("TIPI_Open: β", 0.32)) that has a column in the loadings table, one mixed-sign cloud of its top_words strongest terms, filed under the set's themes/ folder, its legend naming the three places it mattered most.

Source code in src\taters\figures\wordclouds.py

def clouds_for_theme_predictors(loadings: Sequence[Row],
                                mentions: Dict[str, List[Tuple[str, float]]], *,
                                analysis: str, folder: str, set_name: str,
                                top_words: int) -> List[CloudSpec]:
    """
    What each theme that appears in a set's clouds is made of.

    A ridge over topic-model themes draws "Theme_5" as a word, which says
    nothing until the theme's own words are in view. For every theme named
    in ``mentions`` (theme -> where it appeared, as (statistic, value) pairs
    such as ``("TIPI_Open: β", 0.32)``) that has a column in the loadings
    table, one mixed-sign cloud of its ``top_words`` strongest terms, filed
    under the set's ``themes/`` folder, its legend naming the three places
    it mattered most.
    """
    if not loadings or not mentions:
        return []
    out: List[CloudSpec] = []
    for theme in sorted(mentions, key=lambda t: (len(t), t)):
        # a name with no column in the loadings (a control, a component)
        # just gathers no words below, and we skip it there.
        words = []
        for r in loadings:
            v = _num(r.get(theme))
            if not v:
                continue
            label = r.get("term", "")
            if r.get("pos"):
                label = f"{label} ({r['pos']})"
            words.append((label, v))
        words.sort(key=lambda tw: -abs(tw[1]))
        chosen = words[:max(0, int(top_words))]
        if not chosen:
            continue
        seen = sorted(mentions[theme], key=lambda sv: -abs(sv[1]))
        said = "; ".join(f"{stat} {v:+.2f}" for stat, v in seen[:3])
        out.append(CloudSpec(
            analysis=analysis, folder=f"{folder}/{slug(set_name)}/themes",
            set_name=set_name, kind="theme", name=slug(theme),
            title=f"{theme.replace('_', ' ')} — what it is made of ({set_name})",
            legend="Bigger and darker = larger loading; blue loads positively, "
                   "red negatively. Where it mattered most here: " + said
                   + (f" (and {len(seen) - 3} more)" if len(seen) > 3 else "") + ".",
            words=chosen))
    return out

clouds_for_themes

clouds_for_themes(
    rows, *, top_words, min_abs_loading=0.0, shares=None
)

One mixed-sign cloud per theme from a topic model's loadings table (term, [pos,] Theme_1…). A pos column joins the label, so felt (VBD) reads as the tagged term it is.

Source code in src\taters\figures\wordclouds.py

def clouds_for_themes(rows: Sequence[Row], *, top_words: int,
                      min_abs_loading: float = 0.0,
                      shares: Optional[Dict[str, float]] = None) -> List[CloudSpec]:
    """One mixed-sign cloud per theme from a topic model's loadings table
    (``term, [pos,] Theme_1…``). A ``pos`` column joins the label, so
    ``felt (VBD)`` reads as the tagged term it is."""
    if not rows:
        return []
    cols = list(rows[0].keys())
    themes = [c for c in cols if c not in ("term", "pos")]
    out: List[CloudSpec] = []
    for theme in themes:
        words = []
        for r in rows:
            v = _num(r.get(theme))
            if v is None or abs(v) < min_abs_loading or v == 0:
                continue
            label = r.get("term", "")
            if r.get("pos"):
                label = f"{label} ({r['pos']})"
            words.append((label, v))
        words.sort(key=lambda tw: -abs(tw[1]))
        chosen = words[:max(0, int(top_words))]
        share = (shares or {}).get(theme)
        out.append(CloudSpec(
            analysis="Themes",
            name=slug(theme),
            title=f"{theme.replace('_', ' ')} — the words that define it"
                  + (f" ({share:.1f}% of variance)" if share is not None else ""),
            legend="Bigger and darker = larger loading. Blue loads "
                   "positively, red negatively.",
            words=chosen,
            note="" if chosen else f"{theme} has no term loading above "
                                   f"{min_abs_loading:g}."))
    return out

frequency_wordclouds

frequency_wordclouds(
    freq_csv,
    out_dir=None,
    *,
    top_words=100,
    enabled=True,
    overwrite_existing=False,
    encoding="utf-8-sig",
    verbose=True,
    on_progress=None
)

The corpus's most frequent words and phrases as one cloud.

Parameters:

Name	Type	Description	Default
freq_csv	str or Path	The n-gram frequency list (ngram, frequency, …).	required
out_dir	str or Path	Where the picture goes. Defaults to <features folder>/figures/wordclouds/ngram_frequencies.	None
top_words	int	How many terms, by frequency, after phrases that only repeat a more frequent word are left out.	100
enabled	bool	As for :func:stats_wordclouds.	True
overwrite_existing	bool	As for :func:stats_wordclouds.	True
encoding	bool	As for :func:stats_wordclouds.	True
verbose	bool	As for :func:stats_wordclouds.	True
on_progress	bool	As for :func:stats_wordclouds.	True

Returns:

Type	Description
Path	The output folder holding top_words.png.

Source code in src\taters\figures\wordclouds.py

def frequency_wordclouds(
    freq_csv: PathLike,
    out_dir: Optional[PathLike] = None,
    *,
    top_words: int = 100,
    enabled: bool = True,
    overwrite_existing: bool = False,
    encoding: str = "utf-8-sig",
    verbose: bool = True,
    on_progress: Optional[Callable[..., None]] = None,
) -> Path:
    """
    The corpus's most frequent words and phrases as one cloud.

    Parameters
    ----------
    freq_csv : str or Path
        The n-gram frequency list (``ngram``, ``frequency``, …).
    out_dir : str or Path, optional
        Where the picture goes. Defaults to
        ``<features folder>/figures/wordclouds/ngram_frequencies``.
    top_words : int, default 100
        How many terms, by frequency, after phrases that only repeat a more
        frequent word are left out.
    enabled, overwrite_existing, encoding, verbose, on_progress
        As for :func:`stats_wordclouds`.

    Returns
    -------
    Path
        The output folder holding ``top_words.png``.
    """
    freq_csv = Path(freq_csv)
    folder = Path(out_dir) if out_dir else (
        freq_csv.parent / FIGURES_DIR / "ngram_frequencies")
    if not enabled:
        if verbose:
            print("[figures] the frequency word cloud is off for this run.")
        return folder
    missing = pillow_missing_reason()
    if missing:
        if verbose:
            print(f"[figures] {missing}")
        return folder
    specs = clouds_for_frequencies(_rows(freq_csv, encoding), top_words=top_words)
    for s in specs:
        s.inputs = [freq_csv]
    drawn, reused = _draw_all(specs, folder, max_words=top_words,
                              overwrite_existing=overwrite_existing,
                              verbose=verbose, on_progress=on_progress)
    if verbose:
        print(f"[figures] {drawn} frequency cloud(s) drawn, {reused} already "
              f"current -> {folder}")
    return folder

neighbor_wordclouds

neighbor_wordclouds(
    neighbors_csv,
    out_dir=None,
    *,
    top_words=20,
    enabled=True,
    overwrite_existing=False,
    encoding="utf-8-sig",
    verbose=True,
    on_progress=None
)

One word cloud per probe of a word-vector model's neighbors table.

Parameters:

Name	Type	Description	Default
neighbors_csv	str or Path	The *_neighbors.csv the word-vector steps write: probe, rank, word, similarity.	required
out_dir	str or Path	Where the pictures go. Defaults to <features folder>/figures/wordclouds/word_vectors.	None
top_words	int	How many neighbors per probe, by similarity.	20
enabled	bool	As for :func:stats_wordclouds.	True
overwrite_existing	bool	As for :func:stats_wordclouds.	True
encoding	bool	As for :func:stats_wordclouds.	True
verbose	bool	As for :func:stats_wordclouds.	True
on_progress	bool	As for :func:stats_wordclouds.	True

Returns:

Type	Description
Path	The output folder; it also holds an index.md listing the probes.

Source code in src\taters\figures\wordclouds.py

def neighbor_wordclouds(
    neighbors_csv: PathLike,
    out_dir: Optional[PathLike] = None,
    *,
    top_words: int = 20,
    enabled: bool = True,
    overwrite_existing: bool = False,
    encoding: str = "utf-8-sig",
    verbose: bool = True,
    on_progress: Optional[Callable[..., None]] = None,
) -> Path:
    """
    One word cloud per probe of a word-vector model's neighbors table.

    Parameters
    ----------
    neighbors_csv : str or Path
        The ``*_neighbors.csv`` the word-vector steps write:
        ``probe, rank, word, similarity``.
    out_dir : str or Path, optional
        Where the pictures go. Defaults to
        ``<features folder>/figures/wordclouds/word_vectors``.
    top_words : int, default 20
        How many neighbors per probe, by similarity.
    enabled, overwrite_existing, encoding, verbose, on_progress
        As for :func:`stats_wordclouds`.

    Returns
    -------
    Path
        The output folder; it also holds an ``index.md`` listing the probes.
    """
    neighbors_csv = Path(neighbors_csv)
    folder = Path(out_dir) if out_dir else (
        neighbors_csv.parent.parent / FIGURES_DIR / "word_vectors")
    if not enabled:
        if verbose:
            print("[figures] neighbor word clouds are off for this run.")
        return folder
    missing = pillow_missing_reason()
    if missing:
        if verbose:
            print(f"[figures] {missing}")
        return folder
    specs = clouds_for_neighbors(_rows(neighbors_csv, encoding), top_words=top_words)
    for s in specs:
        s.inputs = [neighbors_csv]
    drawn, reused = _draw_all(specs, folder, max_words=top_words,
                              overwrite_existing=overwrite_existing,
                              verbose=verbose, on_progress=on_progress)
    lines = ["# Nearest neighbors", "",
             f"Drawn from `{neighbors_csv.name}`: the {top_words} words closest "
             "to each probe by cosine similarity. Bigger and darker is more "
             "similar. These are the evidence that the model learned what the "
             "study assumes it learned.", ""]
    for s in specs:
        if s.words:
            lines += [f"**{s.title}**  ", f"![{s.title}]({s.name}.png)", ""]
        elif s.note:
            lines += [f"{s.note}", ""]
    _write_if_changed(folder / "index.md", "\n".join(lines))
    if verbose:
        print(f"[figures] {drawn} neighbor cloud(s) drawn, {reused} already "
              f"current -> {folder}")
    return folder

skip_redundant

skip_redundant(phrases)

Drop an n-gram that only repeats words the cloud already shows.

Walking down by frequency, a phrase every one of whose words is already in the cloud on its own -- "of the" beside "of" and "the" -- says nothing new and is skipped; a phrase that brings at least one new word stays ("going to be" when "going" is not there yet). Single words are always kept. A frequency-based test was tried first and kept "of the" because "of" is rarer than "the", which is true and not the point.

Source code in src\taters\figures\wordclouds.py

def skip_redundant(phrases: Sequence[Tuple[str, float]]) -> List[Tuple[str, float]]:
    """
    Drop an n-gram that only repeats words the cloud already shows.

    Walking down by frequency, a phrase every one of whose words is already
    in the cloud on its own -- "of the" beside "of" and "the" -- says nothing
    new and is skipped; a phrase that brings at least one new word stays
    ("going to be" when "going" is not there yet). Single words are always
    kept. A frequency-based test was tried first and kept "of the" because
    "of" is rarer than "the", which is true and not the point.
    """
    kept: List[Tuple[str, float]] = []
    shown: set = set()
    for text, weight in sorted(phrases, key=lambda tw: -tw[1]):
        words = text.split()
        if len(words) > 1 and all(w in shown for w in words):
            continue
        kept.append((text, weight))
        shown.update(words)
    return kept

stats_wordclouds

stats_wordclouds(
    stats_dir="stats_results",
    *,
    max_words=80,
    max_p=0.05,
    component_words=30,
    group_col="",
    enabled=True,
    overwrite_existing=False,
    encoding="utf-8-sig",
    verbose=True,
    on_progress=None
)

Draw word clouds of every statistics result in a folder, and a report section showing them.

Parameters:

Name	Type	Description	Default
stats_dir	str or Path	The statistics output folder. Whatever tables it holds are drawn: ridge and classifier coefficients, correlations, pairwise group differences, PCA loadings, and per-group term frequencies when the analysis table joins a document-term matrix and group_col names the grouping column.	"stats_results"
max_words	int	The most words in any one cloud.	80
max_p	float	Correlations and pairwise differences are drawn only when they pass this, on the adjusted p-value where the table has one.	0.05
component_words	int	How many features to show for each principal component, and how many terms for each topic-model theme a set's clouds name.	30
group_col	str	The metadata column whose levels get a frequency cloud each. Empty draws none.	""
enabled	bool	Off draws nothing and writes nothing -- the switch behind wordclouds: false in a pipeline.	True
overwrite_existing	bool	Redraw pictures that are newer than their tables anyway.	False
encoding	str	The tables' encoding.	"utf-8-sig"
verbose	bool	The usual.	True
on_progress	bool	The usual.	True

Returns:

Type	Description
Path	<stats_dir>/figures/wordclouds. Under it, one folder per analysis and one per feature set inside it -- ridge-regression/dictionary/ -- with components/ and themes/ folders for what that set's predictors are made of.

Notes

Coefficient clouds come from one model per outcome: language fitted with every control when there were controls, language alone otherwise. Controls themselves are not drawn -- age as the biggest word in a cloud about language is true and beside the point. A cloud that would be empty is not drawn; the report says so instead.

Source code in src\taters\figures\wordclouds.py

def stats_wordclouds(
    stats_dir: PathLike = "stats_results",
    *,
    max_words: int = 80,
    max_p: float = 0.05,
    component_words: int = 30,
    group_col: str = "",
    enabled: bool = True,
    overwrite_existing: bool = False,
    encoding: str = "utf-8-sig",
    verbose: bool = True,
    on_progress: Optional[Callable[..., None]] = None,
) -> Path:
    """
    Draw word clouds of every statistics result in a folder, and a report
    section showing them.

    Parameters
    ----------
    stats_dir : str or Path, default "stats_results"
        The statistics output folder. Whatever tables it holds are drawn:
        ridge and classifier coefficients, correlations, pairwise group
        differences, PCA loadings, and per-group term frequencies when the
        analysis table joins a document-term matrix and ``group_col`` names
        the grouping column.
    max_words : int, default 80
        The most words in any one cloud.
    max_p : float, default 0.05
        Correlations and pairwise differences are drawn only when they pass
        this, on the adjusted p-value where the table has one.
    component_words : int, default 30
        How many features to show for each principal component, and how
        many terms for each topic-model theme a set's clouds name.
    group_col : str, default ""
        The metadata column whose levels get a frequency cloud each. Empty
        draws none.
    enabled : bool, default True
        Off draws nothing and writes nothing -- the switch behind
        ``wordclouds: false`` in a pipeline.
    overwrite_existing : bool, default False
        Redraw pictures that are newer than their tables anyway.
    encoding : str, default "utf-8-sig"
        The tables' encoding.
    verbose, on_progress
        The usual.

    Returns
    -------
    Path
        ``<stats_dir>/figures/wordclouds``. Under it, one folder per analysis
        and one per feature set inside it -- ``ridge-regression/dictionary/``
        -- with ``components/`` and ``themes/`` folders for what that set's
        predictors are made of.

    Notes
    -----
    Coefficient clouds come from one model per outcome: language fitted
    with every control when there were controls, language alone otherwise.
    Controls themselves are not drawn -- ``age`` as the biggest word in a
    cloud about language is true and beside the point. A cloud that would be
    empty is not drawn; the report says so instead.
    """
    stats_dir = Path(stats_dir)
    root = stats_dir / FIGURES_DIR
    if not enabled:
        if verbose:
            print("[figures] word clouds are off for this run.")
        return root

    missing = pillow_missing_reason()
    if missing:
        if verbose:
            print(f"[figures] {missing}")
        _write_section_if_changed(stats_dir, _section_md([], skipped=missing))
        return root

    announce(on_progress, "choosing the words for each cloud")
    table_csv = stats_dir / "analysis_table.csv"
    sets, sources = _sidecar(table_csv) if table_csv.is_file() else (None, {})
    components = _components_by_set(stats_dir, encoding)
    allowed = _feature_filter(sets, components)

    specs: List[CloudSpec] = []

    def take(path: Path, made: List[CloudSpec]) -> None:
        for s in made:
            s.inputs = [path]
        specs.extend(made)

    p = stats_dir / "group_differences_pairwise.csv"
    if p.is_file():
        take(p, clouds_for_group_differences(_rows(p, encoding), max_p=max_p))
    for method in ("pearson", "spearman"):
        p = stats_dir / f"correlations_{method}.csv"
        if p.is_file():
            take(p, clouds_for_correlations(_rows(p, encoding), max_p=max_p,
                                            method=method))
    p = stats_dir / "ridge_coefficients.csv"
    if p.is_file():
        take(p, clouds_for_ridge(_rows(p, encoding), allowed))
    p = stats_dir / "classifier_coefficients.csv"
    if p.is_file():
        take(p, clouds_for_classifier(_rows(p, encoding), allowed))
    for p in sorted(stats_dir.glob("*_pca_loadings*.csv")):
        analysis = p.name.split("_pca_loadings")[0]
        take(p, clouds_for_components(_rows(p, encoding), analysis=analysis,
                                      component_words=component_words))
    if sets and group_col and table_csv.is_file():
        take(table_csv, clouds_for_frequencies_by_group(
            _rows(table_csv, encoding), sets, group_col=group_col,
            max_words=max_words))

    # lastly, what the themes a set's clouds name are made of, filed next to
    # them.
    if sets:
        mentions: Dict[Tuple[str, str, str], Dict[str, List[Tuple[str, float]]]] = {}
        for s in specs:
            if s.kind != "cloud" or not s.words or not s.set_name:
                continue
            folder_root = s.folder.rsplit("/", 1)[0]
            bucket = mentions.setdefault((s.analysis, folder_root, s.set_name), {})
            for label, w in s.words:
                bucket.setdefault(label, []).append((s.stat, w))
        loadings_cache: Dict[str, Optional[Path]] = {}
        for (analysis, folder_root, fs), named in mentions.items():
            if fs not in loadings_cache:
                loadings_cache[fs] = _theme_loadings_for(
                    fs, sets.get(fs, []), sources, stats_dir, encoding)
            path = loadings_cache[fs]
            if path is None:
                continue
            made = clouds_for_theme_predictors(
                _rows(path, encoding), named, analysis=analysis,
                folder=folder_root, set_name=fs, top_words=component_words)
            for s in made:
                s.inputs = [path] + [t for t in (stats_dir / "analysis_table.csv",)
                                     if t.is_file()]
            specs.extend(made)

    drawn, reused = _draw_all(specs, root, max_words=max_words,
                              overwrite_existing=overwrite_existing,
                              verbose=verbose, on_progress=on_progress)
    if specs:
        _write_section_if_changed(stats_dir, _section_md(specs))
    if verbose:
        print(f"[figures] {drawn} word cloud(s) drawn, {reused} already current"
              f" -> {root}" if specs else
              "[figures] no statistics tables to draw word clouds from.")
    return root

theme_wordclouds

theme_wordclouds(
    loadings_csv,
    out_dir=None,
    *,
    top_words=30,
    min_abs_loading=0.0,
    enabled=True,
    overwrite_existing=False,
    encoding="utf-8-sig",
    verbose=True,
    on_progress=None
)

One word cloud per theme of a topic model, from its loadings table.

Parameters:

Name	Type	Description	Default
loadings_csv	str or Path	The *_loadings.csv the topic-model step writes: term, an optional pos, then one column per theme.	required
out_dir	str or Path	Where the pictures go. Defaults to <features folder>/figures/wordclouds/topic_model_mem.	None
top_words	int	How many terms per theme, by the size of their loading.	30
min_abs_loading	float	Leave out terms loading less than this in either direction.	0.0
enabled	bool	As for :func:stats_wordclouds.	True
overwrite_existing	bool	As for :func:stats_wordclouds.	True
encoding	bool	As for :func:stats_wordclouds.	True
verbose	bool	As for :func:stats_wordclouds.	True
on_progress	bool	As for :func:stats_wordclouds.	True

Returns:

Type	Description
Path	The output folder; it also holds an index.md listing the themes.

Source code in src\taters\figures\wordclouds.py

def theme_wordclouds(
    loadings_csv: PathLike,
    out_dir: Optional[PathLike] = None,
    *,
    top_words: int = 30,
    min_abs_loading: float = 0.0,
    enabled: bool = True,
    overwrite_existing: bool = False,
    encoding: str = "utf-8-sig",
    verbose: bool = True,
    on_progress: Optional[Callable[..., None]] = None,
) -> Path:
    """
    One word cloud per theme of a topic model, from its loadings table.

    Parameters
    ----------
    loadings_csv : str or Path
        The ``*_loadings.csv`` the topic-model step writes: ``term``, an
        optional ``pos``, then one column per theme.
    out_dir : str or Path, optional
        Where the pictures go. Defaults to
        ``<features folder>/figures/wordclouds/topic_model_mem``.
    top_words : int, default 30
        How many terms per theme, by the size of their loading.
    min_abs_loading : float, default 0.0
        Leave out terms loading less than this in either direction.
    enabled, overwrite_existing, encoding, verbose, on_progress
        As for :func:`stats_wordclouds`.

    Returns
    -------
    Path
        The output folder; it also holds an ``index.md`` listing the themes.
    """
    loadings_csv = Path(loadings_csv)
    folder = Path(out_dir) if out_dir else (
        loadings_csv.parent / FIGURES_DIR / "topic_model_mem")
    if not enabled:
        if verbose:
            print("[figures] theme word clouds are off for this run.")
        return folder
    missing = pillow_missing_reason()
    if missing:
        if verbose:
            print(f"[figures] {missing}")
        return folder

    shares: Dict[str, float] = {}
    eigen = loadings_csv.with_name(
        loadings_csv.name.replace("_loadings", "_eigenvalues"))
    if eigen != loadings_csv and eigen.is_file():
        for r in _rows(eigen, encoding):
            v = _num(r.get("pct_variance"))
            if v is not None:
                shares[r.get("theme", "")] = v
    specs = clouds_for_themes(_rows(loadings_csv, encoding), top_words=top_words,
                              min_abs_loading=min_abs_loading, shares=shares)
    for s in specs:
        s.inputs = [loadings_csv]
    drawn, reused = _draw_all(specs, folder, max_words=top_words,
                              overwrite_existing=overwrite_existing,
                              verbose=verbose, on_progress=on_progress)
    lines = ["# Themes", "", f"Drawn from `{loadings_csv.name}`: the {top_words} "
             "terms loading most strongly on each theme. Blue loads positively, "
             "red negatively; bigger and darker is a larger loading.", ""]
    for s in specs:
        if s.words:
            lines += [f"**{s.title}**  ", f"![{s.title}]({s.name}.png)", ""]
        elif s.note:
            lines += [f"{s.note}", ""]
    _write_if_changed(folder / "index.md", "\n".join(lines))
    if verbose:
        print(f"[figures] {drawn} theme cloud(s) drawn, {reused} already current "
              f"-> {folder}")
    return folder

taters.figures.render

Draw a word cloud: words sized by a weight, colored by its sign and strength.

The layout is ours -- a spiral outward from the center, horizontal words only, no two boxes overlapping -- and it is pure Python, so it is tested without a font engine by handing it a fake measurer. Pillow is used for the two things that need a font: measuring a word at a size, and drawing the PNG. It is imported inside the functions that need it, so importing this module costs nothing and a machine without Pillow gets one plain sentence instead of a traceback.

What the picture means is fixed here, once, for every cloud in Taters: size follows the magnitude of the weight, shade darkens with it, blue is a positive weight and red a negative one. A caller that wants a cloud of one sign filters its words first; a caller with signed loadings passes them as they are and gets both colors in one picture.

Layout dataclass

Layout(placed, dropped, width, height, top_margin)

Where every word landed, and which did not fit at any legible size.

Placement dataclass

Placement(text, weight, size, x, y, w, h)

One word on the canvas: where its ink goes and how big it is.

box property

box

(x0, y0, x1, y1) of the ink, canvas pixels.

color_for

color_for(weight, max_abs)

Blue for a positive weight, red for a negative one, darker the larger it is relative to the cloud's strongest word.

Source code in src\taters\figures\render.py

def color_for(weight: float, max_abs: float) -> RGB:
    """Blue for a positive weight, red for a negative one, darker the larger
    it is relative to the cloud's strongest word."""
    strength = abs(weight) / max_abs if max_abs > 0 else 1.0
    return shade(BLUE_HUE if weight > 0 else RED_HUE, strength)

fit_budget

fit_budget(
    words,
    sizes,
    measure,
    *,
    width,
    height,
    top_margin=0,
    min_font=12,
    fill=FILL
)

Sizes scaled so the words can actually fit on the canvas.

:func:scale_sizes knows the weights and nothing about the canvas; eighty words at up to 96 px need more than 1200 by 800 pixels, and the layout's only answers are shrinking one word at a time or dropping it (36 of 80 dropped, on the first real picture). Measuring the ink at the proposed sizes and scaling every size by the same factor keeps the picture's proportions -- the strongest word is still that much larger -- while making room for all of them.

The same rule grows a sparse cloud: five significant features at the weight-scaled sizes sat in the middle of an otherwise empty picture, so when the ink falls short of the budget every size is scaled up, until the largest word is a quarter of the drawing height. Never below min_font.

Source code in src\taters\figures\render.py

def fit_budget(words: Sequence[Word], sizes: Sequence[int], measure: Measure,
               *, width: int, height: int, top_margin: int = 0,
               min_font: int = 12, fill: float = FILL) -> List[int]:
    """
    Sizes scaled so the words can actually fit on the canvas.

    :func:`scale_sizes` knows the weights and nothing about the canvas; eighty
    words at up to 96 px need more than 1200 by 800 pixels, and the layout's
    only answers are shrinking one word at a time or dropping it (36 of 80
    dropped, on the first real picture). Measuring the ink at the proposed
    sizes and scaling every size by the same factor keeps the picture's
    proportions -- the strongest word is still that much larger -- while
    making room for all of them.

    The same rule grows a sparse cloud: five significant features at the
    weight-scaled sizes sat in the middle of an otherwise empty picture, so
    when the ink falls short of the budget every size is scaled up, until
    the largest word is a quarter of the drawing height. Never below
    ``min_font``.
    """
    if not words:
        return list(sizes)
    ink = 0.0
    for (text, _w), size in zip(words, sizes):
        w, h = measure(text, int(size))
        pad = _pad_of(int(size))
        ink += (w + 2 * pad) * (h + 2 * pad)
    room = fill * width * max(1, height - top_margin)
    if ink <= 0:
        return [int(s) for s in sizes]
    factor = math.sqrt(room / ink)
    if factor > 1.0:
        # when we grow, we aim lower than we'd allow when shrinking. this is
        # because a spiral packs big words less tightly than small ones, and
        # growing to the full budget made the tail shrink back down one word
        # at a time (three seconds a cloud). we also stop once the largest
        # word is a quarter of the drawing height, so that one word doesn't
        # turn into a banner.
        factor = math.sqrt(GROW_FILL * room / (fill * ink))
        cap = max(1, height - top_margin) // 4
        factor = min(factor, cap / max(max(sizes), 1))
        if factor <= 1.0:
            return [int(s) for s in sizes]
    return [max(int(min_font), int(round(s * factor))) for s in sizes]

layout

layout(
    words,
    sizes,
    measure,
    *,
    width,
    height,
    top_margin=0,
    min_font=12,
    shrink=0.85,
    seed=0,
    max_steps=6000
)

Place every word on a width by height canvas below top_margin.

Largest first, each on a spiral out from the center, so the strongest words hold the middle and the rest fill in around them. A word that finds no room at its size shrinks by shrink and tries again, down to min_font; below that it is dropped and named in Layout.dropped rather than drawn illegibly or on top of something. measure(text, size) returns the ink width and height of a word; the renderer passes a real font, the tests pass arithmetic. The only randomness is the angle each spiral starts at, from seed, so the same input draws the same picture every time.

Source code in src\taters\figures\render.py

def layout(words: Sequence[Word], sizes: Sequence[int], measure: Measure, *,
           width: int, height: int, top_margin: int = 0, min_font: int = 12,
           shrink: float = 0.85, seed: int = 0, max_steps: int = 6000) -> Layout:
    """
    Place every word on a ``width`` by ``height`` canvas below ``top_margin``.

    Largest first, each on a spiral out from the center, so the strongest
    words hold the middle and the rest fill in around them. A word that finds
    no room at its size shrinks by ``shrink`` and tries again, down to
    ``min_font``; below that it is dropped and named in ``Layout.dropped``
    rather than drawn illegibly or on top of something. ``measure(text,
    size)`` returns the ink width and height of a word; the renderer passes a
    real font, the tests pass arithmetic. The only randomness is the angle
    each spiral starts at, from ``seed``, so the same input draws the same
    picture every time.
    """
    if len(words) != len(sizes):
        raise ValueError("layout needs one size per word")
    rng = random.Random(seed)
    area = (0, int(top_margin), int(width), int(height))
    cx = width / 2.0
    cy = top_margin + (height - top_margin) / 2.0
    r_max = math.hypot(width, height - top_margin) / 2.0
    order = sorted(range(len(words)), key=lambda i: (-sizes[i], words[i][0]))
    placed: List[Placement] = []
    dropped: List[str] = []
    for i in order:
        text, weight = words[i]
        size = int(sizes[i])
        theta0 = rng.uniform(0.0, 2 * math.pi)
        while True:
            w, h = measure(text, size)
            pad = _pad_of(size)
            bw, bh = w + 2 * pad, h + 2 * pad
            spot = None
            if bw <= width and bh <= height - top_margin:
                spot = _spiral_search(bw, bh, cx, cy, r_max, theta0, area,
                                      placed, max_steps)
            if spot is not None:
                placed.append(Placement(text, weight, size,
                                        spot[0] + pad, spot[1] + pad, w, h))
                break
            if size <= min_font:
                dropped.append(text)
                break
            size = max(int(min_font), int(size * shrink))
    return Layout(_recentred(placed, area), dropped, int(width), int(height),
                  int(top_margin))

pillow_missing_reason

pillow_missing_reason()

None when Pillow can be imported, else the sentence to show instead of a figure.

Source code in src\taters\figures\render.py

def pillow_missing_reason() -> Optional[str]:
    """``None`` when Pillow can be imported, else the sentence to show instead
    of a figure."""
    try:
        _pillow()
    except ImportError as e:
        return str(e)
    return None

render_wordcloud

render_wordcloud(
    words,
    out_png,
    *,
    title,
    legend=None,
    width=1200,
    height=800,
    max_words=80,
    min_font=16,
    max_font=96,
    seed=0,
    font_path=None,
    on_progress=None,
    verbose=False
)

Draw one word cloud to a PNG.

Parameters:

Name	Type	Description	Default
words	sequence of (label, weight) or mapping	The words and their signed weights. Blue for positive, red for negative, size and shade by magnitude; the strongest max_words are drawn.	required
out_png	str or Path	Where the picture goes. Written atomically.	required
title	str	Drawn across the top of the picture, so it still says what it is once pasted into a slide without the report.	required
legend	str	The second caption line. Defaults to the sentence that explains size and color; pass "" for none.	None
width	int	Canvas size in pixels.	1200
height	int	Canvas size in pixels.	1200
max_words	int	How many words at most, and the smallest and largest type sizes.	80
min_font	int	How many words at most, and the smallest and largest type sizes.	80
max_font	int	How many words at most, and the smallest and largest type sizes.	80
seed	int	Fixes the layout, so the same words draw the same picture.	0
font_path	str or Path	A TrueType font to use instead of the bundled DejaVu Sans Bold.	None
on_progress	Optional[Callable[..., None]]	The usual: a phase announcement, and a line naming any words that did not fit.	None
verbose	Optional[Callable[..., None]]	The usual: a phase announcement, and a line naming any words that did not fit.	None

Returns:

Type	Description
Path	out_png.

Raises:

Type	Description
ImportError	When Pillow is not installed, with a sentence that says what to do. Callers that would rather skip the figure check :func:pillow_missing_reason first.

Source code in src\taters\figures\render.py

def render_wordcloud(words: Union[Sequence[Word], Mapping[str, float]],
                     out_png: Union[str, Path], *, title: str,
                     legend: Optional[str] = None, width: int = 1200,
                     height: int = 800, max_words: int = 80, min_font: int = 16,
                     max_font: int = 96, seed: int = 0,
                     font_path: Optional[Union[str, Path]] = None,
                     on_progress: Optional[Callable[..., None]] = None,
                     verbose: bool = False) -> Path:
    """
    Draw one word cloud to a PNG.

    Parameters
    ----------
    words : sequence of (label, weight) or mapping
        The words and their signed weights. Blue for positive, red for
        negative, size and shade by magnitude; the strongest ``max_words``
        are drawn.
    out_png : str or Path
        Where the picture goes. Written atomically.
    title : str
        Drawn across the top of the picture, so it still says what it is
        once pasted into a slide without the report.
    legend : str, optional
        The second caption line. Defaults to the sentence that explains size
        and color; pass ``""`` for none.
    width, height : int
        Canvas size in pixels.
    max_words, min_font, max_font : int
        How many words at most, and the smallest and largest type sizes.
    seed : int
        Fixes the layout, so the same words draw the same picture.
    font_path : str or Path, optional
        A TrueType font to use instead of the bundled DejaVu Sans Bold.
    on_progress, verbose
        The usual: a phase announcement, and a line naming any words that
        did not fit.

    Returns
    -------
    Path
        ``out_png``.

    Raises
    ------
    ImportError
        When Pillow is not installed, with a sentence that says what to do.
        Callers that would rather skip the figure check
        :func:`pillow_missing_reason` first.
    """
    Image, ImageDraw, ImageFont = _pillow()
    out_png = Path(out_png)
    announce(on_progress, f"drawing word cloud: {title}")

    fonts = _Fonts(ImageFont, Path(font_path) if font_path else DEFAULT_FONT)
    pairs = _as_pairs(words, max_words)
    legend = DEFAULT_LEGEND if legend is None else legend

    img = Image.new("RGB", (int(width), int(height)), (255, 255, 255))
    draw = ImageDraw.Draw(img)

    # the caption band comes first, because the cloud gets laid out below it.
    margin = 16
    band = margin
    caption = []
    for text, size, color, max_lines in ((title, 22, (40, 40, 40), 1),
                                          (legend, 15, (110, 110, 110), 3)):
        if not text:
            continue
        font = fonts[size]
        for line in _wrap(font, str(text), width - 2 * margin, max_lines):
            left, top, right, bottom = _ink(font, line)
            caption.append((line, font, color, margin - left, band - top))
            band += (bottom - top) + 6
    band += margin // 2 if caption else 0

    def measure(text: str, size: int) -> Tuple[int, int]:
        left, top, right, bottom = _ink(fonts[size], text)
        return (right - left, bottom - top)

    sizes = scale_sizes([w for _, w in pairs], min_font, max_font)
    sizes = fit_budget(pairs, sizes, measure, width=width, height=height,
                       top_margin=band, min_font=min_font)
    lay = layout(pairs, sizes, measure, width=width, height=height,
                 top_margin=band, min_font=min_font, seed=seed)

    for text, font, color, x, y in caption:
        draw.text((x, y), text, font=font, fill=color)
    max_abs = max((abs(w) for _, w in pairs), default=1.0)
    for p in lay.placed:
        font = fonts[p.size]
        left, top, _r, _b = _ink(font, p.text)
        draw.text((p.x - left, p.y - top), p.text, font=font,
                  fill=color_for(p.weight, max_abs))

    if verbose and lay.dropped:
        print(f"[figures] {len(lay.dropped)} word(s) did not fit in "
              f"{out_png.name} and were left out: "
              f"{', '.join(lay.dropped[:6])}{'…' if len(lay.dropped) > 6 else ''}")
    with atomic_write(out_png, mode="wb") as fh:
        # the scratch file has no extension, so PIL can't guess the format.
        img.save(fh, format="PNG", optimize=True)
    return out_png

scale_sizes

scale_sizes(
    weights, min_font, max_font, *, few=8, few_floor=0.5
)

Font size for each weight: linear in its magnitude, between a floor and max_font.

Linear, not by rank, because the picture is meant to show how much stronger the top word is, not just that it is first. The floor rises to few_floor * max_font when there are few words or fewer, so a cloud of five significant features fills its frame instead of huddling in the middle at the minimum size. One word gets max_font; weights that are all equal share one middle size, since there is nothing to rank.

Source code in src\taters\figures\render.py

def scale_sizes(weights: Sequence[float], min_font: int, max_font: int, *,
                few: int = 8, few_floor: float = 0.5) -> List[int]:
    """
    Font size for each weight: linear in its magnitude, between a floor and
    ``max_font``.

    Linear, not by rank, because the picture is meant to show *how much*
    stronger the top word is, not just that it is first. The floor rises to
    ``few_floor * max_font`` when there are ``few`` words or fewer, so a cloud
    of five significant features fills its frame instead of huddling in the
    middle at the minimum size. One word gets ``max_font``; weights that are
    all equal share one middle size, since there is nothing to rank.
    """
    mags = [abs(float(w)) for w in weights]
    n = len(mags)
    if n == 0:
        return []
    if n == 1:
        return [int(max_font)]
    lo = int(min_font)
    if n <= few:
        lo = max(lo, round(few_floor * max_font))
    wmin, wmax = min(mags), max(mags)
    if wmax == wmin:
        return [(lo + int(max_font)) // 2] * n
    span = wmax - wmin
    return [round(lo + (max_font - lo) * (v - wmin) / span) for v in mags]

shade

shade(hue, strength)

A color of the given hue whose lightness falls as strength rises.

Weak words come out pale and strong ones dark, so the eye reads strength twice -- in the size and in the ink -- and a strong word that had to shrink to fit still reads as strong.

Source code in src\taters\figures\render.py

def shade(hue: float, strength: float) -> RGB:
    """
    A color of the given hue whose lightness falls as ``strength`` rises.

    Weak words come out pale and strong ones dark, so the eye reads strength
    twice -- in the size and in the ink -- and a strong word that had to
    shrink to fit still reads as strong.
    """
    s = min(1.0, max(0.0, float(strength)))
    lightness = 0.72 - s * (0.72 - 0.26)
    r, g, b = colorsys.hls_to_rgb((hue % 360.0) / 360.0, lightness, 0.65)
    return (round(r * 255), round(g * 255), round(b * 255))

taters.figures.charts

Small charts for training reports: line charts, scatter plots, heat tables.

A training report has three pictures a paper wants: a loss curve (a line per fold or per series over epochs), a predicted-versus-observed scatter, and a confusion matrix as a colored table. None of them needs a plotting library -- the axes, ticks, legend and points are a few hundred lines of arithmetic and Pillow, the same dependency the word clouds already use, so a headless machine draws them and a machine without Pillow gets one plain sentence instead of a traceback.

The layout arithmetic (tick choice, data-to-pixel mapping) is pure Python and tested without a font engine; Pillow is imported inside the functions that draw, as in :mod:taters.figures.render.

Axes dataclass

Axes(left, top, right, bottom, x_lo, x_hi, y_lo, y_hi)

The plotting rectangle in pixels and the data range it shows: the one mapping every mark in a chart goes through.

heat_cell_color

heat_cell_color(value, lo, hi)

White at lo through to a saturated blue at hi; monotone in the value, so a darker cell is always a larger number.

Source code in src\taters\figures\charts.py

def heat_cell_color(value: float, lo: float, hi: float) -> RGB:
    """White at ``lo`` through to a saturated blue at ``hi``; monotone in
    the value, so a darker cell is always a larger number."""
    if hi <= lo or not math.isfinite(float(value)):
        strength = 0.0
    else:
        strength = min(1.0, max(0.0, (float(value) - lo) / (hi - lo)))
    r = int(round(255 - strength * (255 - 31)))
    g = int(round(255 - strength * (255 - 119)))
    b = int(round(255 - strength * (255 - 180)))
    return (r, g, b)

heat_table

heat_table(
    rows,
    out_png,
    *,
    row_labels,
    col_labels,
    title="",
    row_title="",
    col_title="",
    cell_format="{:g}",
    font_path=None,
    on_progress=None
)

A matrix as a colored table: each cell shaded by its value and printed with it. Built for confusion matrices (rows observed, columns predicted), where the diagonal should be the dark one.

Source code in src\taters\figures\charts.py

def heat_table(rows: Sequence[Sequence[float]], out_png: PathLike, *,
               row_labels: Sequence[str], col_labels: Sequence[str],
               title: str = "", row_title: str = "", col_title: str = "",
               cell_format: str = "{:g}", font_path: Optional[PathLike] = None,
               on_progress=None) -> Path:
    """
    A matrix as a colored table: each cell shaded by its value and printed
    with it. Built for confusion matrices (rows observed, columns
    predicted), where the diagonal should be the dark one.
    """
    announce(on_progress, f"drawing chart: {title or Path(out_png).name}")
    n_rows, n_cols = len(rows), len(col_labels)
    cell_w, cell_h, margin = 110, 56, 18
    Image, ImageDraw, ImageFont = _pillow()
    fonts = _Fonts(ImageFont, Path(font_path) if font_path else DEFAULT_FONT)
    label_font, cell_font, title_font = fonts[14], fonts[16], fonts[20]
    label_w = int(max((label_font.getlength(str(x)) for x in row_labels), default=40)) + 20
    label_w += 26 if row_title else 0
    header_h = 30 + (26 if col_title else 0)
    title_h = 34 if title else 0
    width = margin * 2 + label_w + cell_w * n_cols
    height = margin * 2 + title_h + header_h + cell_h * n_rows
    img = Image.new("RGB", (width, height), (255, 255, 255))
    draw = ImageDraw.Draw(img)
    y = margin
    if title:
        l_, t_, r_, b_ = _ink(title_font, title)
        draw.text((margin - l_, y - t_), title, font=title_font, fill=_INK)
        y += title_h
    x0 = margin + label_w
    if col_title:
        l_, t_, r_, b_ = _ink(label_font, col_title)
        draw.text((x0 + (cell_w * n_cols) // 2 - (r_ - l_) // 2 - l_, y - t_),
                  col_title, font=label_font, fill=_MUTED)
        y += 26
    for j, label in enumerate(col_labels):
        l_, t_, r_, b_ = _ink(label_font, str(label))
        draw.text((x0 + j * cell_w + cell_w // 2 - (r_ - l_) // 2 - l_, y - t_),
                  str(label), font=label_font, fill=_INK)
    y += 30
    values = [float(v) for r in rows for v in r if v is not None and math.isfinite(float(v))]
    lo, hi = (min(values), max(values)) if values else (0.0, 1.0)
    if row_title:
        tile_font = label_font
        l_, t_, r_, b_ = _ink(tile_font, row_title)
        tile = Image.new("RGBA", (r_ - l_ + 4, b_ - t_ + 4), (255, 255, 255, 0))
        ImageDraw.Draw(tile).text((2 - l_, 2 - t_), row_title, font=tile_font, fill=_MUTED)
        tile = tile.rotate(90, expand=True)
        img.paste(tile, (margin, y + (cell_h * n_rows) // 2 - tile.height // 2), tile)
    for i, (label, row) in enumerate(zip(row_labels, rows)):
        cy = y + i * cell_h
        l_, t_, r_, b_ = _ink(label_font, str(label))
        draw.text((x0 - 10 - (r_ - l_) - l_, cy + cell_h // 2 - (b_ - t_) // 2 - t_),
                  str(label), font=label_font, fill=_INK)
        for j in range(n_cols):
            value = row[j] if j < len(row) else None
            cx = x0 + j * cell_w
            fill = heat_cell_color(value, lo, hi) if value is not None else (245, 245, 245)
            draw.rectangle([cx, cy, cx + cell_w, cy + cell_h], fill=fill,
                           outline=(255, 255, 255), width=2)
            text = "" if value is None else cell_format.format(value)
            if text:
                l_, t_, r_, b_ = _ink(cell_font, text)
                dark = sum(fill) < 380
                draw.text((cx + cell_w // 2 - (r_ - l_) // 2 - l_,
                           cy + cell_h // 2 - (b_ - t_) // 2 - t_), text,
                          font=cell_font, fill=(255, 255, 255) if dark else _INK)
    return _save(img, Path(out_png))

line_chart

line_chart(
    series,
    out_png,
    *,
    title="",
    x_label="",
    y_label="",
    width=900,
    height=560,
    include_zero=False,
    integer_x=True,
    font_path=None,
    on_progress=None
)

One line per series through its (x, y) points, in PALETTE order.

Built for loss curves: {"fold 1 train": [(1, 0.9), (2, 0.7)], "fold 1 validation": [...]}. With integer_x the x ticks fall on whole numbers (epochs). include_zero forces the y axis to start at zero.

Source code in src\taters\figures\charts.py

def line_chart(series: Mapping[str, Sequence[Point]], out_png: PathLike, *,
               title: str = "", x_label: str = "", y_label: str = "",
               width: int = 900, height: int = 560, include_zero: bool = False,
               integer_x: bool = True, font_path: Optional[PathLike] = None,
               on_progress=None) -> Path:
    """
    One line per series through its ``(x, y)`` points, in ``PALETTE`` order.

    Built for loss curves: ``{"fold 1 train": [(1, 0.9), (2, 0.7)], "fold 1
    validation": [...]}``. With ``integer_x`` the x ticks fall on whole
    numbers (epochs). ``include_zero`` forces the y axis to start at zero.
    """
    announce(on_progress, f"drawing chart: {title or Path(out_png).name}")
    img, draw, fonts = _open(width, height, font_path)
    xs = [p[0] for pts in series.values() for p in pts]
    ys = [p[1] for pts in series.values() for p in pts]
    x_lo, x_hi = _range(xs)
    y_lo, y_hi = _range(ys, include_zero=include_zero)
    x_ticks = nice_ticks(x_lo, x_hi, 8)
    if integer_x:
        x_ticks = sorted({int(round(t)) for t in x_ticks if t == int(t)} or {int(x_lo), int(x_hi)})
        if len(x_ticks) == 1:
            x_ticks = [x_ticks[0] - 1, x_ticks[0], x_ticks[0] + 1]
    y_ticks = nice_ticks(y_lo, y_hi, 6)
    legend = [(name, PALETTE[i % len(PALETTE)]) for i, name in enumerate(series)]
    axes = _frame(draw, fonts, width=width, height=height, title=title,
                  x_label=x_label, y_label=y_label, x_ticks=x_ticks,
                  y_ticks=y_ticks, legend=legend if len(series) > 1 else ())
    for i, (name, pts) in enumerate(series.items()):
        color = PALETTE[i % len(PALETTE)]
        pixels = [(axes.px(x), axes.py(y)) for x, y in pts
                  if y is not None and math.isfinite(float(y))]
        if len(pixels) > 1:
            draw.line(pixels, fill=color, width=3, joint="curve")
        for px, py in pixels:
            draw.ellipse([px - 4, py - 4, px + 4, py + 4], fill=color)
    return _save(img, Path(out_png))

nice_ticks

nice_ticks(lo, hi, n=6)

Round tick positions covering [lo, hi]: about n of them, at a step of 1, 2, 2.5 or 5 times a power of ten, the first at or below lo and the last at or above hi. A flat range gets a unit around its value, so a constant series still has an axis.

Source code in src\taters\figures\charts.py

def nice_ticks(lo: float, hi: float, n: int = 6) -> List[float]:
    """
    Round tick positions covering ``[lo, hi]``: about ``n`` of them, at a
    step of 1, 2, 2.5 or 5 times a power of ten, the first at or below
    ``lo`` and the last at or above ``hi``. A flat range gets a unit
    around its value, so a constant series still has an axis.
    """
    lo, hi = float(lo), float(hi)
    if not math.isfinite(lo) or not math.isfinite(hi):
        lo, hi = 0.0, 1.0
    if hi < lo:
        lo, hi = hi, lo
    if hi == lo:
        pad = abs(lo) * 0.1 or 0.5
        lo, hi = lo - pad, hi + pad
    raw = (hi - lo) / max(1, int(n) - 1)
    power = 10 ** math.floor(math.log10(raw))
    step = next((m * power for m in (1, 2, 2.5, 5, 10) if m * power >= raw), 10 * power)
    first = math.floor(lo / step) * step
    ticks = []
    t = first
    while t <= hi + step * 1e-9:
        ticks.append(round(t, 10))
        t += step
    if ticks[-1] < hi:
        ticks.append(round(t, 10))
    return ticks

pillow_missing_reason

pillow_missing_reason()

None when Pillow can be imported, else the sentence to show instead of a figure.

Source code in src\taters\figures\render.py

def pillow_missing_reason() -> Optional[str]:
    """``None`` when Pillow can be imported, else the sentence to show instead
    of a figure."""
    try:
        _pillow()
    except ImportError as e:
        return str(e)
    return None

scatter_chart

scatter_chart(
    points,
    out_png,
    *,
    title="",
    x_label="",
    y_label="",
    identity=True,
    note="",
    width=720,
    height=720,
    font_path=None,
    on_progress=None
)

Points on square axes; with identity a dashed y = x line, so a predicted-versus-observed plot shows at a glance how far from perfect the predictions fall. note (r, R²) is written inside the axes.

Source code in src\taters\figures\charts.py

def scatter_chart(points: Sequence[Point], out_png: PathLike, *,
                  title: str = "", x_label: str = "", y_label: str = "",
                  identity: bool = True, note: str = "", width: int = 720,
                  height: int = 720, font_path: Optional[PathLike] = None,
                  on_progress=None) -> Path:
    """
    Points on square axes; with ``identity`` a dashed y = x line, so a
    predicted-versus-observed plot shows at a glance how far from perfect
    the predictions fall. ``note`` (r, R²) is written inside the axes.
    """
    announce(on_progress, f"drawing chart: {title or Path(out_png).name}")
    img, draw, fonts = _open(width, height, font_path)
    finite = [(float(x), float(y)) for x, y in points
              if x is not None and y is not None
              and math.isfinite(float(x)) and math.isfinite(float(y))]
    lo, hi = _range([v for p in finite for v in p]) if finite else (0.0, 1.0)
    ticks = nice_ticks(lo, hi, 6)
    axes = _frame(draw, fonts, width=width, height=height, title=title,
                  x_label=x_label, y_label=y_label, x_ticks=ticks, y_ticks=ticks)
    if identity:
        a, b = (axes.px(ticks[0]), axes.py(ticks[0])), (axes.px(ticks[-1]), axes.py(ticks[-1]))
        steps = 40
        for i in range(0, steps, 2):
            x0 = a[0] + (b[0] - a[0]) * i / steps
            y0 = a[1] + (b[1] - a[1]) * i / steps
            x1 = a[0] + (b[0] - a[0]) * (i + 1) / steps
            y1 = a[1] + (b[1] - a[1]) * (i + 1) / steps
            draw.line([(x0, y0), (x1, y1)], fill=(170, 170, 170), width=2)
    color = PALETTE[0] + (150,)
    overlay = img.copy().convert("RGBA")
    from PIL import ImageDraw as _ImageDraw

    layer = _ImageDraw.Draw(overlay, "RGBA")
    for x, y in finite:
        px, py = axes.px(x), axes.py(y)
        layer.ellipse([px - 4, py - 4, px + 4, py + 4], fill=color)
    img.paste(overlay.convert("RGB"))
    if note:
        font = fonts[15]
        draw = _ImageDraw.Draw(img)
        left, top, right, bottom = _ink(font, note)
        draw.text((axes.left + 10 - left, axes.top + 8 - top), note, font=font, fill=_INK)
    return _save(img, Path(out_png))