总览 · ← 返回运行 20261003-093415-search-t1-r2-D-s1

节点 n21

lowrank_shape 最小实现：Ledoit-Wolf 型内协方差低秩外推仿射变形；X3 上开/关对照显示机制单调有害（cell_state、covariation 齐降），故默认 α=0（关闭），负结果如实上报。

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

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

调研员的计划

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

对比：这个提交的上一版（种子程序：相对空仓库）。改动的文件：solution/EXECUTION.json +1 −0、solution/METHOD.md +55 −0、solution/run.py +257 −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..8cf4510--- /dev/null+++ b/solution/METHOD.md@@ -0,0 +1,55 @@+lowrank_shape 最小实现：Ledoit-Wolf 型内协方差低秩外推仿射变形；X3 上开/关对照显示机制单调有害（cell_state、covariation 齐降），故默认 α=0（关闭），负结果如实上报。++# lowrank_shape（节点 21，draft）++## 方法（完全按 PLAN family=lowrank_shape 实现）++1. 读全部输入阶段（`view_io.inputs_by_time`，官方面板顺序）。输出细胞 = 最后输入阶段按+   celltype 分层比例抽样到 `target_n_cells`（保组成，不改比例）。+2. ≥2 个输入阶段且 α>0 时启用机制：+   - 两阶段合并选方差 top 2500 HVG（逐块稠密计算基因矩，不整表稠密）；按合并均值/方差 z-score；+   - randomized SVD（`random_state=seed`，确定性）取 30 个 PC 载荷 V（拟合每阶段最多抽 12000 细胞）；+   - 有效类型 = 两阶段都有且各 ≥30 细胞的 celltype 标签（含 "Unknown"；无标签场景未出现，未启用 KMeans 分支）；+   - 每型：`sklearn.covariance.LedoitWolf` 估 30 维 PCA 坐标协方差 Σ0、Σ1；ΔΣ=Σ1−Σ0 特征分解，+     取 |λ| top-10 重建 ΔΣ_lr；Σ_pred = Σ1 + α·ΔΣ_lr，特征值裁剪 ≥1e-4；+     A = Σ_pred^{1/2} Σ1^{−1/2}（eigh）；+   - 输出中该型细胞：z_new = mean1 + A(z−mean1)，x_new = (V^T z_new + r)·sd + mu（r=原残差），夹 ≥0；+   - 非 HVG 基因、无效类型细胞：原样复制最后输入阶段。+3. 单输入阶段（proxy 视图）：无 ΔΣ 可估，机制关闭，输出 = 末阶段分层子样（PLAN 步骤 7 的退路）。++参数：n_hvg=2500, n_pcs=30, k=10, min_cells=30, eigen_clip=1e-4，均按 PLAN 初值；α 见下。+知识来源：无生物学先验被写入程序（纯统计方法，不用通路/谱系/文献知识；prior/、external/ 未读取）。+视图无关：不读视图路径/文件名/board 字段，不用绝对时间（α 为常数），伪装视图下输出逐元素相同。++## 机制生效证据（α=0.3 在 X3 上，stderr JSON）++- 有效类型 4 个：IFT-CM(133/489)、OFT/RV-CM(327/391)、SV-CM(238/176)、Unknown(377/437)；+  实际被修改细胞 448/652（69%，即全部属于有效类型的输出细胞）。+- ΔΣ 前 10 特征值能量占比 0.95–0.99（低秩近似充分，但也提示 ΔΣ 可能被少数方向/噪声主导）。+- 型内协方差 Frobenius 比 1.33–1.96（形状确实被改变，机制在运行）。+- 基因空间最大均值漂移 0.05–1.55 log 单位（PCA 空间均值严格不动；漂移来自夹 ≥0，SV-CM 最大）。+- 四组分变化（X3，A 半，同一 subsample 流程）：+  | 配置 | de_recovery | direction | cell_state | covariation | 总分 |+  |---|---|---|---|---|---|+  | 关（α=0，=copy 末阶段子样） | 44.54 | 49.87 | 51.06 | 49.74 | 48.87 |+  | 开 α=0.1 | 44.17 | 49.87 | 38.09 | 43.38 | 43.61 |+  | 开 α=0.3 | 43.80 | 49.49 | 36.37 | 39.25 | 42.08 |+  | 开 α=0.3 + 稀疏解码（仅改原非零位） | 42.74 | 50.50 | 21.67 | 27.73 | 35.36 |+  proxy10（单输入，机制恒关）：53.24。++## 结论（负结果）++机制按计划实现并确实在运行（上表证据），但方向性预测失败：正 α 外推使型内分布膨胀+（mmd_u 0.033→0.052，variogram 0.00156→0.00227），恰是 PLAN 预期改善的 cell_state、+covariation 两组单调恶化；DE 两组基本不动（均值保持符合设计）；α 越小损失越小但始终为负，+且损害在 α=0.1 已大部分出现（非线性），说明主害不是解码稠密化而是协方差膨胀本身——+两阶段型内协方差差 ΔΣ 在 X3 的细胞数（每型百级）下更可能是采样噪声/技术差异而非发育趋势+（PLAN risks #1 命中）。稀疏解码变体更差，排除“稠密化破坏共变”的解释。+据此默认 α=0（同节点 18 的处理）：提交程序 = 分层子样 copy_last。+复现机制：`LOWRANK_ALPHA=0.3 python run.py ...`；对照：`VEC_MECH_OFF=1`。++## 验证过 / 未验证++- 验证：proxy、X3 两视图跑通 + vec-check ok；seed 0 复跑逐元素一致（determinism）；X3 开/关对照 4 组配置各 1 次查分（共 5 次查询）。+- 未验证：α<0（收缩方向）——不在 PLAN 搜索范围且等于按评分反向调参，未做；k、n_pcs 未扫（α 单调有害，扫其余超参无意义）；final 视图（E8.5+E9.5 双输入，机制会启用）未跑——本节点无权访问，代码路径与 X3 相同（每阶段拟合抽样上限 12000 细胞控制内存/时间）。+- 运行时间：X3 ~40 s、proxy ~10 s，纯 CPU（EXECUTION.json gpu=false），远低于 limits。diff --git a/solution/run.py b/solution/run.pynew file mode 100644index 0000000..54d81e8--- /dev/null+++ b/solution/run.py@@ -0,0 +1,257 @@+"""lowrank_shape: Ledoit-Wolf shrinkage covariance low-rank extrapolation.++Family: lowrank_shape (PLAN). Mechanism: within each valid cell type, estimate+stage-0/stage-1 covariance of the 30-d PCA coordinates (Ledoit-Wolf shrinkage),+take the top-k low-rank part of dSigma = Sigma1 - Sigma0, extrapolate+Sigma_pred = Sigma1 + alpha * dSigma_lr, and map stage-1 cells by the affine+shape transform A = Sigma_pred^(1/2) Sigma1^(-1/2) centered at the type mean.+Means are unchanged; only the within-type distribution shape changes.+Single input stage: no dSigma available -> mechanism off, output = copy of the+last input stage (subsampled).++Control switch: env VEC_MECH_OFF=1 forces alpha=0 (identity transform).+"""++from __future__ import annotations++import argparse+import json+import os+import sys+from pathlib import Path++import numpy as np+from scipy import sparse++sys.path.insert(0, str(Path(__file__).resolve().parent))++from src.task1_temporal import view_io  # noqa: E402  (harness puts modeling/ on PYTHONPATH)++N_HVG = 2500+N_PCS = 30+# PLAN initial alpha = 0.3. On/off control on X3 (weight 2): alpha=0.3 -> 42.08,+# alpha=0.1 -> 43.61, off (alpha=0) -> 48.87; cell_state and covariation (the two+# groups the mechanism targets) both got monotonically worse with alpha.+# Negative result: default 0.0 (same policy as node 18). Re-enable for tests with+# env LOWRANK_ALPHA=0.3; env VEC_MECH_OFF=1 forces the identity control.+ALPHA = 0.0+K_LOW_RANK = 10+MIN_CELLS_PER_TYPE = 30+EIGEN_CLIP = 1e-4+FIT_CAP_PER_STAGE = 12000+++def log(obj):+    print(json.dumps(obj), file=sys.stderr, flush=True)+++def gene_moments(adatas):+    """Per-gene mean and E[x^2] over the concatenation of stages (chunked)."""+    n_genes = adatas[0].shape[1]+    s = np.zeros(n_genes, dtype=np.float64)+    s2 = np.zeros(n_genes, dtype=np.float64)+    n = 0+    for a in adatas:+        X = (a.X if hasattr(a, "X") else a).tocsr()+        for i0 in range(0, X.shape[0], 2048):+            B = X[i0:i0 + 2048]+            Bd = np.asarray(B.todense(), dtype=np.float64)+            s += Bd.sum(axis=0)+            s2 += (Bd * Bd).sum(axis=0)+            n += Bd.shape[0]+    mean = s / n+    var = np.maximum(s2 / n - mean * mean, 0.0)+    return mean, var+++def stratified_subsample(labels, n_out, rng):+    """Proportional stratified sample of row indices, deterministic given rng."""+    n = len(labels)+    if n_out >= n:+        return np.arange(n)+    idx = []+    types = np.unique(labels)+    quotas = {}+    acc = 0.0+    for t in types:+        c = int((labels == t).sum())+        q = int(np.floor(c * n_out / n))+        acc += c * n_out / n - q+        quotas[t] = q+    # distribute leftover by largest fractional remainder, deterministic order+    rem = n_out - sum(quotas.values())+    frac = []+    for t in types:+        c = int((labels == t).sum())+        frac.append((-(c * n_out / n - np.floor(c * n_out / n)), t))+    frac.sort()+    for _, t in frac[:rem]:+        quotas[t] += 1+    for t in types:+        rows = np.flatnonzero(labels == t)+        q = min(quotas[t], len(rows))+        if q > 0:+            idx.append(np.sort(rng.choice(rows, size=q, replace=False)))+    return np.sort(np.concatenate(idx)) if idx else np.arange(0)+++def mat_sqrt_parts(S, floor):+    """eigh-based: return (U, sqrt(eig clipped at floor)) for S PSD-ish."""+    w, U = np.linalg.eigh((S + S.T) / 2.0)+    w = np.maximum(w, floor)+    return U, np.sqrt(w)+++def lowrank_extrapolate(S0, S1, alpha, k):+    dS = S1 - S0+    w, U = np.linalg.eigh((dS + dS.T) / 2.0)+    order = np.argsort(-np.abs(w))+    keep = order[:k]+    dS_lr = (U[:, keep] * w[keep]) @ U[:, keep].T+    S_pred = S1 + alpha * dS_lr+    energy = float(w[keep] ** 2 @ np.ones(len(keep)) / max((w ** 2).sum(), 1e-12))+    return S_pred, energy+++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)+    args = ap.parse_args()++    view = args.data+    rng = np.random.default_rng(args.seed)+    manifest = view_io.load_manifest(view)+    genes = view_io.panel_genes(view, manifest)+    n_panel = len(genes)++    mech_off = os.environ.get("VEC_MECH_OFF", "0") == "1"+    alpha = 0.0 if mech_off else float(os.environ.get("LOWRANK_ALPHA", ALPHA))++    inputs = view_io.inputs_by_time(manifest)+    a_last = view_io.read_stage(view, inputs[-1], genes, missing="fill")+    labels_last = view_io.labels_of(a_last)+    n_out = view_io.target_n_cells(manifest, a_last.n_obs, cap=manifest["max_cells"])+    out_rows = stratified_subsample(labels_last, n_out, rng)++    two_stage = len(inputs) >= 2 and alpha > 0.0+    X_out = a_last.X.tocsr()[out_rows]++    if two_stage:+        a_prev = view_io.read_stage(view, inputs[-2], genes, missing="fill")+        labels_prev = view_io.labels_of(a_prev)+        adatas = [a_prev, a_last]++        # PCA fit on subsamples of both stages+        fit_idx = []+        for a in adatas:+            r = np.random.default_rng(args.seed + 7)+            if a.n_obs > FIT_CAP_PER_STAGE:+                fit_idx.append(r.choice(a.n_obs, size=FIT_CAP_PER_STAGE, replace=False))+            else:+                fit_idx.append(np.arange(a.n_obs))++        gmean, gvar = gene_moments([a.X[fit_idx[i]] for i, a in enumerate(adatas)])+        hvg = np.argsort(-gvar)[:N_HVG]+        hvg = hvg[gvar[hvg] > 1e-12]+        mu = gmean[hvg]+        sd = np.sqrt(gvar[hvg])+        sd = np.maximum(sd, 1e-8)++        # z-scored stacked fit matrix+        Zfit = []+        for i, a in enumerate(adatas):+            B = np.asarray(a.X[fit_idx[i]][:, hvg].todense(), dtype=np.float32)+            Zfit.append((B - mu) / sd)+        Zfit = np.vstack(Zfit).astype(np.float64)+        Zfit -= Zfit.mean(axis=0, keepdims=True)+        # randomized SVD -> loadings V (n_hvg x n_pcs), deterministic+        from sklearn.utils.extmath import randomized_svd+        n_pcs = int(min(N_PCS, Zfit.shape[0] - 1, Zfit.shape[1] - 1))+        _, _, Vt = randomized_svd(Zfit, n_components=n_pcs, n_oversamples=20,+                                  random_state=args.seed)+        V = Vt.T.copy()+        del Zfit, Vt++        def coords(X_hvg_dense):+            Zc = (X_hvg_dense - mu) / sd+            return Zc @ V, Zc++        # valid types: >= MIN_CELLS in both stages+        types = sorted(set(labels_prev) & set(labels_last))+        valid = [t for t in types+                 if int((labels_prev == t).sum()) >= MIN_CELLS_PER_TYPE+                 and int((labels_last == t).sum()) >= MIN_CELLS_PER_TYPE]++        from sklearn.covariance import LedoitWolf++        Z1_all = np.asarray(a_last.X[:, hvg].todense(), dtype=np.float32).astype(np.float64)+        Z0_all = np.asarray(a_prev.X[:, hvg].todense(), dtype=np.float32).astype(np.float64)+        z1_all, Zc1_all = coords(Z1_all)+        z0_all, _ = coords(Z0_all)+        del Z0_all, Z1_all++        new_hvg = np.asarray(X_out[:, hvg].todense(), dtype=np.float32).astype(np.float64)+        out_labels = labels_last[out_rows]+        stats = {"types": {}, "n_cells_changed": 0}++        for t in valid:+            m0 = labels_prev == t+            m1 = labels_last == t+            S0 = LedoitWolf().fit(z0_all[m0]).covariance_+            S1 = LedoitWolf().fit(z1_all_t := z1_all[m1]).covariance_+            mean1 = z1_all_t.mean(axis=0)+            S_pred, energy = lowrank_extrapolate(S0, S1, alpha, K_LOW_RANK)+            Up, sp = mat_sqrt_parts(S_pred, EIGEN_CLIP)+            U1, s1 = mat_sqrt_parts(S1, 1e-6)+            A = (Up * sp) @ Up.T @ (U1 * (1.0 / s1)[None, :]) @ U1.T++            rows_out = np.flatnonzero(out_labels == t)+            if rows_out.size == 0:+                continue+            z = z1_all[rows_out]+            r = ((Zc1_all[rows_out]) - z @ V.T)+            z_new = mean1 + (z - mean1) @ A.T+            Zc_new = z_new @ V.T + r+            recon = np.maximum(Zc_new * sd + mu, 0.0)+            if os.environ.get("DECODE_SPARSE", "0") == "1":+                orig = new_hvg[rows_out]+                recon = np.where(orig > 0, recon, 0.0)+            new_hvg[rows_out] = recon+            stats["n_cells_changed"] += int(rows_out.size)++            # evidence: mean shift in gene space (should be tiny), cov Frobenius ratio+            x_old = np.asarray(a_last.X[rows_out][:, hvg].todense(), dtype=np.float64)+            x_new = new_hvg[rows_out]+            mean_shift = float(np.abs(x_new.mean(0) - x_old.mean(0)).max())+            zt = z - mean1+            cov_emp_after = (zt @ A.T).T @ (zt @ A.T) / max(zt.shape[0] - 1, 1)+            frob_ratio = float(np.linalg.norm(cov_emp_after) / max(np.linalg.norm(S1), 1e-12))+            stats["types"][t] = {+                "n0": int(m0.sum()), "n1": int(m1.sum()),+                "dS_topk_energy": round(energy, 3),+                "max_abs_mean_shift_genespace": round(mean_shift, 4),+                "cov_frob_ratio_vs_S1": round(frob_ratio, 3),+            }+        log({"mechanism": "lowrank_shape", "alpha": alpha, "k": K_LOW_RANK,+             "n_valid_types": len(valid), "stats": stats})++        # rebuild output matrix with modified HVG columns+        hvg_block = sparse.csr_matrix(np.asarray(new_hvg, dtype=np.float32))+        non_hvg_mask = np.ones(n_panel, dtype=bool)+        non_hvg_mask[hvg] = False+        non_hvg_idx = np.flatnonzero(non_hvg_mask)+        combined = sparse.hstack([X_out[:, non_hvg_idx].tocsc(), hvg_block.tocsc()]).tocsc()+        cols = np.concatenate([non_hvg_idx, hvg])+        inv = np.argsort(cols)+        X_out = combined[:, inv].tocsr()+    else:+        log({"mechanism": "lowrank_shape", "alpha": 0.0, "note": "single input or mech off: copy last stage"})++    view_io.write_prediction(X_out.tocsr(), genes, args.out, seed=args.seed)+    log({"n_cells": int(X_out.shape[0]), "nnz": int(X_out.nnz), "out": args.out})+++if __name__ == "__main__":+    main()

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

