总览 · ← 返回运行 20261002-202908-search-t1-scr-B

节点 n6

ot_cfm 族：whiten-PCA(25) 上 minibatch 精确 OT 配对训练 MLP 速度场，自最后输入阻尼 Euler 外推（α=0.65、t 钳 1），PCA 逆变换+逐细胞残差解码；单输入退化为 copy_last。

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

来自 programs.git a31ffb90e9:solution/METHOD.md

调研员的计划

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

对比：这个提交的上一版（种子程序：相对空仓库）。改动的文件：solution/EXECUTION.json +1 −0、solution/METHOD.md +38 −0、solution/run.py +208 −0

diff --git a/solution/EXECUTION.json b/solution/EXECUTION.jsonnew file mode 100644index 0000000..9d5125c--- /dev/null+++ b/solution/EXECUTION.json@@ -0,0 +1 @@+{"gpu": false}diff --git a/solution/METHOD.md b/solution/METHOD.mdnew file mode 100644index 0000000..9691323--- /dev/null+++ b/solution/METHOD.md@@ -0,0 +1,38 @@+ot_cfm 族：whiten-PCA(25) 上 minibatch 精确 OT 配对训练 MLP 速度场，自最后输入阻尼 Euler 外推（α=0.65、t 钳 1），PCA 逆变换+逐细胞残差解码；单输入退化为 copy_last。++## 方法++1. 读 manifest，取两个输入阶段（`view_io.inputs_by_time`），只保留两阶段都实测的基因（`var["covered"]` 交集）。+2. 训练集每阶段子采样 ≤3000 细胞；按合并方差选 2500 HVG；PCA 25 维，**whiten=True**（消融显示 whiten 关键：cell_state 50.8→54.3）。+3. torchcfm `ExactOptimalTransportConditionalFlowMatcher(sigma=0)`，batch=64，MLP [26→64→64→25]（末层 ×0.01 初始化），Adam lr=1e-3，1500 步，MSE(v(xt,t), ut)。+4. 积分：t 从 1 到 t_end = 1 + (target−last)/(last−first)（由 manifest 时间差算出，视图无关），Euler max(5, ceil(5·span)) 步，速度乘阻尼 α=0.65，网络 t 输入钳制 ≤1.0（训练域外不外推 t）。+5. 解码：`clip(PCA⁻¹(z_end) + resid, 0)`，resid = 该细胞原始 HVG 表达 − PCA⁻¹(z_start)；只在原非零位置写入（不抬正零、保持稀疏），非 HVG 基因原样保留。输出细胞数 = `target_n_cells`（最后阶段细胞的 rng 抽样）。+6. 单输入阶段 → copy_last 退路。未用 external/ 与 prior/。++## 关键参数（X3 A 半查分，10 次）++- copy_last（=ablation）47.62；α=0.75 无 whiten 47.18；α=0.35 无 whiten 47.77；whiten α=0.25/0.35/0.5/0.65/0.75/0.85 = 48.41/48.62/48.77/**50.41**/49.50/48.43。默认锁定 whiten + α=0.65 + clamp_t=1.0；固定 numpy 种子后复跑 50.29（OT 配对序列略变，差异在噪声内）。++## 机制生效证据（PLAN mechanism_evidence）++- 速度方向：100 个 E9.0 细胞探针两两余弦均值 ≈ 0.42（<0.95，非常向量）；|v| 均值 4.09、std 1.32（状态依赖）。+- 位移：白化空间逐细胞位移 mean 3.4、std 1.15、max 7.9，std 显著 >0。+- 四组分 vs copy_last(47.62)：cell_state 49.5→59.6，de_recovery 43.1→45.7，direction 49.1→49.1（持平），covariation 48.6→43.6（受损，见弱点）。+- 输出方差/零结构：nnz/cell 1967 vs 输入 1974，基因方差均值 0.105 vs 0.101，无普遍扩散迹象。++## 关闭机制对照（mechanism_off_control）++`--ablation zero_velocity` 将 α 置 0，积分不动，输出 = 最后阶段抽样细胞的复制（PCA 重建+残差在 α=0 时恒等于原值）。X3 查分 47.62，与打开时 50.29 差 +2.67，机制生效。提交默认打开。++## 确定性修复++torchcfm `OTPlanSampler.sample_plan` 内部用 numpy 全局 RNG（`np.random.choice`），导致跨进程不确定；已在训练前 `np.random.seed(seed)`，并 `torch.set_num_threads(1)` + `use_deterministic_algorithms(True)`。同 seed 两次运行输出逐比特相同（已验证）。seed 7 也跑通。++## 验证过 / 未验证++- 已验证：X3 视图（两输入，目标跨 2 个区间）跑通、vec-check ok、~35s（< 30min 限）、内存远小于 28GB、确定性、α/whiten 消融。+- 未验证：final 视图（E8.5+E9.5→E10.5，span=1）与 proxy（单输入退化路径只做过代码路径推演，未实跑）；covariation 受损是否可在后续用位移的邻域平滑/按 PC 分层阻尼修复；训练步数、HVG 数、PCA 维度的进一步搜索。++## 知识来源++未使用任何保留阶段/保留基因型的测量信息或文献数值；只用了通用做法（OT-CFM, Tong et al. 2023 算法本身，经 torchcfm 实现）。α、whiten、clamp_t 全部由 X3 A 半查分选出（视图内数据驱动，不依赖绝对时间或视图身份）。diff --git a/solution/run.py b/solution/run.pynew file mode 100644index 0000000..a9974c2--- /dev/null+++ b/solution/run.py@@ -0,0 +1,208 @@+"""OT-CFM small velocity field: minibatch-OT paired training, damped Euler extrapolation in PCA space.++Family: ot_cfm (PLAN). Mechanism: learn v(z,t) on shared 25-d PCA space with exact-OT minibatch+pairing of stage0->stage1 cells; extrapolate last-stage cells from t=1 to t_end (data-driven from+manifest times) with damped velocity (alpha); decode via PCA inverse + per-cell residual.+Ablation: --ablation zero_velocity zeroes the velocity field -> output equals copy of last stage.+Single input stage -> copy_last fallback. CPU only, deterministic given --seed.+"""+from __future__ import annotations++import argparse+import sys+from pathlib import Path++import numpy as np+from scipy import sparse++from src.task1_temporal import view_io+++def log(*a):+    print(*a, file=sys.stderr, flush=True)+++def copy_last(view, manifest, genes, seed, out):+    rng = np.random.default_rng(seed)+    entries = view_io.inputs_by_time(manifest)+    adata = view_io.read_stage(view, entries[-1], genes, missing="fill")+    n_out = view_io.target_n_cells(manifest, adata.n_obs)+    idx = view_io.sample_rows(adata.n_obs, n_out, rng)+    view_io.write_prediction(adata.X[idx], genes, out, seed=seed)+++def main():+    ap = argparse.ArgumentParser()+    ap.add_argument("--data", required=True)+    ap.add_argument("--out", required=True)+    ap.add_argument("--seed", type=int, required=True)+    ap.add_argument("--ablation", default=None, choices=[None, "zero_velocity"])+    ap.add_argument("--alpha", type=float, default=0.65)+    ap.add_argument("--dim", type=int, default=25)+    ap.add_argument("--steps", type=int, default=1500)+    ap.add_argument("--clamp-t", type=float, default=1.0)+    ap.add_argument("--whiten", action="store_true", default=True)+    ap.add_argument("--no-whiten", dest="whiten", action="store_false")+    args = ap.parse_args()++    view = args.data+    manifest = view_io.load_manifest(view)+    genes = view_io.panel_genes(view, manifest)+    entries = view_io.inputs_by_time(manifest)+    if len(entries) < 2:+        log("single input stage -> copy_last fallback")+        copy_last(view, manifest, genes, args.seed, args.out)+        return++    rng = np.random.default_rng(args.seed)+    import torch++    torch.manual_seed(args.seed)+    torch.set_num_threads(1)+    torch.use_deterministic_algorithms(True)++    a0 = view_io.read_stage(view, entries[0], genes, missing="fill")+    a1 = view_io.read_stage(view, entries[-1], genes, missing="fill")+    t0 = float(entries[0]["time"])+    t1 = float(entries[-1]["time"])+    tT = float(manifest["target"]["time"])+    dt = t1 - t0+    t_end = 1.0 + (tT - t1) / dt if dt > 0 else 1.0+    log(f"times {t0}->{t1} target {tT}; t_end={t_end:.3f}")++    cov = np.asarray(a0.var["covered"], dtype=bool) & np.asarray(a1.var["covered"], dtype=bool)+    cov_idx = np.flatnonzero(cov)+    log(f"covered in both stages: {cov.size} -> {cov_idx.size}")++    X0 = np.asarray(a0.X[:, cov_idx].todense(), dtype=np.float32)+    X1 = np.asarray(a1.X[:, cov_idx].todense(), dtype=np.float32)+    del a0++    # subsample for training+    cap = 3000+    if X0.shape[0] > cap:+        X0 = X0[np.sort(rng.choice(X0.shape[0], cap, replace=False))]+    if X1.shape[0] > cap:+        keep1 = np.sort(rng.choice(X1.shape[0], cap, replace=False))+    else:+        keep1 = np.arange(X1.shape[0])++    # HVG among covered genes: top variance on pooled cells+    pooled = np.concatenate([X0, X1[keep1]], axis=0)+    var = pooled.var(axis=0)+    n_hvg = min(2500, pooled.shape[1])+    hvg = cov_idx[np.argsort(-var)[:n_hvg]]+    hvg_local = np.argsort(-var)[:n_hvg]+    del pooled+    X0h = X0[:, hvg_local]+    X1h = X1[keep1][:, hvg_local]++    # PCA on union+    from sklearn.decomposition import PCA++    n_comp = int(min(args.dim, X0h.shape[0] - 1, X1h.shape[0] - 1, n_hvg))+    Z = np.concatenate([X0h, X1h], axis=0)+    pca = PCA(n_components=n_comp, random_state=args.seed, whiten=args.whiten)+    pca.fit(Z)+    z0 = pca.transform(X0h).astype(np.float32)+    z1_all = pca.transform(X1[:, hvg_local]).astype(np.float32)+    z1 = z1_all[keep1]+    del Z++    # ---- train OT-CFM velocity field ----+    from torchcfm import ExactOptimalTransportConditionalFlowMatcher++    np.random.seed(args.seed)  # torchcfm OTPlanSampler uses numpy global RNG++    fm = ExactOptimalTransportConditionalFlowMatcher(sigma=0.0)+    dim = z0.shape[1]+    torch.manual_seed(args.seed)+    net = torch.nn.Sequential(+        torch.nn.Linear(dim + 1, 64), torch.nn.ReLU(),+        torch.nn.Linear(64, 64), torch.nn.ReLU(),+        torch.nn.Linear(64, dim),+    )+    with torch.no_grad():+        net[-1].weight.mul_(0.01)+        net[-1].bias.mul_(0.01)+    opt = torch.optim.Adam(net.parameters(), lr=1e-3)+    T0 = torch.from_numpy(z0)+    T1 = torch.from_numpy(z1)+    batch = 64+    n0, n1 = T0.shape[0], T1.shape[0]+    g = torch.Generator().manual_seed(args.seed)+    net.train()+    for step in range(args.steps):+        i0 = torch.randint(0, n0, (batch,), generator=g)+        i1 = torch.randint(0, n1, (batch,), generator=g)+        t, xt, ut = fm.sample_location_and_conditional_flow(T0[i0], T1[i1])+        v = net(torch.cat([xt, t[:, None]], dim=1))+        loss = torch.nn.functional.mse_loss(v, ut)+        opt.zero_grad(set_to_none=True)+        loss.backward()+        opt.step()+        if step % 300 == 0 or step == args.steps - 1:+            log(f"step {step} loss {loss.item():.4f}")+    net.eval()++    # ---- mechanism evidence ----+    with torch.no_grad():+        probe = T1[torch.arange(0, n1, max(1, n1 // 100))][:100]+        tp = torch.full((probe.shape[0],), 1.0)+        vp = net(torch.cat([probe, tp[:, None]], dim=1)).numpy()+        norms = np.linalg.norm(vp, axis=1)+        vcos = []+        for i in range(len(vp)):+            c = (vp @ vp[i]) / np.maximum(norms * norms[i], 1e-12)+            vcos.append(np.mean(np.delete(c, i)))+        log(f"velocity: mean pairwise cos {np.mean(vcos):.3f} (degenerate if >0.95), "+            f"|v| mean {norms.mean():.3f} std {norms.std():.3f}")++    # ---- output cells: sample from last stage ----+    n_out = view_io.target_n_cells(manifest, a1.n_obs)+    out_idx = view_io.sample_rows(a1.n_obs, n_out, rng)+    z_start = z1_all[out_idx]+    resid_h = X1[out_idx][:, hvg_local] - pca.inverse_transform(z_start).astype(np.float32)++    # ---- damped Euler integration t: 1 -> t_end ----+    alpha = args.alpha+    if args.ablation == "zero_velocity":+        alpha = 0.0+    span = t_end - 1.0+    n_steps = max(5, int(np.ceil(5 * span)))+    z = z_start.copy()+    h = span / n_steps+    with torch.no_grad():+        for s in range(n_steps):+            tt = 1.0 + s * h+            if args.clamp_t is not None:+                tt = min(tt, args.clamp_t)+            zt = torch.from_numpy(z)+            v = net(torch.cat([zt, torch.full((z.shape[0], 1), tt, dtype=torch.float32)], dim=1)).numpy()+            z = z + alpha * v * h+    disp = np.linalg.norm(z - z_start, axis=1)+    log(f"displacement: mean {disp.mean():.3f} std {disp.std():.3f} max {disp.max():.3f} "+        f"(z scale ~1); steps {n_steps} alpha {alpha}")++    # ---- decode ----+    rec = pca.inverse_transform(z).astype(np.float32)+    pred_h = np.clip(rec + resid_h, 0.0, None)++    Xout = a1.X[out_idx].tocsr().astype(np.float32)+    Xd = np.asarray(Xout.todense(), dtype=np.float32)+    orig_h = Xd[:, hvg]+    new_h = pred_h+    # keep sparsity: do not lift original zeros to positive+    new_h = np.where(orig_h > 0, new_h, 0.0)+    new_h[new_h < 1e-3] = 0.0+    Xd[:, hvg] = new_h+    Xfin = sparse.csr_matrix(Xd)+    log(f"output {Xfin.shape}, nnz/cell mean {Xfin.nnz / Xfin.shape[0]:.1f} "+        f"(input last stage {a1.X.nnz / a1.n_obs:.1f}), var mean {np.var(Xd, axis=0).mean():.4f} "+        f"vs input {np.var(np.asarray(a1.X[out_idx].todense()), axis=0).mean():.4f}")+    view_io.write_prediction(Xfin, genes, args.out, seed=args.seed)+    log("done")+++if __name__ == "__main__":+    main()

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

