One-pass qmap-partial specificity metric, replacing fast and qnorm spec

docs/gcprunner.md

@@ -852,10 +852,6 @@ mounted at `/workspace/cache`. Use this when the weights are actually on disk.
 - `--test` / `--valid` evaluate only the held-out test / validation fold per
   model (from `folds.json`); the default evaluates every fold.
 - `--save` / `--aggregate_genes` need ~60 GB RAM — use `-q l4-large`.
-- `--spec` (specificity) is **auto-pinned to `-q l4-large`** (g2-standard-32,
-  32 vCPU, ~122 GB) regardless of `-q`, since quantile normalization is RAM- and
-  CPU-hungry; it runs with `--ncpus 32`. Your `-q` still governs the lighter
-  coverage-eval jobs.
 - On GCP, `mem` is the VM's hard RAM (the task owns the whole VM), not a soft
   Slurm request — size it via `-q`, not a mem flag. Eval warns if the chosen
   profile is too small.

docs/train.md

@@ -126,9 +126,9 @@ Each epoch, `log.txt` gets a summary line per dataset (`loss`, `r`, `r2`, `spec`
 Tracks are grouped by the targets `group` column (else the `description` prefix). Each group with at least `train.spec_group_min` tracks (default 20, after merging strand pairs) gets:
 
 - `r/<group>`, `r2/<group>`: mean per-track Pearson r and R² over the group's tracks.
-- `spec/<group>`: specificity, the mean per-track Pearson r after dividing every track by its dataset-wide mean and subtracting the group mean at each position. It rewards predicting how a track differs from its group, not the shared signal.
+- `spec/<group>`: specificity, the mean per-track Pearson r of what remains after removing the group's shared signal. Each track's values (strand pairs summed) are quantile-normalized to the group with a fixed per-track lookup table; predictions use the same table as their targets. The group mean at each position is then regressed out of each track, with the track's own slope and intercept. It rewards predicting how a track differs from its group, not the shared signal or the track's depth and signal-to-noise.
 
-`spec` is the unweighted mean of `spec/<group>` over groups. It needs the per-track means that `hound_data` stores in each `examples/*.zarr`; for an older dataset, add them with `python -c "from baskerville import dataset; dataset.write_target_means('data_dir')"`.
+`spec` is the unweighted mean of `spec/<group>` over groups. `hound_eval` reports the same metric per track (`spec` column of `acc.txt`). The tables are built from the exact whole-genome histograms of each track's stored fp16 values, which `hound_data` writes to each `examples/*.zarr` as `target_hist`. For an older dataset, add them with `python -c "from baskerville import dataset; dataset.write_target_hist('data_dir')"`.
 
 `train.stop_stat` chooses the statistic that selects `model_best.pth` and drives `patience`. It is a `{metric: weight}` dict over the keys above plus `loss`, `r_gene`, and `r2_gene`, and the weighted sum is maximized. The presets `"weighted_r_r2"` (default, `{"r": 1, "r2": 0.25}`), `"r"`, `"r2"`, and `"loss"` (`{"loss": -1}`) remain available as strings; these presets score gene-only datasets by −loss. A dict applies only its own keys, e.g. `{"r_gene": 1}` or `{"loss": -1}` for gene-only datasets.
 

docs/train_folds.md

@@ -258,8 +258,6 @@ hound_eval_folds --backend local --rank -o train_folds params.json data_folds
 # Evaluate only test sets (faster)
 hound_eval_folds --backend local --test -o train_folds params.json data_folds
 
-# Include specificity analysis
-hound_eval_folds --backend local --spec -o train_folds params.json data_folds
 ```
 
 ### Understanding Evaluation Output
@@ -269,20 +267,18 @@ After evaluation, each fold directory contains:
 ```
 train_folds/f0c0/
 ├── train/          # Training artifacts
-├── eval/           # Evaluation results
-│   ├── fold0/      # Model evaluated on fold 0
-│   ├── fold1/      # Model evaluated on fold 1
-│   ├── ...
-│   ├── test -> fold0  # Symlink to test fold
-│   └── test.out
-└── spec/           # Specificity analysis (if --spec used)
+└── eval/           # Evaluation results
+    ├── fold0/      # Model evaluated on fold 0
+    ├── fold1/      # Model evaluated on fold 1
+    ├── ...
+    ├── test -> fold0  # Symlink to test fold
+    └── test.out
 ```
 
 Key files:
 
 - `eval/fold{X}/acc.txt`: Accuracy metrics for each target
-- `eval/test/acc.txt`: Test set performance (most important)
-- `spec/acc.txt`: Specificity analysis results
+- `eval/test/acc.txt`: Test set performance (most important), with per-track `spec` (see [train.md](train.md))
 
 ## Analyzing Results
 

pyproject.toml

@@ -80,7 +80,6 @@ hound_eval = "baskerville.scripts.hound_eval:main"
 hound_eval_folds = "baskerville.scripts.hound_eval_folds:main"
 hound_eval_genes = "baskerville.scripts.hound_eval_genes:main"
 hound_eval_genes_folds = "baskerville.scripts.hound_eval_genes_folds:main"
-hound_eval_spec = "baskerville.scripts.hound_eval_spec:main"
 hound_grad = "baskerville.scripts.hound_grad:main"
 hound_grad_folds = "baskerville.scripts.hound_grad_folds:main"
 hound_ism_bed = "baskerville.scripts.hound_ism_bed:main"

src/baskerville/dataset.py

@@ -187,10 +187,10 @@ def __init__(
         # don't persist handles (not fork-safe with DataLoader workers)
         zarr_data = [zarr.open(zf, mode="r") for zf in self.zarr_files]
 
-        # per-track means of the stored targets (write_target_means)
-        self.target_means = None
-        if self.has_coverage and "mean" in zarr_data[0]["target"].attrs:
-            self.target_means = np.array(zarr_data[0]["target"].attrs["mean"])
+        # per-track histograms of the stored targets (write_target_hist)
+        self.target_hist = None
+        if self.has_coverage and "target_hist" in zarr_data[0]:
+            self.target_hist = zarr_data[0]["target_hist"][:]
 
         # count sequences
         zarr_seqs = [zd["sequence"].shape[0] for zd in zarr_data]
@@ -864,47 +864,71 @@ def target_groups(targets_df):
     return np.array(groups)
 
 
-def _target_block_sum(args):
-    """Float64 per-track sum of the stored targets over sequences [start, end)."""
-    zarr_file, start, end = args
-    targets = zarr.open(zarr_file, mode="r")["target"]
-    block_sum = np.zeros(targets.shape[1])
-    for si in range(start, end):
-        block_sum += targets[si].astype("float64").sum(axis=-1)
-    return block_sum
-
+# non-negative finite fp16 bit patterns: 0x0000 (0.0) to 0x7BFF (65504)
+NUM_HIST_BINS = 0x7C00
 
-def write_target_means(data_dir, processes=16, block_seqs=64):
-    """Store exact per-track means of the stored targets in each zarr.
 
-    The mean is over every position of every sequence in all examples/*.zarr,
-    in the stored (transformed) space, and is written to each zarr's
-    ``target.attrs["mean"]``.
+def _target_block_hist(args):
+    """Per-track counts of strand-summed fp16 target bit patterns over
+    sequences [start, end)."""
+    zarr_file, start, end, pair = args
+    targets = zarr.open(zarr_file, mode="r")["target"]
+    # int32: a block's counts per bin are at most (end - start) * length
+    hist = np.zeros((targets.shape[1], NUM_HIST_BINS), dtype=np.int32)
+    for si in range(start, end):
+        y = targets[si].astype(np.float32)
+        if not np.isfinite(y).all() or (y < 0).any():
+            raise ValueError(f"{zarr_file} seq {si}: negative or non-finite targets")
+        # + 0.0 turns -0.0 into 0.0
+        summed = np.minimum(y + y[pair], 65504) + 0.0
+        bits = summed.astype(np.float16).view(np.uint16)
+        for ti, track_bits in enumerate(bits):
+            hist[ti] += np.bincount(track_bits, minlength=NUM_HIST_BINS)
+    return hist
+
+
+def write_target_hist(data_dir, processes=16, block_seqs=256):
+    """Store exact per-track histograms of the stored targets in each zarr.
+
+    Counts are over every position of every sequence in all examples/*.zarr,
+    of the fp16 value y_t + y_pair(t) (strand_pair from targets.txt; unstranded
+    tracks are their own pair, so 2 * y_t), clamped to 65504 and binned by fp16
+    bit pattern. Counts are written to each zarr as ``target_hist``
+    (num_targets, NUM_HIST_BINS) int64,
+    from which SpecPearsonCorrCoef builds its quantile maps.
 
     Args:
-        data_dir: dataset directory with examples/*.zarr.
+        data_dir: dataset directory with targets.txt and examples/*.zarr.
         processes: parallel readers.
         block_seqs: sequences per read task.
 
     Returns:
-        np.ndarray of per-track means.
+        np.ndarray of per-track histograms.
     """
+    targets_df = pd.read_csv(f"{data_dir}/targets.txt", sep="\t", index_col=0)
+    pair = strand_pair_indices(targets_df)
     zarr_files = natsorted(glob.glob(f"{data_dir}/examples/*.zarr"))
     tasks = []
-    num_positions = 0
     for zarr_file in zarr_files:
-        num_seqs, _, length = zarr.open(zarr_file, mode="r")["target"].shape
-        num_positions += num_seqs * length
+        num_seqs = zarr.open(zarr_file, mode="r")["target"].shape[0]
         for start in range(0, num_seqs, block_seqs):
-            tasks.append((zarr_file, start, min(start + block_seqs, num_seqs)))
+            tasks.append((zarr_file, start, min(start + block_seqs, num_seqs), pair))
 
+    hist = 0
     with multiprocessing.get_context("spawn").Pool(processes) as pool:
-        target_sum = np.sum(pool.map(_target_block_sum, tasks), axis=0)
-    target_means = target_sum / num_positions
+        for block_hist in pool.imap_unordered(_target_block_hist, tasks):
+            hist = hist + block_hist.astype(np.int64)
 
     for zarr_file in zarr_files:
-        zarr.open(zarr_file, mode="r+")["target"].attrs["mean"] = target_means.tolist()
-    return target_means
+        root = zarr.open_group(zarr_file, mode="r+")
+        root.create_array(
+            "target_hist",
+            shape=hist.shape,
+            dtype="int64",
+            chunks=(1, NUM_HIST_BINS),
+            overwrite=True,
+        )[:] = hist
+    return hist
 
 
 def targets_prep_strand(targets_df):

src/baskerville/helpers/fold_utils.py

@@ -3,8 +3,10 @@
 Both ``hound_eval_folds`` and ``hound_eval_genes_folds`` need the same logic to
 (a) read the authoritative fold assignment a replicate was trained with,
 (b) decide whether a replicate's weights exist (locally or in GCS marker mode),
-and (c) resolve the ``model_best.pth`` path to pass to the eval command. They're
-collected here so the fold scripts share one implementation.
+(c) resolve the ``model_best.pth`` path to pass to the eval command, and (d)
+link ``test`` to the held-out fold's eval dir. They're collected here so the
+fold scripts (including ``hound_train_folds --eval``) share one
+implementation.
 """
 
 from __future__ import annotations
@@ -60,3 +62,26 @@ def resolve_model_file(
     if gcs_models_dir is not None:
         return f"{GCP_MODELS_MOUNT}/{fold_cross}/train/model_best.pth"
     return f"{container_models_dir}/{fold_cross}/train/model_best.pth"
+
+
+def link_test_metrics(eval_dir, test_fold):
+    """Symlink test -> fold{test_fold} and test.out -> fold{test_fold}.out.
+
+    Relative targets so the links resolve wherever eval_dir lives. Targets may
+    be dangling until eval completes. A stale link (e.g. left by the pre-GCP
+    script that linked test -> fold{fi}, wrong when test_fold != fi for
+    cross > 0) is repointed; a real file/dir at the link path is left alone.
+    """
+    os.makedirs(eval_dir, exist_ok=True)
+    for name, target in (
+        ("test", f"fold{test_fold}"),
+        ("test.out", f"fold{test_fold}.out"),
+    ):
+        link = f"{eval_dir}/{name}"
+        if os.path.islink(link):
+            if os.readlink(link) == target:
+                continue  # already correct
+            os.unlink(link)  # stale link → repoint below
+        elif os.path.lexists(link):
+            continue  # real file/dir squatting the name → don't clobber
+        os.symlink(target, link)