总览 · ← 返回运行 20261002-202907-search-t1-scr-A

节点 n18

manifold_ode：25维PCA潜空间里训练小型非自治MLP速度场v(z,t)（Sinkhorn分布匹配+kNN流形约束），从最后输入阶段阻尼外推(γ)到目标，解码时潜空间位移只加到原细胞非零HVG位置。

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

来自 programs.git 3029342b48:solution/METHOD.md

调研员的计划

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

对比：这个提交的上一版（种子程序：相对空仓库）。改动的文件：solution/EXECUTION.json +1 −0、solution/METHOD.md +45 −0、solution/run.py +272 −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..5b16bbc--- /dev/null+++ b/solution/METHOD.md@@ -0,0 +1,45 @@+manifold_ode：25维PCA潜空间里训练小型非自治MLP速度场v(z,t)（Sinkhorn分布匹配+kNN流形约束），从最后输入阶段阻尼外推(γ)到目标，解码时潜空间位移只加到原细胞非零HVG位置。++## 方法（family_id = manifold_ode，按 PLAN 实现）++- 表示：两输入阶段合并 → top-2500 HVG → PCA(25) → 每维按全局std归一化得 z。+- 速度场：MLP [64,64] tanh，输入 (z,t)，末层权重×0.01 初始化（初始≈恒等）。非自治（t 显式输入）。+- 训练：300 步 Adam(1e-3)，batch 512；RK4 4步从 t=0 积分到 t=1；损失 = debiased Sinkhorn(eps=0.05², 30 iters, 纯 torch 实现——本环境 geomloss 的 multiscale/tensor backend 均损坏) + λ_kin·mean‖v‖² (λ_kin=0.1) + λ_man·0.01·mean‖z_pred − mean(3NN in E9.5 latent)‖² (λ_man=0.5，kNN 流形约束，MIOFlow 启发的机制)。+- 外推：最后阶段细胞从 t=1 积分到 t_end = 1 + (t_target − t_last)/(t_last − t_first)（只用时间差，视图无关；X3 上 t_end=3），位移乘阻尼 γ=0.5。+- 解码：x_final = x_last + clip 到非零掩码的 (V^T·Δz)，clip≥0；非HVG列原样复制最后阶段。输出细胞 = 从最后阶段无放回抽样 max_cells 个（rng(seed)）。+- 单输入阶段退路：copy_last。+- 确定性：np.random.default_rng(seed) + torch.manual_seed(seed)，torch.set_num_threads(1)（多线程下 logsumexp 有严重线程争用，单线程反而快 25 倍）。CPU 运行，EXECUTION.json gpu=false，全程 <2 min。++## 关闭机制对照（mechanism_off_control）++`MECH_OFF=1` 环境变量：v(z,t) 换成常向量 (mean z_last − mean z_first)，λ_man=0，跳过训练，其余管线（积分、γ、解码）完全相同。提交默认打开机制（MECH_OFF 未设）。++X3 A半查分（seed 0）：+| 配置 | board | cell_state | covariation | de_recovery | direction |+|---|---|---|---|---|---|+| ODE γ=0.5, λ_man=0.5（提交默认） | **48.10** | 45.82 | 49.35 | 48.62 | 49.33 |+| ODE γ=0.3 | 48.15 | 47.36 | 49.00 | 47.32 | 49.23 |+| ODE γ=0.7 | 47.71 | 44.21 | 49.66 | 48.62 | 49.45 |+| ODE γ=0.2 | 47.83 | 48.11 | 48.82 | 45.30 | 49.22 |+| off-control 常数位移 γ=0.5 | 47.32 | 48.34 | 48.08 | 43.80 | 49.00 |+| copy_last 参考 | 47.62 | 49.51 | 48.59 | 43.09 | 49.12 |++## 机制生效证据（默认配置，seed 0）++- corr(逐细胞ODE位移, 常位移) = 0.67（<0.99，速度场确为位置依赖）；off-control 下 = 1.0000。+- 位移范数的细胞间 std = 0.557（常数位移为 ~2e-7）。+- frac_closer_manifold（ODE 比常数位移更靠近 E9.5 流形的细胞比例）= 0.47。+- 预测云均值 ODE vs off-control W2 差 = 0.40（>0，两者分布不同）。+- off-control 本身也验证过：其证据量退化为常数（std≈0, corr≈1, W2≈0），确认开关有效。+- 分数上：ODE(48.10) > off-control(47.32) 与 copy_last(47.62)，主要增益在 de_recovery（48.6 vs 43.8/43.1）和 covariation（49.35 vs 48.1/48.6）；差距 ~0.5–0.8 分，**小于 T1 ~2 分的噪声水平**，不能断言显著，只作方向性证据。++## 已验证 / 未验证++- 已验证：X3 视图完整跑通 + vec-check ok；γ∈{0.2,0.3,0.5,0.7} 扫描（A半）；off-control 对照；单线程确定性（同 seed 重跑逐位一致，源于固定 rng/无全局随机依赖）。+- 未验证：λ_kin、λ_man 的扫描（预算内只完成 γ 轴）；final 视图（E8.5+E9.5→E10.5，双官方阶段）上的表现——本节点评分只含 X3；伪装视图重跑（程序只依赖数据与时间差，无路径/绝对时间分支，预期通过）。+- 生物学知识使用：无来自保留阶段/保留基因型的信息；方法为纯数据驱动的动力学外推。未读 external/（X3 视图挂的 external 恰是输入阶段本身，无额外信息）。未读 prior/。+- 弱点：cell_state 略低于 copy_last（外推使细胞云轻微偏离目标边缘分布）；direction 四组分 ~49，与全表其它节点一样卡在地板附近——单区间观测下方向信息本质上弱可辨。++## 下一步建议++1) λ_man 扫描（0/0.1/0.5/1.0）与 γ 交互；2) 解码时对预测细胞做按细胞型的组成保持重采样；3) 双输入时按型内局部 OT 耦合训练场（当前用全局 Sinkhorn，型间混合可能污染速度场）。diff --git a/solution/run.py b/solution/run.pynew file mode 100644index 0000000..9fe9288--- /dev/null+++ b/solution/run.py@@ -0,0 +1,272 @@+"""manifold_ode draft: PCA-latent non-autonomous neural ODE with kNN manifold+constraint, Sinkhorn distribution-matching loss, damped extrapolation.++Family: manifold_ode (PLAN node 18).+Off-control: env MECH_OFF=1 replaces v(z,t) by the constant mean latent+displacement and disables the manifold penalty.+"""+from __future__ import annotations++import argparse+import json+import os+from pathlib import Path++import numpy as np+import scipy.sparse as sp+++# ---------------- config (PLAN defaults) ----------------+N_HVG = 2500+N_PCA = 25+HIDDEN = (64, 64)+STEPS = 300+BATCH = 512+LR = 1e-3+LAM_KIN = float(os.environ.get("LAM_KIN", "0.1"))+LAM_MAN = float(os.environ.get("LAM_MAN", "0.5"))+GAMMA = float(os.environ.get("GAMMA", "0.5"))+KNN_MAN = 3+MECH_OFF = os.environ.get("MECH_OFF", "0") == "1"+BLUR = 0.05+++def read_h5ad(path):+    import anndata+    return anndata.read_h5ad(path)+++def main():+    ap = argparse.ArgumentParser()+    ap.add_argument("--data", required=True)+    ap.add_argument("--out", required=True)+    ap.add_argument("--seed", type=int, default=0)+    args = ap.parse_args()+    view = Path(args.data)+    man = json.loads((view / "manifest.json").read_text())+    genes = [l.strip() for l in (view / man["genes_file"]).read_text().splitlines() if l.strip()]+    n_genes = len(genes)+    inputs = sorted(man["inputs"], key=lambda e: e["time"])+    t_target = float(man["target"]["time"])++    seed = args.seed+    np.random.seed(seed)+    rng = np.random.default_rng(seed)++    last = read_h5ad(view / inputs[-1]["path"])+    Xl = last.X.astype(np.float32)+    Xl = Xl.tocsr()+    n_last = Xl.shape[0]++    n_out = int(man["max_cells"])+    n_out = min(max(n_out, int(man["min_cells"])), max(n_last, int(man["min_cells"])))+    if n_last >= n_out:+        sel = rng.choice(n_last, size=n_out, replace=False)+        sel.sort()+    else:+        sel = rng.choice(n_last, size=n_out, replace=True)++    out = sp.csr_matrix(Xl[sel].copy())++    if len(inputs) < 2:+        # single-stage fallback: copy_last+        write_out(out, genes, args.out, man)+        return++    first = read_h5ad(view / inputs[0]["path"])+    Xf = first.X.astype(np.float32).tocsr()++    dt_in = float(inputs[-1]["time"]) - float(inputs[0]["time"])+    dt_ext = t_target - float(inputs[-1]["time"])+    ratio = dt_ext / dt_in if dt_in > 0 else 1.0+    t_end = 1.0 + max(ratio, 0.0)++    # -------- HVG selection on combined stages --------+    Xc = sp.vstack([Xf, Xl]).tocsr()+    n = Xc.shape[0]+    s1 = np.asarray(Xc.sum(axis=0)).ravel() / n+    s2 = np.asarray(Xc.multiply(Xc).sum(axis=0)).ravel() / n+    var = s2 - s1 ** 2+    hvg_idx = np.argsort(-var)[:N_HVG]+    hvg_idx = np.sort(hvg_idx)++    F = Xc[:, hvg_idx].toarray().astype(np.float32)+    mean_h = F.mean(axis=0)+    Fc = F - mean_h+    # -------- PCA --------+    from sklearn.decomposition import PCA+    pca = PCA(n_components=N_PCA, random_state=seed)+    Z = pca.fit_transform(Fc).astype(np.float32)+    V = pca.components_.astype(np.float32)  # (N_PCA, N_HVG)+    scale = Z.std(axis=0)+    scale[scale < 1e-8] = 1.0+    Zn = Z / scale  # normalize dims for stable sinkhorn blur+    nf = Xf.shape[0]+    z0 = Zn[:nf]+    z1 = Zn[nf:]++    import torch+    torch.set_num_threads(1)+    torch.manual_seed(seed)+    dev = "cpu"+    tz0 = torch.tensor(z0, device=dev)+    tz1 = torch.tensor(z1, device=dev)+    d = N_PCA++    class Vec(torch.nn.Module):+        def __init__(self):+            super().__init__()+            layers = []+            prev = d + 1+            for h in HIDDEN:+                layers += [torch.nn.Linear(prev, h), torch.nn.Tanh()]+                prev = h+            layers += [torch.nn.Linear(prev, d)]+            self.net = torch.nn.Sequential(*layers)+            with torch.no_grad():+                self.net[-1].weight.mul_(0.01)+                self.net[-1].bias.zero_()++        def forward(self, z, t):+            tt = torch.full((z.shape[0], 1), float(t), device=z.device)+            return self.net(torch.cat([z, tt], dim=1))++    const_v = (z1.mean(axis=0) - z0.mean(axis=0))+    const_v = torch.tensor(const_v, device=dev)++    def integrate(z, t0, t1, model, n_steps=4):+        h = (t1 - t0) / n_steps+        t = t0+        for _ in range(n_steps):+            k1 = vel(model, z, t)+            k2 = vel(model, z + 0.5 * h * k1, t + 0.5 * h)+            k3 = vel(model, z + 0.5 * h * k2, t + 0.5 * h)+            k4 = vel(model, z + h * k3, t + h)+            z = z + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)+            t = t + h+        return z++    def vel(model, z, t):+        if MECH_OFF:+            return const_v.expand_as(z)+        return model(z, t)++    # -------- manifold target: last-stage latent cloud --------+    from scipy.spatial import cKDTree+    tree1 = cKDTree(z1)++    model = Vec().to(dev)+    params = [p for p in model.parameters() if p.requires_grad]+    opt = torch.optim.Adam(params, lr=LR) if params else None++    def sinkhorn_cost(a, b, eps=None, n_iter=30):+        # entropic OT cost (debiased), squared-euclidean ground cost+        if eps is None:+            eps = BLUR ** 2+        Caa = torch.cdist(a, a) ** 2+        Cbb = torch.cdist(b, b) ** 2+        Cab = torch.cdist(a, b) ** 2++        def sinkhorn_dual(C):+            f = torch.zeros(C.shape[0], device=C.device)+            g = torch.zeros(C.shape[1], device=C.device)+            for _ in range(n_iter):+                f = -eps * torch.logsumexp((g[None, :] - C) / eps, dim=1)+                g = -eps * torch.logsumexp((f[:, None] - C) / eps, dim=0)+            return f, g++        fab, gab = sinkhorn_dual(Cab)+        faa, gaa = sinkhorn_dual(Caa)+        fbb, gbb = sinkhorn_dual(Cbb)+        Pab = torch.exp((fab[:, None] + gab[None, :] - Cab) / eps)+        Paa = torch.exp((faa[:, None] + gaa[None, :] - Caa) / eps)+        Pbb = torch.exp((fbb[:, None] + gbb[None, :] - Cbb) / eps)+        return (Pab * Cab).sum() - 0.5 * (Paa * Caa).sum() - 0.5 * (Pbb * Cbb).sum()++    def sink(a, b):+        return sinkhorn_cost(a, b)++    n0, n1 = z0.shape[0], z1.shape[0]+    bs = min(BATCH, n0, n1)+    for step in range(0 if MECH_OFF else STEPS):+        opt.zero_grad()+        i0 = torch.randint(0, n0, (bs,))+        i1 = torch.randint(0, n1, (bs,))+        zb = tz0[i0]+        # small time jitter for non-autonomous coverage+        traj_z = integrate(zb, 0.0, 1.0, model, n_steps=4)+        l_main = sink(traj_z, tz1[i1])+        v_mid = vel(model, traj_z.detach() * 0 + zb.detach(), 0.5)+        l_kin = (v_mid ** 2).sum(dim=1).mean()+        loss = l_main + LAM_KIN * l_kin+        if (not MECH_OFF) and LAM_MAN > 0:+            zp = traj_z.detach().cpu().numpy()+            _, idx = tree1.query(zp, k=KNN_MAN)+            tgt = torch.tensor(z1[idx].mean(axis=1), device=dev, dtype=torch.float32)+            l_man = ((traj_z - tgt) ** 2).sum(dim=1).mean()+            loss = loss + LAM_MAN * 0.01 * l_man+        loss.backward()+        opt.step()++    # -------- extrapolate last-stage cells to target --------+    zlast = torch.tensor(z1, device=dev)+    with torch.no_grad():+        zpred_full = integrate(zlast, 1.0, t_end, model, n_steps=4)+        disp = zpred_full - zlast+        zpred = zlast + GAMMA * disp++        # off-control comparison / mechanism evidence+        zc = zlast + GAMMA * const_v * (t_end - 1.0)+        ev = {}+        dn = disp.cpu().numpy()+        cn = (const_v * (t_end - 1.0)).cpu().numpy()+        norms = np.linalg.norm(dn, axis=1)+        ev["disp_norm_std"] = float(norms.std())+        if norms.std() > 1e-8:+            cc = np.corrcoef(dn.ravel(), np.broadcast_to(cn, dn.shape).ravel())[0, 1]+            ev["corr_vs_const"] = float(cc)+        else:+            ev["corr_vs_const"] = 1.0+        zp_np = zpred.cpu().numpy()+        zc_np = zc.cpu().numpy()+        d_ode, _ = tree1.query(zp_np, k=1)+        d_con, _ = tree1.query(zc_np, k=1)+        ev["frac_closer_manifold"] = float((d_ode < d_con).mean())+        from scipy.spatial.distance import cdist+        ev["w2_ode_vs_control"] = float(np.sqrt(((zp_np.mean(0) - zc_np.mean(0)) ** 2).sum()))+        print("MECHANISM_EVIDENCE", json.dumps(ev))++    # -------- decode --------+    Zp = (zpred.cpu().numpy() * scale).astype(np.float32)+    dH = Zp @ V  # (n1, N_HVG) latent-space delta vs reconstruction of z1+    Z1u = (z1 * scale).astype(np.float32)+    dH = dH - Z1u @ V++    Xl_sel = Xl[sel].tocsr()+    Hd = Xl_sel[:, hvg_idx].toarray()+    mask = Hd > 0+    Xnew = np.clip(Hd + np.where(mask, dH[sel], 0.0), 0, None).astype(np.float32)+    rows, cols = np.nonzero(Xnew)+    A = sp.coo_matrix((Xnew[rows, cols], (rows, cols)), shape=mask.shape).tocsr()+    hvg_mask = np.zeros(n_genes, dtype=bool)+    hvg_mask[hvg_idx] = True+    B = Xl_sel[:, ~hvg_mask].tocsr()+    out = sp.hstack([A, B]).tocsr()+    order = np.concatenate([hvg_idx, np.nonzero(~hvg_mask)[0]])+    out = out[:, np.argsort(order)].tocsr()+    out.eliminate_zeros()+    write_out(out, genes, args.out, man)+++def write_out(X, genes, path, man):+    import anndata+    X = sp.csr_matrix(X)+    adata = anndata.AnnData(X=X.astype(np.float32))+    adata.obs_names = [f"cell{i}" for i in range(X.shape[0])]+    adata.var_names = genes+    Path(path).parent.mkdir(parents=True, exist_ok=True)+    adata.write_h5ad(path)+++if __name__ == "__main__":+    main()

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