用到的知识库条目

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

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

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

分析员

ANALYSIS.json written

工程师

完成。

**做了什么**：从零实现 PLAN 的 ot_cfm 方法族 —— 共享 whiten-PCA(25) 空间上用 torchcfm 精确 OT minibatch 配对训练小 MLP 速度场（1500 步），从最后输入阶段做阻尼 Euler 外推（α=0.65，网络 t 输入钳制在训练域 ≤1，t_end 由 manifest 时间差算出，视图无关），解码 = PCA 逆变换 + 逐细胞残差、保持稀疏。单输入退化为 copy_last；`--ablation zero_velocity` 对照已实现（关闭时 = copy_last）。修复了 torchcfm 内部 numpy 全局 RNG 导致的跨进程不确定性（`np.random.seed` + 单线程 + deterministic algorithms，已验证逐比特复现）。

**查分结果**（X3 A 半，10/20 次）：copy_last 对照 47.62；无 whiten 各 α 均 ≤47.8；whiten 后 α 扫描 0.25→0.85 = 48.41/48.62/48.77/**50.41**/49.50/48.43，锁定 α=0.65；确定性修复后默认配置复跑 **50.29**（cell_state 59.6、de_recovery 45.7、direction 49.1、covariation 43.6）。机制证据：速度两两余弦 0.42（非常向量）、位移 std 1.15 > 0。