用到的知识库条目

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

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

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

分析员

ANALYSIS.json written

工程师

完成。节点 21（draft，family=lowrank_shape）总结：

**做了什么**：按 PLAN 完整实现 Ledoit-Wolf 型内协方差低秩外推仿射变形（HVG2500→PCA30→每型 Σ0/Σ1 收缩估计→ΔΣ top-10 低秩→A=Σ_pred^½Σ₁^−½，均值不动），含单输入退路和 `VEC_MECH_OFF`/`LOWRANK_ALPHA` 对照开关；proxy、X3 均跑通 vec-check，seed 确定性验证通过，纯 CPU。

**查分（5/20 次）**：X3 关(α=0)=48.87，开 α=0.1=43.61，α=0.3=42.08，α=0.3+稀疏解码=35.36；proxy10=53.24。机制证据齐全（448/652 细胞被改、Frobenius 比 1.33–1.96、ΔΣ top-10 能量 0.95+），但正 α 使 cell_state、covariation 单调恶化——PLAN risks#1 命中：百级细胞数下 ΔΣ 更像采样噪声而非发育趋势。已如实记录负结果，默认 α=0（同节点 18 先例），预估节点分 ≈ (53.2+2×48.9)/3 ≈ 50.3。

**下一步最值得试**：本方向可放弃正 α 协方差膨胀；树里更有效的仍是组成/选择类杠杆（节点 7/16 的 composition_trend），X3 上 copy 类基线 cell_state≈51 已近地板，想超过需在型内成熟度选择或跨数据集（Qiu）标签对齐上做文章，而非表达解码。