用到的知识库条目

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

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

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

分析员

ANALYSIS.json written

工程师

完成。总结：

**做了什么**：从零实现 manifold_ode 家族 draft（node 18）：25维 PCA 潜空间的非自治 MLP 速度场 v(z,t)，自写 debiased Sinkhorn（geomloss 后端在本环境损坏）+ kNN 流形约束（λ_man=0.5）训练 300 步，RK4 积分，γ=0.5 阻尼外推，位移只加到最后阶段细胞的非零 HVG 位置；含 `MECH_OFF=1` 常数位移对照、单输入 copy_last 退路、CPU 确定性运行（<2min）。

**查分（X3 A半，7 次）**：ODE γ=0.5 **48.10** / γ=0.3 48.15 / γ=0.2 47.83 / γ=0.7 47.71；off-control 47.32；copy_last 47.62。机制证据：逐细胞位移 vs 常位移相关 0.67（对照为 1.0）、位移范数 std 0.56、增益主要在 de_recovery（48.6 vs 43.1）——但总差距 <2 分噪声，仅方向性成立。

**下一步最值得试**：λ_man/λ_kin 扫描与 γ 交互；按细胞型局部分布匹配代替全局 Sinkhorn（避免型间混合污染速度场）；cell_state 是当前最弱组分，可试解码端组成保持重采样。