**下一步最值得试**：covariation 是唯一受损组分（48.6→43.6）——可对位移场做邻域平滑或按 PC 序分层阻尼（低阶 PC 大步、高阶 PC 小步），有望在保住 cell_state 增益的同时收回 2–3 分；其次在 final 视图形态（span=1）上验证 α=0.65 是否仍最优。

调研员

e-3, 训练 1500 步（约 2-3 分钟）\n   - 若 torchcfm 不可用则手写：t~U(0,1), xt=(1-t)x0+t*x1, ut=x1-x0, 用 POT ot.emd 做 batch 内配对\n4. 推理/积分：\n   - 对 E9.5 全部细胞，在 PCA 空间从 t=1 积分到 t=2\n   - 固定 Euler，5 步，每步乘以阻尼 α=0.75（即有效外推 0.75 天而非 1 天）\n   - 外推时速度场仍用训练好的 v(x,t)，t 输入范围 [1,2]（训练只覆盖 [0,1]，外推依赖网络泛化）\n5. 解码：\n   - 用 PCA components_ 投影回基因空间得到预测表达矩阵的 PCA 重建部分\n   - 对每个细胞加回其自身的残差（原始表达 - PCA 重建），保持细胞间协变结构\n   - clip 负值为 0，保持稀疏性（不将零抬正）\n6. 关闭机制对照（ablation）：\n   - 同一程序 --velocity zero：将速度场输出置零再积分 → 输出应等于 copy_last（E9.5 原样输出）\n   - 预期差别：打开机制后预测应偏离 E9.5，方向组与 DE 组应有变化；若输出与 copy_last 逐比特相同则机制未生效\n7. vec-score 查分：先在小规模（500 细胞、500 步训练）跑通确认无报错，再全量训练查分。若查分 < 48（低于 copy_last 47.92 + 噪声），检查速度场是否退化为常向量。\n8. 关键参数搜索范围（仅在时间允许时）：α ∈ {0.5, 0.75, 1.0}；PCA 维度 ∈ {20, 25, 30}；训练步数 ∈ {1000, 1500, 2000}。优先固定 α=0.75, dim=25, steps=1500 跑一次。",
  "expected_groups": ["de_recovery", "direction", "cell_state"],
  "risks": "1. 两个时间点只提供一个区间的配对，速度场在 t>1 外推区域无训练信号，可能过冲或方向错误 → 用阻尼 α=0.75 缓解；Engineer 应检查预测的方差是否显著大于输入（普遍扩散迹象）。2. 速度场可能退化为近似常向量（每个区域方向相同）→ 报告速度在 PCA 空间不同区域的方向余弦相似度分布；若 >0.95 则退化为常数位移。3. torchcfm 的 ExactOTCFM 对 batch 内配对用精确 OT，若 batch 太小（<32）配对噪声大 → 设 batch_size=64。4. 30 分钟时限：训练本身很快（<5 分钟），风险在数据加载和 PCA；若 HVG 选择耗时可固定取方差前 2500 基因。5. X3 是心脏数据、可能细胞数较少 → 若某阶段 <500 细胞则不子采样全用。",
  "family_id": "ot_cfm",
  "mechanism": "在共享 PCA 空间中，用 minibatch 精确 OT 配对 E8.5→E9.5 细胞对，训练一个时间条件的小 MLP 速度场 v(x,t)；推理时从 E9.5 细胞出发做阻尼 Euler 积分外推，解码时加回逐细胞残差。",
  "vs_constant_shift": "常数位移对同一类型/簇的所有细胞施加相同向量；OT-CFM 的速度是细胞状态 x 的非线性函数——即使同类型内，处于不同 PCA 位置的细胞获得不同方向和大小的速度。此外速度随积分步更新（x 变了 → v 变了），不是固定一步位移。",
  "mechanism_evidence": "Engineer 应报告：(1) 速度场在 25 维 PCA 空间中不同位置的方向分布——计算 100 个随机细胞对的余弦相似度均值，若 >0.95 则退化为常向量；(2) 位移幅度的细胞间标准差（应显著大于 0）；(3) 四组分各自相对 copy_last 的变化量；(4) 预测分布的方差/零值比例与 E9.5 输入的对比（检测普遍扩散）。",
  "mechanism_off_control": "同一 run.py 加 --ablation zero_velocity 参数：将训练好的速度网络输出乘以 0 再积分。此时 Euler 积分不改变起点，输出应逐比特等于 copy_last（E9.5 原样复制加残差）。若打开/关闭输出有任何差异，说明有非速度场的副作用需排查。预期差别：打开时四组分中 direction 和 de_recovery 应有 ≥2 分变化，关闭时回落到 copy_last 水平（~47-48）。",
  "sources": []
}

审查员

REVIEW.json written