调研员

3 为初始步长。对 Σ_pred 做特征值裁剪：任何 < 1e-4 的特征值替换为 1e-4，确保正定。\n5. 仿射变换矩阵：A = Σ_pred^{1/2} @ Σ₁^{−1/2}（均为 30×30，通过特征分解计算矩阵平方根和逆平方根）。对 E9.5 该型的每个细胞：z_new = z_mean₁ + A @ (z − z_mean₁)，其中 z_mean₁ 是该型在阶段 1 的 PCA 均值。基因空间输出：x_new_zscore = V^T z_new + r（加回原始残差），再反标准化回原始表达空间。\n6. 非 HVG 基因和无效类型的细胞：直接复制最后一个输入阶段的值。\n7. 单输入退路（proxy10）：只有一个输入阶段时，无法估计 ΔΣ，变换取单位阵（A=I），输出 = 输入，机制完全关闭。\n\n关键参数初值与搜索范围：\n- n_hvg = 2500（固定）\n- n_pcs = 30（固定）\n- α（外推步长）= 0.3，搜索 [0.1, 0.5]，步长 0.1；若 X3 分数下降则减小到 0.1\n- k（低秩截断）= 10，搜索 [5, 15]\n- min_cells = 30（每阶段每型最少细胞数）\n- eigen_clip = 1e-4\n\nvec-score 快速筛选：先在小样本（每型最多 500 细胞、500 HVG）上跑通流程，确认协方差确实被修改（Frobenius 范数变化 > 1%）。然后全量跑，用 vec-score 查分。比较 α=0（关闭）和 α=0.3（开启）的四组分差异，重点看 covariation 和 de_recovery。若总分差异 < 2 分，重复查分 3 次取均值确认。",
  "expected_groups": ["covariation", "cell_state"],
  "risks": "1. 协方差变化主要是采样噪声而非发育趋势：两阶段间型内协方差差异可能来自有限的细胞数。缓解：Ledoit-Wolf 收缩 + 低秩截断（只取前 10 个方向）+ 小步长 α=0.3。Engineer 应在小样本上先检查 ΔΣ 的特征值谱是否有明显信号（前几个特征值远大于其余），若全为噪声量级则机制可能无效。2. 仿射变换产生极端表达值：裁剪特征值后应检查输出表达的范围，若超出输入范围 3 倍标准差则触发警告。3. proxy10 上机制完全关闭，分数 ≈ copy_last（~52.75），节点分数主要依赖 X3 表现；若 X3 无改善则总分不会超过基线。4. 细胞类型分组不准确（标签缺失或聚类不当）导致协方差估计混杂不同细胞状态。Engineer 应检查每个有效类型的细胞数分布。5. PCA 残差加回后，非主成分方向的共变结构不变，可能限制了改善幅度。这是低秩约束的固有代价。",
  "family_id": "lowrank_shape",
  "mechanism": "在低维 PCA 空间用 Ledoit-Wolf 收缩估计两阶段型内协方差，取协方差差异的低秩分量外推一小步，以受限仿射变换改变中心化残差（形状），均值不变。",
  "vs_constant_shift": "常数位移给每个型的所有细胞加同一个向量，改变的是群体均值位置，不影响型内分布形状；lowrank_shape 不移动均值（变换以型均值为中心），改变的是型内协方差结构（分布的宽度、方向、共变模式）。常数位移主要影响 DE 指标（基因排序），lowrank_shape 主要影响 variogram（共变结构）和 mmd_u（分布形状），理论上不改变 DE。",
  "mechanism_evidence": "Engineer 应检查：1）变换前后每个有效类型的均值变化（应 < 0.01 标准差，证明均值未移动）；2）变换前后型内协方差矩阵的 Frobenius 范数比（应偏离 1.0，证明形状确实改变）；3）ΔΣ 前 10 个特征值占总能量的比例（若 > 60% 说明低秩近似合理）；4）四组分各自变化：预期 covariation 改善、de_recovery 基本不变、direction 基本不变、cell_state 可能小幅变化；5）实际被修改的细胞数和比例（应为所有属于有效类型的细胞）。",
  "mechanism_off_control": "同一程序设 α=0（或等价地令 A=I 单位阵），此时 Σ_pred = Σ₁，变换为恒等映射，输出与输入完全相同。预期差别：关闭时输出 = 输入（逐元素相同），四组分与 copy_last 一致；开启时输出应不同，covariation 分组分数应有变化（改善或恶化）。若两者输出逐元素相同，说明机制未运行。"
}
```