diff --git a/aaanalysis/feature_engineering/_backend/num_feat/feature_matrix.py b/aaanalysis/feature_engineering/_backend/num_feat/feature_matrix.py new file mode 100644 index 00000000..d36b9b1c --- /dev/null +++ b/aaanalysis/feature_engineering/_backend/num_feat/feature_matrix.py @@ -0,0 +1,138 @@ +""" +This is a script for the backend of ``NumericalFeature.feature_matrix``. + +``get_feature_matrix_num_`` reconstructs the ``(n_samples, n_features)`` value +matrix ``X`` from pre-sliced per-part numerical tensors (``dict_num_parts``), +the numerical-mode analog of :func:`utils_feature.get_feature_matrix_`. + +The value of feature ``PART-SPLIT-SCALE`` for sample ``i`` is: + + ``round(nanmean( arr[part][i, split_positions_i, scale_idx] ), 5)`` + +where ``split_positions_i`` are the **0-based residue indices** obtained by +applying the ``SPLIT`` to ``arange(L_part_i)`` (``L_part_i`` = the sample's real +residue count for that part), and ``scale_idx`` is the D-axis index of ``SCALE`` +in ``df_scales.columns``. This is the SAME reconstruction ``CPP.run_num`` uses in +``recompute_feature_matrix`` (position buffer built from ``arange`` + per-split +``np.nanmean`` + ``np.round(_, 5)``), so the columns of ``X`` are byte-identical +to the values ``run_num`` computed for the same feature ids — verified against +``recompute_feature_matrix`` for uniform and variable-length parts. + +Note that the ``positions`` column of ``run_num``'s ``df_feat`` is a *display* +numbering in TMD-JMD coordinate space (JMD-offset, e.g. ``21..30`` for a TMD), +NOT an index into the ``(L, D)`` tensor; value reconstruction therefore re-applies +the SPLIT to the residue axis rather than reading ``positions``. +""" +import numpy as np +from joblib import Parallel, delayed + +import aaanalysis.utils as ut +from ..cpp._split import SplitVec +from ..cpp.utils_feature import _get_split_info + + +# I Helper Functions +def _feature_column_num(feature=None, dict_num_parts=None, part_lens=None, + scale_to_idx=None, sp_vec=None): + """Compute the ``(n_samples,)`` value column for one ``PART-SPLIT-SCALE`` feature.""" + part_upper, split_str, scale = ut.split_feat_id(feat_id=feature) + part = part_upper.lower() + split_type, split_kwargs = _get_split_info(split=split_str) + f_split = getattr(sp_vec, split_type.lower()) + arr_2d = dict_num_parts[part][:, :, scale_to_idx[scale]] # (n, L_max) float64 + lens = part_lens[part] + n_samples = arr_2d.shape[0] + + # Fast path: every sample shares the same real length for this part, so the + # split positions are identical across samples -> one vectorized gather + + # ``nanmean(axis=1)`` (element-wise identical to the per-sample nanmean). + if n_samples > 0 and bool(np.all(lens == lens[0])): + positions = np.asarray( + f_split(seq=np.arange(int(lens[0]), dtype=np.int64), **split_kwargs), + dtype=np.int64, + ) + if positions.size == 0: + return np.round(np.full(n_samples, np.nan, dtype=np.float64), 5) + col = np.nanmean(arr_2d[:, positions], axis=1) + return np.round(col, 5) + + # General path: positions depend on each sample's real length. + col = np.empty(n_samples, dtype=np.float64) + for i in range(n_samples): + positions = np.asarray( + f_split(seq=np.arange(int(lens[i]), dtype=np.int64), **split_kwargs), + dtype=np.int64, + ) + col[i] = np.nanmean(arr_2d[i, positions]) if positions.size else np.nan + return np.round(col, 5) + + +def _compute_features_num(features=None, dict_num_parts=None, part_lens=None, + scale_to_idx=None): + """Fill an ``(n_samples, n_features)`` matrix for a chunk of features.""" + sp_vec = SplitVec(type_str=False) # one shared instance, hoisted out of the loop + n_samples = next(iter(dict_num_parts.values())).shape[0] + X = np.empty((n_samples, len(features)), dtype=np.float64) + for j, feature in enumerate(features): + X[:, j] = _feature_column_num( + feature=feature, dict_num_parts=dict_num_parts, part_lens=part_lens, + scale_to_idx=scale_to_idx, sp_vec=sp_vec, + ) + return X + + +# II Main Functions +def get_feature_matrix_num_(features=None, dict_num_parts=None, part_lens=None, + df_scales=None, n_jobs=None): + """Build the ``(n_samples, n_features)`` numerical-mode feature matrix ``X``. + + Numerical-mode analog of ``utils_feature.get_feature_matrix_``: instead of an + AA -> scale lookup, per-residue values come from the pre-sliced + ``dict_num_parts`` tensors. Byte-identical to the values ``CPP.run_num`` + computes for the same feature ids (verified against ``recompute_feature_matrix``). + + Parameters + ---------- + features : list of str + Feature ids (``PART-SPLIT-SCALE``); output column order matches. + dict_num_parts : dict[str, np.ndarray] + Per-part ``(n_samples, L_part_max, D)`` NaN-padded numerical tensors, as + produced by ``NumericalFeature.get_parts``. + part_lens : dict[str, np.ndarray] + Per-part ``(n_samples,)`` int64 real residue counts, from + ``_derive_dict_part_lens(df_parts)`` — the SAME length source ``CPP.run_num`` + uses (non-gap character count of the string ``df_parts``). Applying each + split to ``arange(part_lens[part][i])`` therefore lands on exactly the + residues ``run_num`` selected. + df_scales : pd.DataFrame + Names the D dimensions; ``df_scales.columns`` order defines the SCALE -> + D-index mapping. + n_jobs : int, None, or -1 + CPU cores for the feature loop (resolved via ``ut.resolve_n_jobs``). + + Returns + ------- + X : np.ndarray, shape (n_samples, n_features) float64 + Per-sample per-feature values, rounded to 5 decimals. + """ + features = [features] if isinstance(features, str) else list(features) + scale_to_idx = {s: i for i, s in enumerate(list(df_scales.columns))} + + n_jobs = ut.resolve_n_jobs(n_jobs=n_jobs, n_work=len(features)) + if n_jobs == 1 or len(features) <= 1: + return _compute_features_num( + features=features, dict_num_parts=dict_num_parts, + part_lens=part_lens, scale_to_idx=scale_to_idx, + ) + + feature_chunks = np.array_split(features, n_jobs) + + def _mp(chunk): + return _compute_features_num( + features=list(chunk), dict_num_parts=dict_num_parts, + part_lens=part_lens, scale_to_idx=scale_to_idx, + ) + + with Parallel(n_jobs=n_jobs) as parallel: + results = parallel(delayed(_mp)(chunk) for chunk in feature_chunks) + return np.concatenate(results, axis=1) diff --git a/aaanalysis/feature_engineering/_numerical_feature.py b/aaanalysis/feature_engineering/_numerical_feature.py index 61c5cdb6..954b2684 100644 --- a/aaanalysis/feature_engineering/_numerical_feature.py +++ b/aaanalysis/feature_engineering/_numerical_feature.py @@ -9,10 +9,13 @@ import aaanalysis.utils as ut -from ._backend.check_feature import check_df_scales, expand_pos_anchors_ +from ._backend.check_feature import (check_df_scales, expand_pos_anchors_, + check_match_df_scales_features) from ._backend.feature_filter import filter_correlation_ from ._backend.num_feat.extend_alphabet import extend_alphabet_ +from ._backend.num_feat.feature_matrix import get_feature_matrix_num_ from ._backend.cpp._filters._assign import assign_dict_num_to_parts +from ._cpp import _derive_dict_part_lens # I Helper Functions @@ -66,6 +69,37 @@ def check_match_df_seq_dict_num(df_seq=None, dict_num=None) -> None: raise ValueError("'dict_num[*]' has D=0; should be >= 1.") +def check_dict_num_parts(dict_num_parts=None) -> Tuple[int, int]: + """Validate ``dict_num_parts`` (the value source for ``feature_matrix``) and return ``(n_samples, D)``. + + Each value must be a 3-D ``(n_samples, L_part_max, D)`` ndarray, with a + consistent ``n_samples`` (axis 0) and ``D`` (axis 2) across all parts — the + contract emitted by :meth:`NumericalFeature.get_parts`. + """ + ut.check_dict(name="dict_num_parts", val=dict_num_parts, accept_none=False) + if len(dict_num_parts) == 0: + raise ValueError("'dict_num_parts' should not be empty; it must map each part name " + "to a (n_samples, L_part_max, D) tensor (see NumericalFeature.get_parts).") + n_seen, d_seen = set(), set() + for part, arr in dict_num_parts.items(): + if not hasattr(arr, "shape") or getattr(arr, "ndim", None) != 3: + raise ValueError( + f"'dict_num_parts[{part!r}]' should be a 3-D ndarray (n_samples, L_part_max, D); " + f"got ndim={getattr(arr, 'ndim', None)}." + ) + n_seen.add(arr.shape[0]) + d_seen.add(arr.shape[2]) + if len(n_seen) > 1: + raise ValueError(f"'dict_num_parts' has inconsistent n_samples (axis 0) across parts: " + f"{sorted(n_seen)} — all parts must share the same number of samples.") + if len(d_seen) > 1: + raise ValueError(f"'dict_num_parts' has inconsistent D (axis 2) across parts: " + f"{sorted(d_seen)} — all parts must share the same dimensionality.") + if d_seen and 0 in d_seen: + raise ValueError("'dict_num_parts[*]' has D=0; should be >= 1.") + return n_seen.pop(), d_seen.pop() + + # II Main Functions class NumericalFeature: """ @@ -75,7 +109,8 @@ class NumericalFeature: It provides numeric helpers for the :class:`CPP` feature engineering pipeline: extending the amino acid alphabet of a scale DataFrame, slicing per-residue tensors into sequence parts (numerical analog of :meth:`SequenceFeature.get_df_parts`), - and removing redundant features by Pearson correlation. + reconstructing the model matrix ``X`` from :meth:`CPP.run_num`-selected features + (:meth:`feature_matrix`), and removing redundant features by Pearson correlation. .. versionadded:: 0.1.3 """ @@ -277,6 +312,160 @@ def get_parts(df_seq: pd.DataFrame, # no need to expose them as a separate output (keeps the API one-shape). return df_parts, dict_part_vals + @staticmethod + def feature_matrix(features: Union[ut.ArrayLike1D, pd.DataFrame], + dict_num_parts: Dict[str, np.ndarray], + df_parts: pd.DataFrame, + df_scales: Optional[pd.DataFrame] = None, + n_jobs: Union[int, None] = 1, + ) -> np.ndarray: + """ + Create the numerical-mode feature matrix ``X`` for given feature ids and per-part tensors. + + Numerical analog of :meth:`SequenceFeature.feature_matrix`: for each sample and each + feature id, reconstructs the feature value from the pre-sliced per-residue tensors in + ``dict_num_parts`` (PLM embeddings, structure, annotations, ...) instead of an + amino-acid-to-scale lookup, so per-residue context is preserved. The result is the + numerical input ``X`` for a downstream model or for + :meth:`NumericalFeature.filter_correlation`. This is the missing step that turns + :meth:`CPP.run_num`-selected features back into a model matrix. + + .. versionadded:: 1.1.0 + + Parameters + ---------- + features : array-like, shape (n_features,) or pd.DataFrame + Ids of features (``'PART-SPLIT-SCALE'``) for which the matrix of feature values should + be created. Alternatively, a ``df_feat`` DataFrame (e.g. from :meth:`CPP.run_num`), in + which case its ``'feature'`` column is used. + dict_num_parts : dict[str, np.ndarray] + Per-part NaN-padded numerical tensors, as produced by :meth:`NumericalFeature.get_parts` + and consumed by :meth:`CPP.run_num`. Each value has shape ``(n_samples, L_part_max, D)``, + row-aligned across parts; must cover every part referenced in ``features``. + df_parts : pd.DataFrame, shape (n_samples, n_parts) + The string ``df_parts`` returned **alongside** ``dict_num_parts`` by + :meth:`NumericalFeature.get_parts` (row-aligned with it). Supplies each part's real + residue length via the exact same helper :meth:`CPP.run_num` uses internally (non-gap + character count), so every split lands on the residues ``run_num`` selected. Its columns + must cover every part in ``dict_num_parts``. + df_scales : pd.DataFrame, shape (n_letters, n_scales), optional + DataFrame whose columns **name the D dimensions** of ``dict_num_parts`` (the same + ``df_scales`` used to construct the :class:`CPP` for :meth:`CPP.run_num`); its column + order defines the ``SCALE`` -> D-index mapping. Its row (amino acid) values are unused in + numerical mode. Default from :meth:`load_scales` unless specified in + ``options['df_scales']`` — pass your own when the D axis is a custom (e.g. embedding) + space, since the default AA-scale set will not match a custom ``D``. + n_jobs : int, None, or -1, default=1 + Number of CPU cores (>=1) used for multiprocessing. If ``None``, the number is optimized + automatically. If ``-1``, the number is set to all available cores. Overridden by + ``options['n_jobs']`` when set. + + Returns + ------- + X : array-like, shape (n_samples, n_features) + Feature matrix containing feature values for samples. Column ``i`` corresponds to + feature ``i`` in ``features``; values are byte-identical to those :meth:`CPP.run_num` + computed for the same feature ids — both take per-part lengths from the same ``df_parts``, + so they select identical residues for every split. + + Raises + ------ + ValueError + If ``dict_num_parts`` is empty or has inconsistent / non-3D tensors, if ``df_parts`` is + missing a part column or its row count / real lengths do not match ``dict_num_parts``, if + the ``D`` axis does not match ``len(df_scales.columns)``, if a feature id is malformed or + references a part/scale absent from ``dict_num_parts`` / ``df_scales``, or if a feature's + split selects only padded (all-NaN) residues (producing a ``NaN`` value). + + Notes + ----- + * **Mapping ``positions`` back to ``dict_num_parts`` — do not index with it.** The + ``'positions'`` column of a :meth:`CPP.run_num` ``df_feat`` is a **display numbering** in + TMD-JMD coordinate space: positions are offset by the JMD length and use the ``start`` / + ``tmd_len`` / ``jmd_n_len`` / ``jmd_c_len`` shown at run time (e.g. a ``jmd_n_len=10`` TMD + is numbered ``21..30``, decoupled from the actual per-sample TMD length). These numbers do + **not** directly index the ``(L, D)`` per-part array. ``feature_matrix`` therefore does not + read ``positions``; it reconstructs each value the same way :meth:`CPP.run_num` does — by + re-applying the ``SPLIT`` encoded in the feature id to the part's **0-based residue axis** + (``arange(L_part)``, where ``L_part`` is the sample's real residue count taken from + ``df_parts``), selecting the ``SCALE`` column along ``D``, and averaging the selected + residues (``nanmean``, rounded to 5 decimals). Because the split is applied to the residue + axis and not to the display numbering, ``X`` lines up with ``run_num`` by construction. + * **Real per-part lengths come from ``df_parts``**, via the very same length rule + (non-gap character count of ``df_parts``) that :meth:`CPP.run_num` applies + internally — not inferred from the tensor's NaN padding. Passing the ``df_parts`` that + :meth:`NumericalFeature.get_parts` returns alongside ``dict_num_parts`` therefore makes + ``X`` identical to ``run_num`` in every case, including when a *real* residue is all-NaN + across ``D`` (an unresolved structure position or a masked embedding): its length is still + counted from the string, exactly as ``run_num`` counts it. + * A feature whose split selects only padded (all-NaN) residues yields a ``NaN`` value; this + raises a ``ValueError`` rather than returning a silently-``NaN`` column. + * Unlike :meth:`SequenceFeature.feature_matrix` there is no ``accept_gaps`` option: a numeric + tensor carries no gap characters, and NaN padding is ignored via ``nanmean`` automatically. + + See Also + -------- + * :meth:`SequenceFeature.feature_matrix`: sequence-mode (per-AA scale) equivalent. + * :meth:`NumericalFeature.get_parts`: produces the ``dict_num_parts`` consumed here. + * :meth:`CPP.run_num`: selects the feature ids whose values this method reconstructs. + + Examples + -------- + .. include:: examples/nf_feature_matrix.rst + """ + # Load defaults + if df_scales is None: + df_scales = ut.load_default_scales() + # Check input + _n_samples, D = check_dict_num_parts(dict_num_parts=dict_num_parts) + ut.check_df(name="df_parts", df=df_parts, accept_none=False, + cols_required=list(dict_num_parts.keys())) + if len(df_parts) != _n_samples: + raise ValueError( + f"'df_parts' has {len(df_parts)} rows but 'dict_num_parts' has {_n_samples} samples. " + f"Pass the 'df_parts' returned alongside 'dict_num_parts' by NumericalFeature.get_parts." + ) + check_df_scales(df_scales=df_scales) + n_scales = len(df_scales.columns) + if D != n_scales: + raise ValueError( + f"'dict_num_parts' D={D} should equal len(df_scales.columns)={n_scales}. " + f"'df_scales' names the D dimensions in numerical mode — pass the same 'df_scales' " + f"used to construct the CPP for 'run_num'." + ) + features = ut.check_features(features=features, list_parts=list(dict_num_parts.keys()), + list_scales=list(df_scales.columns)) + check_match_df_scales_features(df_scales=df_scales, features=features) + n_jobs = ut.check_n_jobs(n_jobs=n_jobs) + # Real per-part lengths from the SAME source CPP.run_num uses (the string 'df_parts', + # non-gap character count) rather than the tensor's NaN padding, so each split lands on + # exactly the residues run_num selected — identical X even when a real residue is all-NaN. + part_lens = _derive_dict_part_lens(df_parts=df_parts) + for part, arr in dict_num_parts.items(): + if part_lens[part].size and int(part_lens[part].max()) > arr.shape[1]: + raise ValueError( + f"'df_parts[{part!r}]' has a real residue length exceeding the padded tensor " + f"length L_part_max={arr.shape[1]} in 'dict_num_parts'. Pass 'df_parts' and " + f"'dict_num_parts' from the same NumericalFeature.get_parts call." + ) + # Build the feature matrix (byte-identical to CPP.run_num's value reconstruction) + try: + X = get_feature_matrix_num_(features=features, dict_num_parts=dict_num_parts, + part_lens=part_lens, df_scales=df_scales, n_jobs=n_jobs) + except IndexError as error: + raise ValueError( + "A feature's 'SPLIT' references a residue position beyond a part's length in " + "'dict_num_parts'. Ensure 'features' were generated for these parts (e.g. via " + "CPP.run_num on the same 'dict_num_parts')." + ) from error + if not np.isfinite(X).all(): + raise ValueError( + "'feature_matrix' produced NaN feature values: at least one feature's split selects " + "only padded (all-NaN) residues in 'dict_num_parts'. Drop that feature or widen the " + "part/split so it covers real residues." + ) + return X + @staticmethod def extend_alphabet(df_scales: pd.DataFrame, new_letter: str, diff --git a/docs/source/index/release_notes.rst b/docs/source/index/release_notes.rst index ed7d1583..a20d5f69 100644 --- a/docs/source/index/release_notes.rst +++ b/docs/source/index/release_notes.rst @@ -69,6 +69,14 @@ Added missing / non-canonical residues — the sequence's mean profile in scale-space (the scale-based analogue of amino-acid composition). The no-positional-split baseline to compare against ``feature_matrix`` / CPP; optional ``return_df=True`` for a labeled frame. +- :meth:`~aaanalysis.NumericalFeature.feature_matrix`: Turns :meth:`~aaanalysis.CPP.run_num`-selected + features back into a model matrix ``X`` — the numerical analog of + :meth:`~aaanalysis.SequenceFeature.feature_matrix`. Reconstructs each ``PART-SPLIT-SCALE`` value + from the per-residue tensors in ``dict_num_parts``, with per-part lengths taken from ``df_parts`` + (the same length source :meth:`~aaanalysis.CPP.run_num` uses), re-applying the split to the part's + residue axis rather than the JMD-offset ``positions`` display numbering. ``X`` is therefore + byte-identical to the values :meth:`~aaanalysis.CPP.run_num` computed and preserves the per-residue + context that per-AA-averaged sequence features discard. - **SequenceFeature.get_df_parts_from_windows**: Assemble a reference ``df_parts`` from per-part window sets (e.g. ``AAWindowSampler.sample_synthetic`` output). - **SequenceFeature.get_seq_kws**: Return one protein's ``{jmd_n_seq, tmd_seq, jmd_c_seq}`` diff --git a/examples/feature_engineering/nf_feature_matrix.ipynb b/examples/feature_engineering/nf_feature_matrix.ipynb new file mode 100644 index 00000000..1b0a2409 --- /dev/null +++ b/examples/feature_engineering/nf_feature_matrix.ipynb @@ -0,0 +1,704 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e6e83659", + "metadata": {}, + "source": [ + "``NumericalFeature.feature_matrix`` turns :meth:`CPP.run_num`-selected features back into a model matrix ``X``. It reconstructs each feature value from the per-residue tensors in ``dict_num_parts`` (the numerical analog of :meth:`SequenceFeature.feature_matrix`), preserving the per-residue context that per-AA-averaged sequence features discard.\n", + "\n", + "We first build a small per-residue ``dict_num`` (here a synthetic ``[0, 1]``-normalized embedding), slice it into parts with :meth:`NumericalFeature.get_parts`, and run :meth:`CPP.run_num` to obtain a ``df_feat`` of selected features." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c7568921", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T15:30:13.246936Z", + "iopub.status.busy": "2026-07-05T15:30:13.246842Z", + "iopub.status.idle": "2026-07-05T15:30:17.266301Z", + "shell.execute_reply": "2026-07-05T15:30:17.266101Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "DataFrame shape: (15, 13)\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
 featurecategorysubcategoryscale_namescale_descriptionabs_aucabs_mean_difmean_difstd_teststd_refp_val_mann_whitneyp_val_fdr_bhpositions
