fix(pp): index top peaks by region in normalize_peaks

[nkempynck]

A bit of a concerning indexing bug that luckily does not have huge effects.

Sorting the thresholded column gives positions within that column, which are
offset from region indices by the number of regions dropped ahead of them. Those
positions were used to index the full region matrix, so the Gini scores deciding
the low-Gini subset were read off unrelated regions whenever a region fell at or
below peak_threshold.

That in effect skipped the broad/specific test: the peaks averaged for the weight
were the right top peaks, but were never checked for being broad. Conversely the
same wrong indices were stored, so the returned region DataFrame listed genuinely
broad regions that were not top peaks.

Carry the kept region indices through the top-k selection instead.

Weights move ~1% on average (2.8% max) on the mouse cortex tutorial data.

Second commit: normalize_peaks now raises when a cell type ends up with no
selected peaks, instead of giving it an inf weight (or nan for all of them, if
it happens to every cell type) and writing that into .X. Reachable by raising
gini_std_threshold past the point where no top peak still counts as broad, which
is dataset dependent — 2.5 on the mouse cortex tutorial data, right next to the
2.0 in the docstring example.

[nkempynck]

thanks to @UCDNJJ for noticing!!

[nkempynck]

still have to update the tutorial notebook

[Copilot]

Pull request overview

Fixes correctness issues in crested.pp.normalize_peaks where peak indices were mis-mapped after thresholding (causing Gini-based broad/specific filtering and returned peak lists to be computed from unrelated regions), and hardens normalization by raising when a cell type ends up with no selected peaks (instead of silently writing inf/nan into .X).

Changes:

• Correctly carries original region indices through top-k selection when peak_threshold filters regions prior to sorting.
• Raises a clear ValueError when no peaks pass top-k + Gini selection for one or more cell types.
• Adds closed-form tests covering the indexing bug and the new error behavior; documents both changes in the changelog.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 1 comment.

File	Description
src/crested/pp/_normalization.py	Fixes index mapping for top-k selection; adds an explicit error when normalization weights would be undefined.
tests/test_pp.py	Adds analytic regression tests for correct top-peak/Gini coupling and for raising on empty selections.
docs/changelog.md	Documents the indexing bugfix, the new error behavior, and the added regression test.

---

💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

[nkempynck]

ready to merge if all agree

[LukasMahieu]

avoid exact paths to our storage. Just remove notebook output cells that include these.

[casblaauw]

I find the code a bit confusing to follow and the comments are very verbose as always with generated stuff, so I'll clean up the comments a bit

[casblaauw]

Had a look over it and adjusted some stuff, here's a summary:

• Separated it out into sections with more comments and renamed variables to be clearer to me
• Replaced saving as 0 (and then trying to check that with readying arbitrary ints as booleans, which feels risky) for failed classes with just saving the class name
• Reduced gini to 1D vector once rather than doing it separately every time
• Drop redundant per-class gini calculating (since gini is an independent property of every region)
• Dropped a redundant duplicated sort/argsort call (sorting filtered_col for both sorted_col and top_indices)
• Dropped creating of lists in the loop that aren't used elsewhere
• Also raise an error if adata.obsm['weights'] already exists, rather than trying to run the normalization on already-normalized data

Here's a version of my reworked inner filtering step with annotations of which indices belong to which, so that you can check my logic:

    for i in range(target_matrix.shape[1]):
        # Apply peak_threshold: minimum peak height filtering for this cell type
        filtered_indices = np.where(target_matrix[:, i] > peak_threshold)[0] # filtered_indices indices into target_matrix
        filtered_col = target_matrix[filtered_indices, i]

        # Apply top_k_percent: Get top k of the values that pass threshold
        top_k_index = int(len(filtered_indices) * top_k_percent)
        sorted_filtered_indices = np.argsort(filtered_col)[::-1] # sorted_filtered_indices indices into filtered_[col|indices]
        top_indices = filtered_indices[sorted_filtered_indices[:top_k_index]] # top_indices is a subset of filtered_indices, so indices into target_matrix
        top_col = target_matrix[:, i][top_indices] # aka filtered_col[sorted_filtered_indices[:top_k_index]]

        # Apply gini_threshold: Retrieve gini scores for the regions passing previous filters and keep those below threshold
        gini_scores = overall_gini_scores[top_indices]
        low_gini_mask = gini_scores < gini_threshold # indices into top_[col|indices]

        # Save final scores
        if low_gini_mask.sum() > 0:
            top_k_mean = np.mean(top_col[low_gini_mask])
            all_low_gini_indices.update(top_indices[low_gini_mask])
            top_k_percent_means.append(top_k_mean)
        else:
            # No regions left after all filtering, saving cell type name to raise later
            failed_classes.append(adata.obs_names[i])

I checked my changes and they return identical values to Niklas' rework:

# Main branch, no changes
[[1.275371 ]
 [1.0847164]
 [2.2603295]
 [3.0955002]
 [2.2730987]
 [2.19935  ]
 [2.2757301]
 [2.3690026]
...]
# Niklas' branch, no Cas changes
[[1.2876699]
 [1.116849 ]
 [2.307081 ]
 [3.0326943]
 [2.3098855]
 [2.2099555]
 [2.3125596]
 [2.3567338]
...]
# Niklas branch with Cas changes
[[1.2876699]
 [1.116849 ]
 [2.307081 ]
 [3.0326943]
 [2.3098855]
 [2.2099555]
 [2.3125596]
 [2.3567338]
...]

[casblaauw]

Changes made manually

[nkempynck]

thanks cas!