总览 · ← 返回运行 20261001-032558-search-t1-spark-c3

节点 n6

按类型单快照伪时间（根=最增殖细胞）回归逐基因斜率并 EB 收缩，以线性空间倍数 exp(s·slope·w_g)（s=0.15）乘法位移抽样细胞、零结构保留；两阶段时改用 α=0.25 收缩 delta 并按与斜率的余弦门控；增殖重加权 γ 测无收益置 0。

方法说明?节点程序自带的 METHOD.md：这个程序做了什么、为什么。

来自 programs.git 7efa4b504b:solution/METHOD.md

调研员的计划

代码改动?这个节点的程序和父节点程序的逐行差别：绿色是新增，红色是删除。

对比：父节点版本 afbc939bfc。改动的文件：solution/METHOD.md +31 −0、solution/run.py +249 −11

diff --git a/solution/METHOD.md b/solution/METHOD.mdnew file mode 100644index 0000000..acb3d50--- /dev/null+++ b/solution/METHOD.md@@ -0,0 +1,31 @@+按类型单快照伪时间（根=最增殖细胞）回归逐基因斜率并 EB 收缩，以线性空间倍数 exp(s·slope·w_g)（s=0.15）乘法位移抽样细胞、零结构保留；两阶段时改用 α=0.25 收缩 delta 并按与斜率的余弦门控；增殖重加权 γ 测无收益置 0。++# 生长率重加权 + 单快照伪时间乘法位移（EB 收缩）+ 收缩两阶段 delta++单阶段（proxy）：对每个细胞类型在 SVD 坐标上建 kNN 图，以最增殖细胞为根做 dijkstra 伪时间，逐基因回归斜率并按 t²/(1+t²) 经验贝叶斯收缩，再以线性空间倍数 exp(s·slope·w_g)（s=0.15，zeros 保持为零）小步位移抽样细胞。两阶段（final）：类型伪批量 delta 做逐基因 EB 收缩并乘 α=0.25，与伪时间斜率方向余弦门控（cos<0.3 减半，cos<0 或新类型改用 s·slope）。增殖重加权 γ 默认 0（未测出收益，代码保留）。++## 方法细节++- **P0 生长率重加权（γ，默认 0）**：从 view `prior/reactome/gene_sets.gmt` 取名字含 "Cell Cycle"、大小 10–500 的集合并 ∩ 面板（读到 302 基因；读不到时退回通用细胞周期符号表）。每类型分 = ≤2000 随机细胞的 log1p 均值；w_c = clip((p_c/median)^γ, 0.2, 5)，按 w 无放回抽样。γ=0 精确恢复父行为（sample_rows）。+- **P1 单快照伪时间位移（s，默认 0.15，乘法模式）**：每类型（n≥50，共 18/18）：TruncatedSVD 30 维 → kNN(k=15) 高斯核（带宽=中位距离）→ 对称化 → dijkstra（权=1/(aff+1e-8)），根 = 类型内增殖分最高细胞；dpt 按 p99 归一。slope_g = cov(x_g,dpt)/var(dpt)（稀疏矩量法，≤4000 细胞），EB 权重 w_g = t²/(1+t²)，t² = n·r²/(1−r²)。**乘法应用**：x ← log1p(expm1(x)·exp(clip(s·slope·w_g,±3)))，零结构完全保留（加法模式会稠密化、covariation 掉 8.6 分，已否决）。+- **P2 两阶段分支（final，α=0.25）**：type_deltas 逐基因 EB（se² = 两阶段合并类型内方差/n）；cos_c = cos(delta_c, slope_c)：cos≥0.3 → α·delta·wδ；0≤cos<0.3 → 0.5α·delta·wδ；cos<0 或 prev 缺该类型 → s·slope·w。α 未在 proxy 上调（proxy 无两阶段），取 k018 教训（α=1 → 48.6 < 地板）的保守先验。++## 已验证（T1:val proxy，9 次查分）++| 配置 | 榜分 | 备注 |+|---|---|---|+| γ0 s0（=父） | 49.77 | 理智检查通过 |+| γ0.5 s0 / γ1 s0 | 49.84 / 49.00 | 增殖重加权无效，弃 |+| γ0 s0.25 加法 / s0.5 加法 | 49.24 / 43.37 | direction +4.7 但 covariation −8.6，稠密化致命 |+| γ0 s0.25 乘法 | 51.43 / 50.98 (seed0/1) | 四组全升 |+| **γ0 s0.15 乘法（默认）** | **51.58 / 50.90 (seed0/1)** | 均值 51.24，最优 |+| γ0.5 s0.25 乘法 | 51.33 | cov +1.5 但 cs −1.4，不加 γ |++- 运行 ~85s、峰值 ~2.5GB（限 30min/28GB）。final 两阶段分支（EB delta、cos 门控、乘法应用）已用合成数据单元自测：无 NaN、形状正确、零结构保留；proxy 上无法实测该分支。+- 确定性：仅用 np.random.default_rng(seed)；TruncatedSVD 固定 random_state=seed。++## 未验证 / 风险++- final 端 P2 的 α=0.25、cos 门控阈值均未在真实两阶段数据上验证。+- DPT 轴可能被细胞周期主导（根 = 最增殖细胞会强化）；s 已保守取 0.15，且乘法模式保零结构，covariation/cell_state 实测未受损。+- 组成（cell_state）仍是弱组（49.3 vs 节点 2 的 56.9）：增殖重加权在 proxy 上无效，说明 E8.5→E9.5 的组成变化不是增殖驱动的，需要别的信号（如节点 2 的组织学先验，但须从数据导出）。diff --git a/solution/run.py b/solution/run.pyindex f3a0f25..3051539 100644--- a/solution/run.py+++ b/solution/run.py@@ -1,22 +1,34 @@ #!/usr/bin/env python3-"""pseudobulk_shift: latest stage + per-cell-type pseudobulk delta of the last step.+"""Growth-rate reweighted sampling + single-snapshot pseudotime shift. -The delta is mean(last|type) - mean(prev|type) over the two latest inputs,-computed on the full stages and added once to a subsample of the latest stage-(clipped at 0). Types missing from the earlier stage are copied unchanged.+P0: per-cell proliferation score from a reactome "Cell Cycle" gene set found in+    the view's prior/ (fallback: generic cell-cycle symbols). Per-type weight+    w_c = clip((p_c / median(p))**gamma, 0.2, 5); cells sampled without+    replacement with prob proportional to w_c. gamma=0 -> uniform (parent).+P1: per cell type, build a kNN graph on truncated-SVD coords, diffusion-style+    pseudotime (dijkstra) rooted at the most proliferative cell, per-gene slope+    with empirical-Bayes shrinkage w_g = t^2/(1+t^2); shift cells by+    s * slope * w_g (clipped at 0). s=0 -> no shift.+P2 (>=2 input stages): per-type pseudobulk delta with per-gene EB shrinkage,+    scaled by alpha=0.25 (or 0.5*alpha when delta/slope directions disagree,+    cos<0.3); types with cos<0 or missing from prev get the P1 slope shift. -With a single input stage (T1 proxy: E8.5 only) there is no step to take a-delta from, so this falls back to copy_last with the same sampling. The proxy-therefore cannot tell this seed from copy_last; that gap is expected.+All biology inputs come from the view data + view prior/; nothing stage-+specific is hardcoded. Deterministic given --seed. """  from __future__ import annotations  import argparse+from pathlib import Path  import numpy as np+from scipy import sparse+from scipy.sparse import csgraph+from scipy.spatial import cKDTree+from sklearn.decomposition import TruncatedSVD -from src.task1_temporal.baselines import shift_rows, type_deltas+from src.task1_temporal.baselines import type_deltas from src.task1_temporal.view_io import (     inputs_by_time,     labels_of,@@ -28,12 +40,197 @@ from src.task1_temporal.view_io import (     write_prediction, ) +FALLBACK_CYCLE = [+    "Mki67", "Top2a", "Cdk1", "Ccnb1", "Ccna2", "Aurka", "Aurkb", "Plk1",+    "Bub1", "Rrm2", "Pcna", "Mcm2", "Mcm3", "Mcm4", "Mcm5", "Mcm6", "Mcm7",+    "Cdc20", "Kif11",+]+GAMMA = 0.0+S_SHIFT = 0.15+ALPHA = 0.25+W_LO, W_HI = 0.2, 5.0+++def col_stats(X: sparse.csr_matrix):+    """(mean_g, var_g, n) with unbiased var; NaN -> 0."""+    n = X.shape[0]+    s1 = np.asarray(X.sum(axis=0)).ravel().astype(np.float64)+    sq = X.copy()+    sq.data = sq.data.astype(np.float64) ** 2+    s2 = np.asarray(sq.sum(axis=0)).ravel()+    mean = s1 / max(n, 1)+    var = (s2 - n * mean**2) / max(n - 1, 1)+    var = np.where(np.isfinite(var), var, 0.0)+    var = np.maximum(var, 0.0)+    return mean, var, n+++def load_cycle_genes(view, manifest, panel) -> np.ndarray:+    idx_of = {g: i for i, g in enumerate(panel)}+    gmt = None+    for e in manifest.get("prior", []) or []:+        if e.get("id") == "reactome":+            p = Path(view) / e["path"] / "gene_sets.gmt"+            if p.exists():+                gmt = p+                break+    found: set[str] = set()+    if gmt is not None:+        for line in gmt.read_text(encoding="utf-8").splitlines():+            parts = line.rstrip("\n").split("\t")+            if len(parts) < 3:+                continue+            name = parts[1].lower()+            if "cell cycle" not in name:+                continue+            members = {m for m in parts[2:] if m}+            if not (10 <= len(members) <= 500):+                continue+            found |= {m for m in members if m in idx_of}+    if not found:+        found = {g for g in FALLBACK_CYCLE if g in idx_of}+    if not found:+        return np.zeros(0, dtype=np.int64)+    return np.sort(np.fromiter((idx_of[g] for g in found), dtype=np.int64))+++def prolif_scores(X: sparse.csr_matrix, cycle_idx: np.ndarray) -> np.ndarray:+    if cycle_idx.size == 0:+        return np.zeros(X.shape[0])+    sub = X[:, cycle_idx]+    return np.asarray(sub.sum(axis=1)).ravel() / float(cycle_idx.size)+++def type_weights(labels: np.ndarray, scores: np.ndarray, gamma: float, rng) -> np.ndarray:+    types = np.unique(labels)+    p = np.empty(types.size)+    for i, t in enumerate(types):+        idx = np.flatnonzero(labels == t)+        if idx.size > 2000:+            idx = rng.choice(idx, size=2000, replace=False)+        p[i] = scores[idx].mean()+    med = np.median(p)+    if not np.isfinite(med) or med <= 1e-12:+        w = np.ones(types.size)+    else:+        w = np.clip((np.maximum(p, 1e-12) / med) ** gamma, W_LO, W_HI)+        w = np.where(np.isfinite(w), w, 1.0)+    t_to_w = {str(t): w[i] for i, t in enumerate(types)}+    return np.array([t_to_w[str(l)] for l in labels])+++def type_slopes(X: sparse.csr_matrix, labels: np.ndarray, scores: np.ndarray,+                rng, seed: int, max_cells: int = 4000):+    """slope_g * w_g (EB-shrunk) per type, from single-snapshot pseudotime."""+    out: dict[str, np.ndarray] = {}+    G = X.shape[1]+    for t in np.unique(labels):+        idx = np.flatnonzero(labels == t)+        n_c = idx.size+        if n_c < 50:+            continue+        Xc = X[idx]+        sc = scores[idx]+        if n_c > max_cells:+            sub = rng.choice(n_c, size=max_cells, replace=False)+            Xc, sc = Xc[sub], sc[sub]+        n = Xc.shape[0]+        comps = int(min(30, n - 2, Xc.shape[1] - 1))+        if comps < 2:+            continue+        coords = TruncatedSVD(n_components=comps, random_state=seed).fit_transform(Xc)+        k = int(min(15, n - 1))+        dist, nbr = cKDTree(coords).query(coords, k=k + 1)+        d, nb = dist[:, 1:], nbr[:, 1:]+        sigma = np.median(d)+        if not np.isfinite(sigma) or sigma <= 1e-12:+            continue+        aff = np.exp(-(d**2) / (2.0 * sigma**2))+        rows = np.repeat(np.arange(n), k)+        A = sparse.coo_matrix((aff.ravel(), (rows, nb.ravel())), shape=(n, n)).tocsr()+        A = A.maximum(A.T)+        A.data = 1.0 / (A.data + 1e-8)+        root = int(np.argmax(sc))+        dpt = csgraph.dijkstra(A, directed=False, indices=root)+        fin = np.isfinite(dpt)+        if fin.sum() < 30:+            continue+        tv = dpt[fin]+        p99 = np.percentile(tv, 99)+        if p99 <= 1e-12:+            continue+        tv = np.clip(tv / p99, 0.0, 1.0)+        var_t = tv.var()+        if var_t <= 1e-8:+            continue+        nf = int(fin.sum())+        tm = tv - tv.mean()+        Xf = Xc[fin]+        mean_g, var_g, _ = col_stats(Xf)+        cov = (np.asarray(Xf.T.dot(tm)).ravel() - nf * mean_g * tm.mean()) / (nf - 1)+        slope = cov / var_t+        with np.errstate(divide="ignore", invalid="ignore"):+            r2 = cov**2 / np.maximum(var_g * var_t, 1e-30)+        r2 = np.clip(np.where(np.isfinite(r2), r2, 0.0), 0.0, 0.999999)+        t2 = nf * r2 / (1.0 - r2)+        w_g = t2 / (1.0 + t2)+        v = slope * w_g+        v = np.where(np.isfinite(v), v, 0.0).astype(np.float32)+        out[str(t)] = v+    return out+++def apply_shifts(X: sparse.csr_matrix, labels: np.ndarray,+                 shifts: dict[str, np.ndarray], mult: bool = True) -> sparse.csr_matrix:+    blocks, order = [], []+    for t in np.unique(labels):+        idx = np.flatnonzero(labels == t)+        order.append(idx)+        v = shifts.get(str(t))+        if v is not None and np.any(v != 0):+            D = np.asarray(X[idx].todense(), dtype=np.float64)+            if mult:+                fold = np.exp(np.clip(v.astype(np.float64), -3.0, 3.0))+                np.expm1(D, out=D)+                D *= fold+                np.log1p(D, out=D)+            else:+                D += v.astype(np.float64)+            np.maximum(D, 0.0, out=D)+            blocks.append(sparse.csr_matrix(D.astype(np.float32)))+        else:+            blocks.append(X[idx])+    out = sparse.vstack(blocks, format="csr")+    inv = np.empty(X.shape[0], dtype=np.int64)+    inv[np.concatenate(order)] = np.arange(X.shape[0])+    return out[inv]+++def eb_delta_weights(deltas: dict[str, np.ndarray], prev_X, prev_lab, last_X, last_lab):+    """Per-gene EB shrinkage w = t^2/(1+t^2) for each two-stage delta."""+    out = {}+    for t, d in deltas.items():+        ip = np.flatnonzero(prev_lab == t)+        il = np.flatnonzero(last_lab == t)+        _, vp, np_ = col_stats(prev_X[ip])+        _, vl, nl_ = col_stats(last_X[il])+        se2 = vp / max(np_, 1) + vl / max(nl_, 1)+        with np.errstate(divide="ignore", invalid="ignore"):+            t2 = d.astype(np.float64) ** 2 / np.maximum(se2, 1e-30)+        t2 = np.where(np.isfinite(t2), t2, 0.0)+        out[t] = (t2 / (1.0 + t2)).astype(np.float32)+    return out+  def main() -> None:     parser = argparse.ArgumentParser()     parser.add_argument("--data", required=True)     parser.add_argument("--out", required=True)     parser.add_argument("--seed", type=int, default=0)+    parser.add_argument("--gamma", type=float, default=GAMMA)+    parser.add_argument("--s", type=float, default=S_SHIFT)+    parser.add_argument("--alpha", type=float, default=ALPHA)+    parser.add_argument("--mult", type=int, default=1)     args = parser.parse_args()      manifest = load_manifest(args.data)@@ -41,13 +238,54 @@ def main() -> None:     stages = inputs_by_time(manifest)     last = read_stage(args.data, stages[-1], genes)     rng = np.random.default_rng(args.seed)-    rows = sample_rows(last.n_obs, target_n_cells(manifest, last.n_obs), rng)+    labels = labels_of(last)++    cycle_idx = load_cycle_genes(args.data, manifest, genes)+    scores = prolif_scores(last.X, cycle_idx)++    n_out = target_n_cells(manifest, last.n_obs)+    if args.gamma == 0:+        rows = sample_rows(last.n_obs, n_out, rng)+    else:+        w_cells = type_weights(labels, scores, args.gamma, rng)+        p = w_cells / w_cells.sum()+        rows = np.sort(rng.choice(last.n_obs, size=n_out, replace=False, p=p))     X = last.X[rows]+    lab = labels[rows]++    slopes = {}+    if args.s > 0 or len(stages) >= 2:+        slopes = type_slopes(last.X, labels, scores, rng, args.seed)++    shifts: dict[str, np.ndarray] = {}     if len(stages) >= 2:         prev = read_stage(args.data, stages[-2], genes)-        deltas = type_deltas(prev.X, labels_of(prev), last.X, labels_of(last))+        prev_lab = labels_of(prev)+        deltas = type_deltas(prev.X, prev_lab, last.X, labels)+        wd = eb_delta_weights(deltas, prev.X, prev_lab, last.X, labels)         del prev-        X = shift_rows(X, labels_of(last)[rows], deltas)+        for t in set(list(deltas.keys()) + list(slopes.keys())):+            d, sl = deltas.get(t), slopes.get(t)+            if d is not None:+                cos = -1.0+                if sl is not None:+                    nd, ns = np.linalg.norm(d), np.linalg.norm(sl)+                    if nd > 0 and ns > 0:+                        cos = float(d @ sl / (nd * ns))+                if cos >= 0.3:+                    shifts[t] = (args.alpha * d * wd[t]).astype(np.float32)+                elif cos >= 0.0:+                    shifts[t] = (0.5 * args.alpha * d * wd[t]).astype(np.float32)+                elif sl is not None:+                    shifts[t] = (args.s * sl).astype(np.float32)+            elif sl is not None:+                shifts[t] = (args.s * sl).astype(np.float32)+    else:+        for t, sl in slopes.items():+            shifts[t] = (args.s * sl).astype(np.float32)++    if shifts:+        X = apply_shifts(X, lab, shifts, mult=bool(args.mult))     write_prediction(X, genes, args.out, seed=args.seed)  