1JMD_N_TMD_N-Pat...ern(C,2,6)-dim3Embeddingblock0dim3dim30.5000000.558000-0.5580000.0940000.1690000.0209211.00000015,19
2JMD_N_TMD_N-Pat...(C,2,6,10)-dim3Embeddingblock0dim3dim30.5000000.416000-0.4160000.1280000.1100000.0209210.06790911,15,19
3JMD_N_TMD_N-Pat...n(N,10,14)-dim6Embeddingblock1dim6dim60.5000000.4070000.4070000.1590000.1080000.0209210.06681310,14
4TMD-Pattern(C,9,12)-dim2Embeddingblock0dim2dim20.5000000.4050000.4050000.1030000.1080000.0209210.07204219,22
5TMD_C_JMD_C-Pat...rn(N,9,12)-dim2Embeddingblock0dim2dim20.5000000.4050000.4050000.1030000.1080000.0209210.07531729,32
6JMD_N_TMD_N-Pat...ern(N,2,6)-dim3Embeddingblock0dim3dim30.5000000.380000-0.3800000.0840000.1380000.0209210.1090112,6
7TMD_C_JMD_C-Seg...ent(10,12)-dim0Embeddingblock0dim0dim00.5000000.3780000.3780000.1020000.1580000.0209210.10487236
8TMD_C_JMD_C-Pat...ern(N,4,7)-dim3Embeddingblock0dim3dim30.5000000.372000-0.3720000.1420000.0980000.0209210.10103524,27
9TMD-Pattern(C,11,15)-dim7Embeddingblock1dim7dim70.5000000.3660000.3660000.1150000.1080000.0209210.09981716,20
10TMD_C_JMD_C-Pat...rn(N,6,10)-dim7Embeddingblock1dim7dim70.5000000.3660000.3660000.1150000.1080000.0209210.09205426,30
\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import aaanalysis as aa\n", + "import aaanalysis.utils as ut\n", + "aa.options[\"verbose\"] = False\n", + "\n", + "# Tiny example: 8 proteins, a D=8 per-residue embedding, two classes\n", + "n_samples, D = 8, 8\n", + "rng = np.random.default_rng(0)\n", + "seqs = [\"ACDEFGHIKLMNPQRSTVWY\" * 3] * n_samples # length 60\n", + "df_seq = pd.DataFrame({\"entry\": [f\"P{i}\" for i in range(n_samples)], \"sequence\": seqs,\n", + " \"tmd_start\": 11, \"tmd_stop\": 50})\n", + "labels = [1, 1, 1, 1, 0, 0, 0, 0]\n", + "dict_num = {e: rng.random((60, D)) for e in df_seq[\"entry\"]} # already [0, 1]-normalized\n", + "\n", + "# df_scales NAMES the D dimensions (its amino-acid values are unused in numerical mode);\n", + "# df_cat assigns each dimension a category for the CPP redundancy filter.\n", + "dim_names = [f\"dim{i}\" for i in range(D)]\n", + "df_scales = pd.DataFrame(np.zeros((20, D)), index=list(\"ACDEFGHIKLMNPQRSTVWY\"), columns=dim_names)\n", + "df_cat = pd.DataFrame({ut.COL_SCALE_ID: dim_names, ut.COL_CAT: [\"Embedding\"] * D,\n", + " ut.COL_SUBCAT: [f\"block{i // 4}\" for i in range(D)],\n", + " ut.COL_SCALE_NAME: dim_names, ut.COL_SCALE_DES: dim_names})\n", + "\n", + "nf = aa.NumericalFeature()\n", + "df_parts, dict_num_parts = nf.get_parts(df_seq=df_seq, dict_num=dict_num)\n", + "cpp = aa.CPP(df_parts=df_parts, df_scales=df_scales, df_cat=df_cat, verbose=False)\n", + "df_feat = cpp.run_num(dict_num_parts=dict_num_parts, labels=labels, n_filter=15, n_jobs=1)\n", + "aa.display_df(df_feat, n_rows=10, show_shape=True)" + ] + }, + { + "cell_type": "markdown", + "id": "7effe297", + "metadata": {}, + "source": [ + "The required inputs are ``features`` (the ``'feature'`` ids), ``dict_num_parts`` (the per-part tensors from :meth:`get_parts`), and ``df_parts`` (the string parts returned alongside them, which supply each part's real residue length exactly as :meth:`CPP.run_num` derives it); ``df_scales`` names the ``D`` dimensions so the ``SCALE`` token of each feature id maps to the right column of the tensor. ``X`` has shape ``(n_samples, n_features)`` and its values are byte-identical to those :meth:`CPP.run_num` computed for the same feature ids:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "0c4fb9f6", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T15:30:17.267416Z", + "iopub.status.busy": "2026-07-05T15:30:17.267351Z", + "iopub.status.idle": "2026-07-05T15:30:17.272983Z", + "shell.execute_reply": "2026-07-05T15:30:17.272771Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "n samples: 8\n", + "n features: 15\n", + "Shape of X: (8, 15)\n", + "DataFrame shape: (8, 15)\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
 JMD_N_TMD_N-Pattern(C,2,6)-dim3JMD_N_TMD_N-Pattern(C,2,6,10)-dim3JMD_N_TMD_N-Pattern(N,10,14)-dim6TMD-Pattern(C,9,12)-dim2TMD_C_JMD_C-Pattern(N,9,12)-dim2JMD_N_TMD_N-Pattern(N,2,6)-dim3TMD_C_JMD_C-Segment(10,12)-dim0TMD_C_JMD_C-Pattern(N,4,7)-dim3TMD-Pattern(C,11,15)-dim7TMD_C_JMD_C-Pattern(N,6,10)-dim7TMD_C_JMD_C-Segment(4,11)-dim2TMD_C_JMD_C-Segment(5,14)-dim2TMD_C_JMD_C-Segment(5,15)-dim2TMD_C_JMD_C-Pattern(C,8,11)-dim2TMD_C_JMD_C-Pattern(N,5,9)-dim2
entry               
P00.3674700.5704700.7816600.7379200.7379200.1703200.6831400.0931800.8153000.8153000.5602700.5602700.5602700.7271500.676630
P10.1137300.2261700.8004600.5610600.5610600.1584700.6176000.3596300.9539200.9539200.8286700.8286700.8286700.6203800.538700
P20.1757000.3687700.4607700.5553500.5553500.3305400.8846000.4854400.7594800.7594800.5806900.5806900.5806900.7586200.760320
P30.2158800.2997100.4878400.4510900.4510900.1047200.6624000.2866200.6321700.6321700.6148000.6148000.6148000.4031400.532370
P40.9588800.9612800.2972400.1072600.1072600.5610600.2477400.6785800.5269400.5269400.2595300.2595300.2595300.2653000.090950
P50.5243100.6644600.2108400.1059900.1059900.3563400.5038900.7618800.5281600.5281600.3202000.3202000.3202000.2350300.404810
P60.7247000.7444800.3385700.1129300.1129300.7293400.1191300.7540700.3680700.3680700.2838200.2838200.2838200.2188100.149100
P70.8961300.7570800.0554800.3583300.3583300.6391500.4637800.5180700.2745600.2745600.2814400.2814400.2814400.3682200.455140
\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "features = df_feat[\"feature\"].to_list()\n", + "X = nf.feature_matrix(features=features, dict_num_parts=dict_num_parts, df_parts=df_parts, df_scales=df_scales)\n", + "print(f\"n samples: {X.shape[0]}\")\n", + "print(f\"n features: {X.shape[1]}\")\n", + "print(f\"Shape of X: {X.shape}\")\n", + "aa.display_df(pd.DataFrame(X, index=df_parts.index, columns=features), n_rows=10, show_shape=True)" + ] + }, + { + "cell_type": "markdown", + "id": "b5a30c0d", + "metadata": {}, + "source": [ + "Instead of a list of feature ids, the ``df_feat`` DataFrame can be passed directly as ``features`` — its ``'feature'`` column is used automatically, so ``features=df_feat`` is equivalent to ``features=list(df_feat['feature'])``:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "e3e1d6dd", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T15:30:17.274079Z", + "iopub.status.busy": "2026-07-05T15:30:17.273999Z", + "iopub.status.idle": "2026-07-05T15:30:17.276724Z", + "shell.execute_reply": "2026-07-05T15:30:17.276358Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Shape of X: (8, 15)\n", + "Identical to list-of-ids result: True\n" + ] + } + ], + "source": [ + "X_from_df_feat = nf.feature_matrix(features=df_feat, dict_num_parts=dict_num_parts, df_parts=df_parts, df_scales=df_scales)\n", + "print(f\"Shape of X: {X_from_df_feat.shape}\")\n", + "print(f\"Identical to list-of-ids result: {np.array_equal(X, X_from_df_feat)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "116a4a6c", + "metadata": {}, + "source": [ + "``df_scales`` names the ``D`` dimensions of ``dict_num_parts`` and must be the **same** ``df_scales`` used to construct the :class:`CPP` for :meth:`run_num`; its column order defines the ``SCALE`` → dimension mapping. When omitted it defaults to the bundled amino-acid scales, which only fits when ``D`` matches that default — for a custom embedding space, always pass your own ``df_scales`` (a wrong ``D`` raises a clear error):" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "cec8bc92", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T15:30:17.278184Z", + "iopub.status.busy": "2026-07-05T15:30:17.278067Z", + "iopub.status.idle": "2026-07-05T15:30:17.280895Z", + "shell.execute_reply": "2026-07-05T15:30:17.280670Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "features: ['TMD_C_JMD_C-Segment(10,12)-dim0']\n", + "Shape of X: (8, 1)\n" + ] + } + ], + "source": [ + "# Selecting a single embedding dimension by name still reproduces the run_num values\n", + "one_dim_feats = [f for f in features if f.endswith(\"dim0\")] or features[:1]\n", + "X_one = nf.feature_matrix(features=one_dim_feats, dict_num_parts=dict_num_parts, df_parts=df_parts, df_scales=df_scales)\n", + "print(f\"features: {one_dim_feats}\")\n", + "print(f\"Shape of X: {X_one.shape}\")" + ] + }, + { + "cell_type": "markdown", + "id": "d4c0cdcb", + "metadata": {}, + "source": [ + "Multiprocessing is controlled by ``n_jobs`` (set to the maximum when ``n_jobs=None``, all cores when ``n_jobs=-1``). It is only worthwhile for large feature sets (>~1000 features per core) because of process-management overhead; the result is identical regardless of ``n_jobs``:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "0401dad8", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T15:30:17.281947Z", + "iopub.status.busy": "2026-07-05T15:30:17.281885Z", + "iopub.status.idle": "2026-07-05T15:30:17.284936Z", + "shell.execute_reply": "2026-07-05T15:30:17.284602Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Identical to default (n_jobs=1) result: True\n" + ] + } + ], + "source": [ + "X_serial = nf.feature_matrix(features=features, dict_num_parts=dict_num_parts,\n", + " df_parts=df_parts, df_scales=df_scales, n_jobs=1)\n", + "print(f\"Identical to default (n_jobs=1) result: {np.array_equal(X, X_serial)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "60d5b322", + "metadata": {}, + "source": [ + "**Mapping ``positions`` back to the tensor.** The ``'positions'`` column of ``df_feat`` is a *display* numbering in TMD-JMD coordinate space (JMD-offset, e.g. ``21..30`` for a TMD), **not** an index into the ``(L, D)`` per-part array. ``feature_matrix`` therefore never reads ``positions``; it reconstructs each value by re-applying the ``SPLIT`` in the feature id to the part's 0-based residue axis — exactly as :meth:`CPP.run_num` does — so ``X`` lines up with ``run_num`` by construction:" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "c64a9f19", + "metadata": { + "execution": { + "iopub.execute_input": "2026-07-05T15:30:17.286467Z", + "iopub.status.busy": "2026-07-05T15:30:17.286352Z", + "iopub.status.idle": "2026-07-05T15:30:17.288913Z", + "shell.execute_reply": "2026-07-05T15:30:17.288677Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "feature: JMD_N_TMD_N-Pattern(C,2,6)-dim3\n", + "df_feat 'positions' (display numbering, NOT tensor indices): 15,19\n" + ] + } + ], + "source": [ + "sample_feat = features[0]\n", + "sample_pos = df_feat.loc[df_feat[\"feature\"] == sample_feat, \"positions\"].iloc[0]\n", + "print(f\"feature: {sample_feat}\")\n", + "print(f\"df_feat 'positions' (display numbering, NOT tensor indices): {sample_pos}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.14.0" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tests/unit/numerical_feature_tests/test_nf_feature_matrix.py b/tests/unit/numerical_feature_tests/test_nf_feature_matrix.py new file mode 100644 index 00000000..31255d3e --- /dev/null +++ b/tests/unit/numerical_feature_tests/test_nf_feature_matrix.py @@ -0,0 +1,386 @@ +"""Tests for NumericalFeature.feature_matrix(). + +``nf.feature_matrix(features, dict_num_parts, df_scales=...)`` reconstructs the +``(n_samples, n_features)`` model matrix ``X`` from the per-part numerical tensors +``dict_num_parts`` (numerical analog of ``SequenceFeature.feature_matrix``). The +values are byte-identical to those ``CPP.run_num`` computes for the same feature +ids — the golden-value tests below pin that against ``recompute_feature_matrix`` +(the exact backend ``run_num`` uses) and against hand-computed means. +""" +import warnings + +import numpy as np +import pandas as pd +import pytest + +import aaanalysis as aa +import aaanalysis.utils as ut +from aaanalysis.feature_engineering._backend.cpp._filters._recompute import recompute_feature_matrix +from aaanalysis.feature_engineering._cpp import _derive_dict_part_lens + +aa.options["verbose"] = False +warnings.filterwarnings("ignore") + + +# I Helper functions +def _name_dims(D): + """Synthetic df_scales / df_cat naming D embedding dimensions.""" + dim_names = [f"dim{i}" for i in range(D)] + df_scales = pd.DataFrame(np.zeros((20, D)), index=list("ACDEFGHIKLMNPQRSTVWY"), columns=dim_names) + df_cat = pd.DataFrame({ut.COL_SCALE_ID: dim_names, ut.COL_CAT: ["Emb"] * D, + ut.COL_SUBCAT: [f"b{i // 2}" for i in range(D)], + ut.COL_SCALE_NAME: dim_names, ut.COL_SCALE_DES: dim_names}) + return df_scales, df_cat + + +def _run_num_fixture(n=8, D=6, seed=0, stops=None): + """Deterministic (df_parts, dict_num_parts, df_scales, df_cat, features) via get_parts + run_num.""" + rng = np.random.default_rng(seed) + seqs = ["ACDEFGHIKLMNPQRSTVWY" * 3] * n # length 60 + if stops is None: + stops = [50] * n + df_seq = pd.DataFrame({"entry": [f"P{i}" for i in range(n)], "sequence": seqs, + "tmd_start": [11] * n, "tmd_stop": stops}) + dict_num = {e: rng.random((60, D)) for e in df_seq["entry"]} + labels = [1] * (n // 2) + [0] * (n - n // 2) + df_scales, df_cat = _name_dims(D) + nf = aa.NumericalFeature() + df_parts, dict_num_parts = nf.get_parts(df_seq=df_seq, dict_num=dict_num) + cpp = aa.CPP(df_parts=df_parts, df_scales=df_scales, df_cat=df_cat, verbose=False) + df_feat = cpp.run_num(dict_num_parts=dict_num_parts, labels=labels, n_filter=15, n_jobs=1) + return dict(df_parts=df_parts, dict_num_parts=dict_num_parts, df_scales=df_scales, + df_cat=df_cat, df_feat=df_feat, features=df_feat["feature"].to_list(), + cpp=cpp, labels=labels) + + +# II Per-parameter tests +class TestFeatureMatrix: + """Per-parameter positive and negative coverage.""" + + # features + def test_valid_features_list(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert X.shape == (len(fx["df_parts"]), len(fx["features"])) + + def test_valid_features_df_feat(self): + """A df_feat DataFrame is accepted: its 'feature' column is used.""" + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X_ids = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + X_df = nf.feature_matrix(features=fx["df_feat"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert np.array_equal(X_ids, X_df) + + def test_valid_features_single_str(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"][0], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert X.shape == (len(fx["df_parts"]), 1) + + def test_invalid_features_empty(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=[], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + + def test_invalid_features_malformed(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=["TMD-Segment(1,2)"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + + def test_invalid_features_unknown_part(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=["NOPART-Segment(1,2)-dim0"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + + def test_invalid_features_unknown_scale(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=["TMD-Segment(1,2)-nope"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + + # dict_num_parts + def test_valid_dict_num_parts(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert np.isfinite(X).all() + + def test_invalid_dict_num_parts_none(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=None, df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + + def test_invalid_dict_num_parts_empty(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts={}, df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + + def test_invalid_dict_num_parts_2d(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + bad = {p: a[:, :, 0] for p, a in fx["dict_num_parts"].items()} # drop D axis -> 2D + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=bad, df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + + def test_invalid_dict_num_parts_inconsistent_n(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + bad = dict(fx["dict_num_parts"]) + first = list(bad)[0] + bad[first] = bad[first][:-1] # different n_samples + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=bad, df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + + def test_invalid_dict_num_parts_inconsistent_D(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + bad = dict(fx["dict_num_parts"]) + first = list(bad)[0] + bad[first] = bad[first][:, :, :-1] # different D + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=bad, df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + + # df_parts + def test_valid_df_parts_lengths_drive_splits(self): + """df_parts supplies the real lengths; passing get_parts' df_parts yields finite X.""" + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], + df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + assert np.isfinite(X).all() + + def test_invalid_df_parts_row_mismatch(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], + df_parts=fx["df_parts"].iloc[:-1], df_scales=fx["df_scales"]) + + def test_invalid_df_parts_missing_part_column(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + part0 = list(fx["dict_num_parts"])[0] + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], + df_parts=fx["df_parts"].drop(columns=[part0]), df_scales=fx["df_scales"]) + + def test_invalid_df_parts_length_exceeds_tensor(self): + """A df_parts real length beyond the padded tensor length is a get_parts mismatch.""" + fx = _run_num_fixture() + nf = aa.NumericalFeature() + part0 = list(fx["dict_num_parts"])[0] + l_max = fx["dict_num_parts"][part0].shape[1] + bad = fx["df_parts"].copy() + bad[part0] = ["A" * (l_max + 1)] * len(bad) # real length exceeds tensor L_part_max + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], + df_parts=bad, df_scales=fx["df_scales"]) + + # df_scales + def test_valid_df_scales_custom(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert X.shape[1] == len(fx["features"]) + + def test_valid_df_scales_default_when_D_matches(self): + """Omitting df_scales falls back to the bundled scales; works only if D matches.""" + df_default = ut.load_default_scales() + D = len(df_default.columns) + fx = _run_num_fixture(D=D) + # Feature scales come from our synthetic dim-names, not the default columns, + # so features must reference default scale ids for the default path. Build a + # trivial feature on a default scale id instead. + scale0 = list(df_default.columns)[0] + feat = f"TMD-Segment(1,1)-{scale0}" + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=[feat], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"]) + assert X.shape == (len(fx["df_parts"]), 1) + + def test_invalid_df_scales_D_mismatch(self): + fx = _run_num_fixture(D=6) + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"].iloc[:, :4]) # D=4 != 6 + + def test_invalid_df_scales_nan(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + bad = fx["df_scales"].copy() + bad.iloc[0, 0] = np.nan + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], df_scales=bad) + + # n_jobs + def test_valid_n_jobs(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X_ref = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"], n_jobs=1) + for n_jobs in [1, -1, None]: + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"], n_jobs=n_jobs) + assert np.array_equal(X, X_ref), f"n_jobs={n_jobs} changed the result" + + def test_invalid_n_jobs(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"], n_jobs=0) + + +# III Complex / contract tests +class TestFeatureMatrixComplex: + """Shape, dtype, column ordering, and cross-form consistency.""" + + def test_output_shape_and_dtype(self): + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert isinstance(X, np.ndarray) + assert X.dtype == np.float64 + assert X.shape == (len(fx["df_parts"]), len(fx["features"])) + + def test_column_order_matches_features(self): + """Reordering the feature list permutes X's columns the same way.""" + fx = _run_num_fixture() + nf = aa.NumericalFeature() + feats = fx["features"] + X = nf.feature_matrix(features=feats, dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + perm = list(reversed(range(len(feats)))) + feats_perm = [feats[i] for i in perm] + X_perm = nf.feature_matrix(features=feats_perm, dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + assert np.array_equal(X_perm, X[:, perm]) + + def test_variable_tmd_length(self): + """Ragged (variable-TMD) parts are handled and stay run_num-consistent.""" + fx = _run_num_fixture(stops=[50, 50, 50, 50, 45, 45, 45, 45]) + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + X_ref = recompute_feature_matrix( + dict_part_vals=fx["dict_num_parts"], + dict_part_lens=_derive_dict_part_lens(df_parts=fx["df_parts"]), + list_scales=list(fx["df_scales"].columns), features=fx["features"], + split_kws=fx["cpp"].split_kws, + ) + assert np.array_equal(X, X_ref, equal_nan=True) + + +# IV Golden-value tests +class TestFeatureMatrixGoldenValues: + """Pin exact values against hand computation and run_num's own backend.""" + + def test_consistency_with_run_num_recompute(self): + """feature_matrix == the exact matrix run_num built for the same feature ids.""" + fx = _run_num_fixture() + nf = aa.NumericalFeature() + X = nf.feature_matrix(features=fx["features"], dict_num_parts=fx["dict_num_parts"], df_parts=fx["df_parts"], + df_scales=fx["df_scales"]) + X_ref = recompute_feature_matrix( + dict_part_vals=fx["dict_num_parts"], + dict_part_lens=_derive_dict_part_lens(df_parts=fx["df_parts"]), + list_scales=list(fx["df_scales"].columns), features=fx["features"], + split_kws=fx["cpp"].split_kws, + ) + assert np.array_equal(X, X_ref, equal_nan=True) + + def test_all_nan_real_residue_matches_run_num(self): + """A genuine (non-padding) residue that is all-NaN across D still matches run_num. + + The length is taken from the string ``df_parts`` (as run_num does), not inferred + from the tensor's NaN pattern, so the split boundaries do not shift. Uses whole / + half-part segments that span the blanked residue plus others, so ``nanmean`` stays + finite (a split selecting only the blanked residue would legitimately raise). + """ + fx = _run_num_fixture() + nf = aa.NumericalFeature() + dict_num_parts = dict(fx["dict_num_parts"]) + tmd = dict_num_parts["tmd"].copy() + tmd[:, 5, :] = np.nan # blank a real interior TMD residue for every sample + dict_num_parts["tmd"] = tmd + feats = ["TMD-Segment(1,1)-dim0", "TMD-Segment(1,2)-dim1"] + X = nf.feature_matrix(features=feats, dict_num_parts=dict_num_parts, + df_parts=fx["df_parts"], df_scales=fx["df_scales"]) + X_ref = recompute_feature_matrix( + dict_part_vals=dict_num_parts, + dict_part_lens=_derive_dict_part_lens(df_parts=fx["df_parts"]), + list_scales=list(fx["df_scales"].columns), features=feats, + split_kws=fx["cpp"].split_kws, + ) + assert np.array_equal(X, X_ref, equal_nan=True) + + def test_golden_whole_part_segment(self): + """Segment(1,1) = whole-part mean of the chosen dimension (hand-computed).""" + n, L, D = 4, 10, 3 + rng = np.random.default_rng(7) + arr = rng.random((n, L, D)) + dict_num_parts = {"tmd": arr} + df_parts = pd.DataFrame({"tmd": ["A" * L] * n}) # real length L (dense, no padding) + df_scales, _ = _name_dims(D) + nf = aa.NumericalFeature() + feat = "TMD-Segment(1,1)-dim2" + X = nf.feature_matrix(features=[feat], dict_num_parts=dict_num_parts, df_parts=df_parts, df_scales=df_scales) + expected = np.round(arr[:, :, 2].mean(axis=1), 5) + assert np.allclose(X[:, 0], expected) + + def test_golden_first_half_segment(self): + """Segment(1,2) = mean over the first half of the part (int floor split).""" + n, L, D = 4, 10, 3 + rng = np.random.default_rng(8) + arr = rng.random((n, L, D)) + dict_num_parts = {"tmd": arr} + df_parts = pd.DataFrame({"tmd": ["A" * L] * n}) # real length L (dense, no padding) + df_scales, _ = _name_dims(D) + nf = aa.NumericalFeature() + feat = "TMD-Segment(1,2)-dim0" + X = nf.feature_matrix(features=[feat], dict_num_parts=dict_num_parts, df_parts=df_parts, df_scales=df_scales) + expected = np.round(arr[:, 0:5, 0].mean(axis=1), 5) # first half of L=10 + assert np.allclose(X[:, 0], expected) + + def test_golden_pattern_positions(self): + """Pattern(N,1,3) selects 0-based residues 0 and 2 (list_pos - 1).""" + n, L, D = 4, 10, 3 + rng = np.random.default_rng(9) + arr = rng.random((n, L, D)) + dict_num_parts = {"tmd": arr} + df_parts = pd.DataFrame({"tmd": ["A" * L] * n}) # real length L (dense, no padding) + df_scales, _ = _name_dims(D) + nf = aa.NumericalFeature() + feat = "TMD-Pattern(N,1,3)-dim1" + X = nf.feature_matrix(features=[feat], dict_num_parts=dict_num_parts, df_parts=df_parts, df_scales=df_scales) + expected = np.round(arr[:, [0, 2], 1].mean(axis=1), 5) + assert np.allclose(X[:, 0], expected) + + def test_nan_split_raises(self): + """A split that selects only all-NaN residues yields a NaN value and raises a clear ValueError.""" + n, L, D = 3, 6, 2 + arr = np.full((n, L, D), np.nan) # every selected residue is NaN + dict_num_parts = {"tmd": arr} + df_parts = pd.DataFrame({"tmd": ["A" * L] * n}) # real length L (dense, no padding) + df_scales, _ = _name_dims(D) + nf = aa.NumericalFeature() + with pytest.raises(ValueError): + nf.feature_matrix(features=["TMD-Segment(1,1)-dim0"], dict_num_parts=dict_num_parts, + df_parts=df_parts, df_scales=df_scales)