调研员

are not conclusive -- use multiple queries if borderline.", "family_id": "manifold_ode", "mechanism": "A small non-autonomous MLP vector field v(z,t) in 25-dim PCA space, trained with Sinkhorn distribution-matching loss between integrated E8.5 and observed E9.5, plus a kNN manifold penalty that pulls predictions toward the observed E9.5 data manifold; extrapolation from E9.5 to target uses damped integration (gamma<1).", "vs_constant_shift": "Constant shift applies one global displacement vector to all cells regardless of state. This ODE's velocity depends on both position z and time t, so different regions of latent space move differently; the manifold constraint further bends trajectories toward observed data density rather than following straight lines. The off-control (constant mean displacement) directly tests whether the position-dependent velocity + manifold pull produce different predictions.", "mechanism_evidence": "Engineer should report: (1) pairwise correlation between per-cell displacement vectors from ODE vs constant shift -- if >0.99, mechanism is effectively constant; (2) std of displacement norms across cells -- constant shift gives zero variance, ODE should give nonzero; (3) fraction of predicted cells whose nearest E9.5 neighbor distance is smaller under ODE than under constant shift (manifold pull effect); (4) four-group score breakdown for ODE vs constant-shift control vs copy_last.", "mechanism_off_control": "Same program, same data pipeline, same decoding; replace v(z,t) with a constant vector (mean of E9.5 latent minus mean of E8.5 latent) and set lambda_manifold=0. Integrate with same solver and gamma. Expected difference: if the manifold ODE mechanism is active, the ODE run should produce a different spatial distribution in latent space (measurable by Wasserstein distance between the two predicted clouds > 0) and different four-group scores, particularly in direction and covariation. If outputs are identical, the mechanism is not running."}