调研来源?调研员查到并用到的知识条目和文献检索结果（只列标题和编号）。

用到的知识库条目

分析结果?分析员写的 ANALYSIS.json：改了什么、各组分数怎么变、假设是否成立、经验和下一步建议。

对话摘要?每个角色和大模型对话的统计：轮数、工具调用、用时、token 数和最后的回答摘录；原始记录只给路径。

只给统计和最后回答的摘录；完整对话请到原始记录位置里列出的文件看。

分析员

ANALYSIS.json written

工程师

**做了什么**：在父节点 pseudobulk_shift 上实现了 P0 增殖重加权 + P1 单快照伪时间位移 + P2 两阶段收缩 delta（α=0.25，cos 门控，final 分支已合成数据单测）。关键发现：加法位移稠密化毁掉 covariation（−8.6），改为线性空间倍数乘法 exp(s·slope·w_g) 后零结构保留、四组全升；增殖重加权 γ 在 proxy 无效（49.84/49.00），置 0 保留代码。

**查分**（11/20 次）：默认配置 γ0 s0.15 乘法 = **51.58（seed0）/ 50.90（seed1）**，均值 51.24 vs 父 49.77；s0.25 均值 51.21 打平。方向组 50.16→52.5，de_recovery 50.0→52.5，covariation 51.5→52.6，cell_state 48.1→49.3。默认输出与已打分文件逐字节一致，vec-check 通过。

**下一步最值得试**：cell_state 仍是最大弱组（49.3 vs 节点 2 的 56.9）——组成杠杆未被数据驱动方法捕获，建议从 E8.5 内部信号（如类型间增殖/分化梯度、邻域结构）导出组成先验并与节点 2 的手调权重对照；以及在 final 视图上实测 P2 分支的 α 与 cos 门控。

调研员

PLAN.json written