v2431/FormulaShowcase
0
1# modules.py2# -*- coding: utf-8 -*-3from __future__ import annotations4 5from dataclasses import dataclass, field6from pathlib import Path7import pickle8from typing import Iterable, Optional, Set, Dict, Any, Tuple, List, Sequence9import collections10 11import networkx as nx12 13import netmedpy14import pandas as pd # type: ignore15import numpy as np16from collections import deque17from tqdm import tqdm18from concurrent.futures import ThreadPoolExecutor, as_completed19 20import pulp21import math22from pyscipopt import Model, quicksum23 24# from utils import _sample_preserving_degrees_fixed25# netmedpy.NetMedPy._sample_preserving_degrees = _sample_preserving_degrees_fixed26 27# ----------------------------------------------------------------------28# Proximity 部分（Step 2）29# ----------------------------------------------------------------------30@dataclass31class ProximityResult:32 """33 存放网络邻近度计算结果的数据类（支持单模块和多模块场景）。34 35 Attributes36 ----------37 module_name : str38 模块名称或场景名称39 full_stats : dict40 NetMedPy screening 返回的完整结果字典41 """42 43 module_name: str44 full_stats: Dict[str, Any]45 thres_z: Optional[float] = None46 47 def __post_init__(self):48 """自动构建 z_scores_by_module"""49 z_df = self.full_stats.get("z_score")50 if z_df is None or not isinstance(z_df, pd.DataFrame):51 raise ValueError("full_stats 中未找到 'z_score' DataFrame")52 self._z_scores_by_module = z_df.apply(pd.to_numeric, errors="coerce")53 54 @property55 def z_scores_by_module(self) -> pd.DataFrame:56 """返回多模块的 z-score 矩阵（行=候选实体，列=模块）"""57 return self._z_scores_by_module58 59 def aggregate_z_scores(self, mode: str = "min") -> pd.Series:60 """61 聚合多模块的 z-score 为场景级 z-score。62 63 Parameters64 ----------65 mode : str, default "min"66 聚合模式：'min', 'mean', 'max'67 68 Returns69 -------70 pd.Series71 聚合后的 z-score（index=候选实体id）72 """73 if mode == "min":74 z_agg = self._z_scores_by_module.min(axis=1)75 elif mode == "mean":76 z_agg = self._z_scores_by_module.mean(axis=1)77 elif mode == "max":78 z_agg = self._z_scores_by_module.max(axis=1)79 else:80 raise ValueError(f"不支持的聚合模式: {mode}，请使用 'min', 'mean', 或 'max'")81 82 z_agg.name = self.module_name83 return z_agg84 85 def filter_by_proximity(86 self,87 ns_targets: Dict[str, Iterable[str]],88 thres_z: float = -1.6,89 agg_mode: str = "min",90 max_candidates: Optional[int] = None,91 ) -> Dict[str, Iterable[str]]:92 """93 基于聚合 z-score 筛选天然来源。94 95 Parameters96 ----------97 ns_targets : dict98 天然来源 id -> 靶标蛋白列表99 thres_z : float, default -1.6100 z-score 阈值（越小越好）101 agg_mode : str, default "min"102 聚合模式：'min', 'mean', 'max'103 max_candidates : int or None, default None104 最大候选数量限制。如果不为 None，在阈值筛选后按 z-score 排序取前 N 个105 106 Returns107 -------108 dict109 筛选后的天然来源字典110 """111 self.thres_z = thres_z112 z_agg = self.aggregate_z_scores(mode=agg_mode)113 selected = z_agg[z_agg <= thres_z].index.tolist()114 115 # 如果指定了 max_candidates，进一步限制数量（保留 z 最小的前 N 个）116 if max_candidates is not None and len(selected) > max_candidates:117 z_selected = z_agg.loc[selected].sort_values()118 selected = z_selected.index[:max_candidates].tolist()119 print(f"[ProximityResult] Filtered by z<={thres_z} ({agg_mode}): {len(z_agg[z_agg <= thres_z])} candidates, limited to top {max_candidates}")120 else:121 print(f"[ProximityResult] Filtered by z<={thres_z} ({agg_mode}): {len(selected)} candidates")122 123 # 同步更新 z_scores_by_module，只保留筛选后的候选124 self._z_scores_by_module = self._z_scores_by_module.loc[selected].copy()125 126 return {ns: ns_targets[ns] for ns in selected if ns in ns_targets}127 128 @property129 def significant_hits(self) -> Dict[str, List[str]]:130 """返回每个候选实体显著命中的 module 列表。131 Returns132 -------133 Dict[str, List[str]]134 {candidate_id: [module_name1, module_name2, ...]}135 """136 z_df = self.z_scores_by_module137 result: Dict[str, List[str]] = {}138 139 for cid in z_df.index:140 hit_modules = []141 for module_name in z_df.columns:142 z = z_df.loc[cid, module_name]143 if pd.notna(z) and z <= self.thres_z:144 hit_modules.append(module_name)145 if hit_modules:146 result[cid] = hit_modules147 148 return result149 150 151class ProximityAnalyzer:152 """153 基于 NetMedPy 的网络邻近度计算器。154 155 主要职责：156 ----------157 1. 预计算并缓存 PPI 的距离矩阵 (all_pair_distances)158 2. 使用 screening 统一计算一批候选实体相对于指定模块的 proximity z-score159 """160 161 def __init__(162 self,163 ppi_graph: nx.Graph,164 *,165 distance_metric: str = "shortest_path",166 null_model: str = "soft_log_binning",167 repeat_allowed: bool = False,168 n_iter: int = 1000,169 n_procs: int = 1,170 bin_size: int = 100,171 ) -> None:172 """173 Parameters174 ----------175 ppi_graph : networkx.Graph176 PPI 子图177 distance_metric : str, default "shortest_path"178 NetMedPy 支持的距离度量：179 - "shortest_path"180 - "random_walk"181 - "biased_random_walk"182 - "communicability"183 等184 null_model : str, default "soft_log_binning"185 空模型类型（与 LCC 一致）186 repeat_allowed : bool, default False187 空模型采样时是否允许重复选择节点，True则只允许松弛一次bin，False则反复扩大bin188 n_iter : int, default 1000189 空模型随机采样次数190 n_procs : int, default 1191 并行进程数192 """193 194 self.ppi_graph = ppi_graph195 self.distance_metric = distance_metric196 self.null_model = null_model197 self.repeat_allowed = repeat_allowed198 self.n_iter = n_iter199 self.n_procs = n_procs200 self.bin_size = bin_size201 202 # 距离矩阵缓存203 self._distance_matrix = None204 205 # ----------------------------206 # 内部：懒加载距离矩阵207 # ----------------------------208 def _ensure_distance_matrix(self):209 """210 懒加载 PPI 的全对距离矩阵，仅计算一次。211 """212 if self._distance_matrix is None:213 self._distance_matrix = netmedpy.all_pair_distances(214 self.ppi_graph,215 distance=self.distance_metric,216 n_processors=self.n_procs,217 )218 return self._distance_matrix219 220 def load_distance_matrix(self, path: str):221 """222 从指定路径加载预计算的距离矩阵。223 224 Parameters225 ----------226 path : str227 距离矩阵文件的路径 （pkl）228 """229 self._distance_matrix = pickle.load(open(path, "rb"))230 231 @property232 def distance_matrix(self):233 """234 获取距离矩阵。235 236 Returns237 -------238 pd.DataFrame239 PPI 的全对距离矩阵240 """241 if self._distance_matrix is None:242 self._distance_matrix = self._ensure_distance_matrix()243 return self._distance_matrix244 245 # ----------------------------246 # 核心 API：单模块247 # ----------------------------248 def compute_for_sources(249 self,250 source_sets: Dict[str, Iterable[str]],251 module_nodes: Iterable[str],252 module_name: str = "module",253 ) -> ProximityResult:254 """255 计算一批候选实体相对于某模块的 proximity z-score（单模块版）。256 257 Parameters258 ----------259 source_sets : dict260 key: 候选实体 id（如 natural_source id）261 value: 对应的蛋白靶标列表 / 可迭代262 module_nodes : Iterable[str]263 模块的蛋白节点集合（建议使用 LCCResult.lcc_nodes）264 module_name : str, default "module"265 模块名称（用于 NetMedPy 输出的行名）266 267 Returns268 -------269 ProximityResult270 """271 graph_nodes = set(self.ppi_graph.nodes())272 273 # 1) 清洗模块节点：只保留存在于 PPI 的节点274 module_nodes_set = set(module_nodes) & graph_nodes275 if len(module_nodes_set) == 0:276 raise ValueError("module_nodes 在 PPI 中没有任何有效节点，无法计算 proximity。")277 278 # 2) 清洗 source_sets：只保留至少有 1 个节点在 PPI 中的候选实体279 cleaned_sources: Dict[str, list] = {}280 for sid, nodes in source_sets.items():281 nodes_in_graph = list(set(nodes) & graph_nodes)282 if len(nodes_in_graph) > 0:283 cleaned_sources[sid] = nodes_in_graph284 285 if not cleaned_sources:286 raise ValueError("source_sets 中没有任何候选实体在 PPI 中有有效靶标。")287 288 # 3) 准备 NetMedPy 所需输入格式289 distance_matrix = self.distance_matrix290 module_dict = {module_name: list(module_nodes_set)} # 列：模块291 292 # 4) 调用 screening 计算 proximity z-score293 screen_data = netmedpy.screening(294 cleaned_sources,295 module_dict,296 self.ppi_graph,297 distance_matrix,298 score="proximity",299 properties=[300 "z_score",301 "p_value_single_tail",302 "p_value_double_tail",303 "raw_amspl",304 ],305 null_model=self.null_model,306 repeat_allowed=self.repeat_allowed,307 n_iter=self.n_iter,308 bin_size=self.bin_size,309 n_procs=self.n_procs,310 )311 312 z_df: pd.DataFrame = screen_data["z_score"]313 314 if module_name not in z_df.columns:315 raise RuntimeError(f"NetMedPy screening result does not contain column: {module_name}")316 317 # 取出这一列，就是每个候选实体的 z-score318 z_series: pd.Series = z_df[module_name]319 z_series = pd.to_numeric(z_series, errors="coerce")320 321 return ProximityResult(322 module_name=module_name,323 z_scores=z_series,324 full_stats=screen_data,325 )326 327 # ----------------------------328 # 核心 API：多模块329 # ----------------------------330 def compute_for_modules(331 self,332 source_sets: Dict[str, Iterable[str]],333 modules: Dict[str, Iterable[str]],334 scenario_name: str = "scenario",335 ) -> Dict[str, Any]:336 """337 多模块版本的 proximity 计算，直接返回 NetMedPy 的 screen_data。338 339 Parameters340 ----------341 source_sets : dict342 key: 候选实体 id343 value: 对应的蛋白靶标列表344 modules : dict345 key: 模块名346 value: 模块蛋白集合347 scenario_name : str, default "scenario"348 场景名称（未使用，保留用于日志）349 350 Returns351 -------352 dict353 NetMedPy screening 的完整结果字典354 """355 graph_nodes = set(self.ppi_graph.nodes())356 357 # 清洗模块358 cleaned_modules: Dict[str, List[str]] = {}359 for mname, nodes in modules.items():360 nodes_in_graph = list(set(nodes) & graph_nodes)361 if nodes_in_graph:362 cleaned_modules[mname] = nodes_in_graph363 364 if not cleaned_modules:365 raise ValueError("modules 中没有任何模块在 PPI 中有有效节点。")366 367 # 清洗候选实体368 cleaned_sources: Dict[str, List[str]] = {}369 for sid, nodes in source_sets.items():370 nodes_in_graph = list(set(nodes) & graph_nodes)371 if len(nodes_in_graph) > 0:372 cleaned_sources[sid] = nodes_in_graph373 374 if not cleaned_sources:375 raise ValueError("source_sets 中没有任何候选实体在 PPI 中有有效靶标。")376 377 distance_matrix = self.distance_matrix378 379 screen_data = netmedpy.screening(380 cleaned_sources,381 cleaned_modules,382 self.ppi_graph,383 distance_matrix,384 score="proximity",385 properties=[386 "z_score",387 "p_value_single_tail",388 "p_value_double_tail",389 "raw_amspl",390 ],391 null_model=self.null_model,392 n_iter=self.n_iter,393 n_procs=self.n_procs,394 )395 396 print('[ProximityAnalyzer] Proximity screening finished')397 return screen_data398 399 400# ----------------------------------------------------------------------401# Separation 部分（Step 3）402# ----------------------------------------------------------------------403class SeparationAnalyzer:404 """405 基于 NetMedPy 的网络分离度（separation）计算器。406 407 用途（Step 3）：408 -------------409 - 给定一组候选实体（草药 / 食物） -> protein 靶标集合，410 计算它们之间的 pairwise separation（默认 source_sets 和 target_sets 相同），411 并返回 separation 的矩阵。412 """413 414 def __init__(415 self,416 ppi_graph: nx.Graph,417 distance_matrix,418 n_procs: int = 4,419 ) -> None:420 """421 Parameters422 ----------423 ppi_graph : nx.Graph424 与 Step 1/2 相同的 PPI 网络。425 distance_matrix :426 netmedpy.all_pair_distances 返回的 DistanceMatrix 对象。427 n_procs : int428 并行进程数。429 """430 self.ppi_graph = ppi_graph431 self.distance_matrix = distance_matrix432 self.n_procs = n_procs433 434 def _clean_sets(435 self,436 sets: Dict[str, Iterable[str]],437 ) -> Dict[str, list]:438 """439 辅助函数：把每个实体的 protein 集合限制在 PPI 节点里，440 丢掉在 PPI 中完全没有靶标的实体。441 """442 graph_nodes = set(self.ppi_graph.nodes())443 cleaned: Dict[str, list] = {}444 for sid, nodes in sets.items():445 nodes_in_graph = list(set(nodes) & graph_nodes)446 if nodes_in_graph:447 cleaned[sid] = nodes_in_graph448 return cleaned449 450 def compute_pairwise(451 self,452 source_sets: Dict[str, Iterable[str]],453 target_sets: Optional[Dict[str, Iterable[str]]] = None,454 ):455 """456 计算候选实体之间（或 source vs target）的 pairwise separation。457 458 Parameters459 ----------460 source_sets : dict461 {实体id -> [protein 节点列表]}，例如 {herb1: [P1, P2, ...], herb2: [...]}。462 target_sets : dict or None463 如果为 None，则默认 target_sets = source_sets，464 即计算所有候选之间的两两 separation。465 466 Returns467 -------468 dict469 NetMedPy screening 的原始结果。470 """471 # 1) 清洗输入472 cleaned_sources = self._clean_sets(source_sets)473 if not cleaned_sources:474 raise ValueError("source_sets 中没有任何实体在 PPI 中有有效靶标。")475 476 if target_sets is None:477 cleaned_targets = cleaned_sources478 else:479 cleaned_targets = self._clean_sets(target_sets)480 if not cleaned_targets:481 raise ValueError("target_sets 中没有任何实体在 PPI 中有有效靶标。")482 483 # 2) 直接调用 NetMedPy 的 screening(score=\"separation\")484 screen_data = netmedpy.screening(485 cleaned_sources,486 cleaned_targets,487 self.ppi_graph,488 self.distance_matrix,489 score="separation",490 properties=["raw_separation"],491 n_procs=self.n_procs,492 )493 494 return screen_data495 496# ----------------------------------------------------------------------497# 组方优化 部分（Step 4）- Unified FormulaOptimizer (SCIP / CBC / GRASP)498# ----------------------------------------------------------------------499 500@dataclass501class FormulaSolution:502 """一个候选组方解（支持同构/异构 size）。"""503 selected: List[str]504 objective: float505 506 # size 信息507 size: int508 size_by_group: Optional[Dict[str, int]] = None509 510 # 一些可读指标（不一定等同于优化目标）511 mean_z: float = float("nan")512 mean_separation: float = float("nan")513 coverage: float = float("nan") # [0,1] 或 weighted 版本514 515 # 目标分解（便于 debug/报告）516 obj_components: Dict[str, float] = field(default_factory=dict)517 518 # 其他元信息（例如：method, seed, scheme_id, raw_obj）519 meta: Dict[str, Any] = field(default_factory=dict)520 521 522class FormulaOptimizer:523 """524 Step4：多味组方组合优化器（SCIP / CBC(CBC) / GRASP）。525 - scip: PySCIPOpt + SCIP，线性加性目标 + 近优域 + overlap 多样化526 - cbc : PuLP + CBC，模块级 proximity(min-linearization) + no-good(overlap) 多样化527 - grasp: 子模 GRASP（log1p prox + facility-location diversity + coverage）528 529 说明：530 - optimize() 只保留“公共参数”；求解器/多样性等参数通过 optimizer_setup() 预先配置。531 - target_size 支持：532 * int 或 list[int]：同构总味数533 * dict[str,int] 或 list[dict]：异构配额（需 source_groups: candidate -> set[group]）534 """535 536 def __init__(537 self,538 *,539 prox_result: Any,540 sep_result: Dict[str, Any],541 ns_targets: Dict[str, Iterable[str]],542 module_nodes: Set[str],543 modules: Optional[Dict[str, Set[str]]] = None,544 agg_mode: str = "min",545 ) -> None:546 self.prox_result = prox_result547 self.sep_result = sep_result548 self.ns_targets = {k: list(v) for k, v in ns_targets.items()}549 self.module_nodes: Set[str] = set(module_nodes)550 self.modules = modules551 self.agg_mode = agg_mode552 553 # --- z-score matrix ---554 z_df = getattr(prox_result, "z_scores_by_module", None)555 if z_df is None:556 raise ValueError("prox_result 必须包含 z_scores_by_module 属性（DataFrame）。")557 if not isinstance(z_df, pd.DataFrame):558 raise TypeError("prox_result.z_scores_by_module 必须是 pandas.DataFrame。")559 self.z_scores_by_module: pd.DataFrame = z_df560 561 # --- separation matrix ---562 if "raw_separation" not in sep_result:563 raise ValueError("sep_result 必须包含 key='raw_separation' 的 DataFrame。")564 self.sep_df: pd.DataFrame = sep_result["raw_separation"]565 if not isinstance(self.sep_df, pd.DataFrame):566 raise TypeError("sep_result['raw_separation'] 必须是 pandas.DataFrame。")567 568 # --- solver config ---569 self._method: str = "scip" # default570 self._cfg: Dict[str, Any] = {571 # multi-solution quality/diversity572 "quality_rel_gap": 0.10,573 "quality_abs_gap": None,574 "min_hamming": 2,575 576 # scip/cbc577 "time_limit": 30.0,578 "mip_rel_gap": 0.01,579 "threads": 0,580 "verbose": False,581 582 # grasp583 "grasp_pool_factor": 50,584 "rcl_size": 10,585 "random_seed": 2025,586 587 # cbc (CBC no-good overlap ratio)588 "diversity_ratio": 0.0, # 0.0 -> 至少差 1 味；0.3 -> 至少 30% 不同589 }590 591 # ============================592 # Setup API593 # ============================594 595 def optimizer_setup(self, *, method: str = "scip", **solver_params: Any) -> None:596 """597 在调用 optimize() 之前配置求解器与其参数。598 599 Parameters600 ----------601 method:602 "scip" | "cbc" | "grasp"603 - 兼容旧写法：method="ilp" 会被映射到 "scip"604 solver_params:605 - 通用：quality_rel_gap, quality_abs_gap, min_hamming, verbose606 - scip/cbc：time_limit, mip_rel_gap, threads607 - grasp：grasp_pool_factor, rcl_size, random_seed608 - cbc：diversity_ratio609 """610 m = str(method).lower().strip()611 if m == "ilp":612 m = "scip"613 if m not in {"scip", "cbc", "grasp"}:614 raise ValueError("method 仅支持 'scip' / 'cbc' / 'grasp'（兼容 'ilp'->'scip'）。")615 self._method = m616 617 for k, v in solver_params.items():618 self._cfg[k] = v619 620 # ============================621 # Public API622 # ============================623 624 def optimize(625 self,626 *,627 candidate_ids: Optional[Sequence[str]] = None,628 target_size: Optional[Any] = None,629 source_groups: Optional[Dict[str, Set[str]]] = None,630 n_solutions: int = 3,631 lambda_prox: float = 1.0,632 lambda_sep: float = 1.0,633 lambda_neg_sep: float = 0.0,634 lambda_cov: float = 0.0,635 weighted_coverage: bool = False,636 node_weights: Optional[Dict[str, float]] = None,637 module_weights: Optional[Dict[str, float]] = None,638 ) -> List[FormulaSolution]:639 """返回多样解列表（长度 <= n_solutions * len(size_schemes or 1)）。"""640 641 # 提前读取 setup config，保证 verbose 等变量在后续任何地方可用642 method = self._method643 cfg = self._cfg644 verbose = bool(cfg.get("verbose", False))645 646 quality_rel_gap = float(cfg.get("quality_rel_gap", 0.10))647 quality_abs_gap = cfg.get("quality_abs_gap", None)648 quality_abs_gap = float(quality_abs_gap) if quality_abs_gap is not None else None649 min_hamming = int(cfg.get("min_hamming", 2))650 651 time_limit = float(cfg.get("time_limit", 30.0))652 mip_rel_gap = float(cfg.get("mip_rel_gap", 0.01))653 threads = int(cfg.get("threads", 0))654 655 grasp_pool_factor = int(cfg.get("grasp_pool_factor", 50))656 rcl_size = int(cfg.get("rcl_size", 10))657 random_seed = int(cfg.get("random_seed", 2025))658 659 diversity_ratio = float(cfg.get("diversity_ratio", 0.0))660 661 if candidate_ids is None:662 candidate_ids = list(self.ns_targets.keys())663 candidate_ids = list(candidate_ids)664 if not candidate_ids:665 raise ValueError("candidate_ids 为空。")666 667 # 保证候选存在于矩阵中668 candidate_ids = [cid for cid in candidate_ids if cid in self.z_scores_by_module.index]669 if not candidate_ids:670 raise ValueError("candidate_ids 与 z_scores_by_module 无交集。")671 672 # separation 可能只包含子集673 row_ids = set(self.sep_df.index)674 col_ids = set(self.sep_df.columns)675 candidate_ids = [cid for cid in candidate_ids if (cid in row_ids) or (cid in col_ids)]676 if not candidate_ids:677 raise ValueError("候选草药在 separation 矩阵中不存在有效的行/列。")678 679 # ---------- 解析 target_size（支持 int / list[int] / dict / list[dict]） ----------680 use_group_sizes = False681 schemes: List[Optional[Dict[str, int]]] = []682 size_list: List[int] = []683 684 if isinstance(target_size, dict):685 use_group_sizes = True686 schemes = [self._normalize_size_dict(target_size)]687 elif isinstance(target_size, (list, tuple)) and target_size and all(isinstance(x, dict) for x in target_size):688 use_group_sizes = True689 schemes = [self._normalize_size_dict(s) for s in target_size]690 691 if use_group_sizes:692 if source_groups is None:693 raise ValueError("使用 dict 形式的 target_size 时必须提供 source_groups 映射。")694 else:695 if isinstance(target_size, (list, tuple, set, np.ndarray)):696 size_list = [int(s) for s in target_size if int(s) > 0]697 size_list = sorted(set(size_list))698 elif isinstance(target_size, (int, np.integer)):699 if int(target_size) > 0:700 size_list = [int(target_size)]701 else:702 raise ValueError("target_size 必须是 int, list[int], dict 或 list[dict]。")703 704 if not size_list:705 raise ValueError("无法从 target_size 确定需要优化的味数。")706 707 schemes = [None] * len(size_list)708 709 # ---------- 保险机制：调整 size 不超过候选数量 ----------710 if use_group_sizes:711 # 异构配额：统计每个 group 的候选数量712 group_candidate_count: Dict[str, int] = {}713 for cid in candidate_ids:714 groups = source_groups.get(cid, set())715 for g in groups:716 group_candidate_count[g] = group_candidate_count.get(g, 0) + 1717 718 # 调整每个 scheme 的配额719 adjusted_schemes: List[Dict[str, int]] = []720 for scheme in schemes:721 adjusted_scheme: Dict[str, int] = {}722 total_adjusted = 0723 for g, requested in scheme.items():724 available = group_candidate_count.get(g, 0)725 adjusted = min(requested, available)726 if adjusted > 0:727 adjusted_scheme[g] = adjusted728 total_adjusted += adjusted729 if adjusted < requested and verbose:730 print(f"[FormulaOptimizer] Safety mechanism: group '{g}' requested {requested} items, but only {available} candidates available, adjusted to {adjusted}")731 732 if total_adjusted > 0:733 adjusted_schemes.append(adjusted_scheme)734 elif verbose:735 print(f"[FormulaOptimizer] Safety mechanism: scheme {scheme} adjusted to total of 0, skipped")736 737 schemes = adjusted_schemes738 if not schemes:739 print("[FormulaOptimizer] All schemes adjusted to no available candidates, returning empty solution.")740 return []741 else:742 # 同构总味数：直接用 min(size, len(candidate_ids))743 n_candidates = len(candidate_ids)744 adjusted_size_list: List[int] = []745 for s in size_list:746 adjusted = min(s, n_candidates)747 if adjusted > 0:748 adjusted_size_list.append(adjusted)749 if adjusted < s and verbose:750 print(f"[FormulaOptimizer] Safety mechanism: requested {s} items, but only {n_candidates} candidates available, adjusted to {adjusted}")751 752 size_list = adjusted_size_list753 if not size_list:754 print("[FormulaOptimizer] All sizes adjusted to 0, returning empty solution.")755 return []756 757 schemes = [None] * len(size_list)758 759 # 构建 hit 映射760 herb_to_hits, gene_to_herbs = self._build_herb_gene_maps(candidate_ids)761 762 module_size = len(self.module_nodes)763 if module_size <= 0:764 module_size = max(1, len(gene_to_herbs))765 766 # separation lookup767 pos_sep, neg_sep = self._build_sep_lookup(candidate_ids)768 769 # 聚合 z（用于可读指标）770 agg_z = self._aggregate_z_scores(candidate_ids, mode=self.agg_mode)771 772 all_solutions: List[FormulaSolution] = []773 774 for scheme_idx, scheme in enumerate(schemes):775 k_val = None if use_group_sizes else (size_list[scheme_idx] if scheme_idx < len(size_list) else None)776 777 if method == "scip":778 sols = self._solve_scip_multi(779 candidate_ids=candidate_ids,780 scheme=scheme,781 k=k_val,782 n_solutions=n_solutions,783 quality_rel_gap=quality_rel_gap,784 quality_abs_gap=quality_abs_gap,785 min_hamming=min_hamming,786 lambda_prox=lambda_prox,787 lambda_sep=lambda_sep,788 lambda_neg_sep=lambda_neg_sep,789 lambda_cov=lambda_cov,790 weighted_coverage=weighted_coverage,791 node_weights=node_weights,792 module_weights=module_weights,793 herb_to_hits=herb_to_hits,794 gene_to_herbs=gene_to_herbs,795 module_size=module_size,796 pos_sep=pos_sep,797 neg_sep=neg_sep,798 agg_z=agg_z,799 source_groups=source_groups,800 time_limit=time_limit,801 mip_rel_gap=mip_rel_gap,802 threads=threads,803 verbose=verbose,804 scheme_idx=scheme_idx,805 )806 all_solutions.extend(sols)807 808 elif method == "cbc":809 sols = self._solve_cbc_multi(810 candidate_ids=candidate_ids,811 scheme=scheme,812 k=k_val,813 n_solutions=n_solutions,814 quality_rel_gap=quality_rel_gap,815 quality_abs_gap=quality_abs_gap,816 diversity_ratio=diversity_ratio,817 lambda_prox=lambda_prox,818 lambda_sep=lambda_sep,819 lambda_neg_sep=lambda_neg_sep,820 lambda_cov=lambda_cov,821 weighted_coverage=weighted_coverage,822 node_weights=node_weights,823 module_weights=module_weights,824 herb_to_hits=herb_to_hits,825 gene_to_herbs=gene_to_herbs,826 module_size=module_size,827 pos_sep=pos_sep,828 neg_sep=neg_sep,829 agg_z=agg_z,830 source_groups=source_groups,831 time_limit=time_limit,832 mip_rel_gap=mip_rel_gap,833 threads=threads,834 verbose=verbose,835 scheme_idx=scheme_idx,836 )837 all_solutions.extend(sols)838 839 elif method == "grasp":840 rng = np.random.default_rng(int(random_seed))841 sols = self._solve_grasp_multi(842 candidate_ids=candidate_ids,843 scheme=scheme,844 k=k_val,845 n_solutions=n_solutions,846 quality_rel_gap=quality_rel_gap,847 quality_abs_gap=quality_abs_gap,848 min_hamming=min_hamming,849 lambda_prox=lambda_prox,850 lambda_sep=lambda_sep,851 lambda_cov=lambda_cov,852 weighted_coverage=weighted_coverage,853 node_weights=node_weights,854 module_weights=module_weights,855 herb_to_hits=herb_to_hits,856 gene_to_herbs=gene_to_herbs,857 module_size=module_size,858 agg_z=agg_z,859 pos_sep=pos_sep,860 source_groups=source_groups,861 rng=rng,862 pool_factor=grasp_pool_factor,863 rcl_size=rcl_size,864 verbose=verbose,865 scheme_idx=scheme_idx,866 )867 all_solutions.extend(sols)868 869 else:870 raise RuntimeError(f"Unknown method '{method}'")871 872 return all_solutions873 874 # ============================875 # Common helpers876 # ============================877 878 def _normalize_size_dict(self, d: Dict[str, int]) -> Dict[str, int]:879 out: Dict[str, int] = {}880 for k, v in d.items():881 vv = int(v)882 if vv > 0:883 out[str(k)] = vv884 if not out:885 raise ValueError("size scheme 为空或全为非正数。")886 return out887 888 def _aggregate_z_scores(self, candidate_ids: List[str], mode: str = "min") -> pd.Series:889 sub = self.z_scores_by_module.loc[candidate_ids].copy()890 sub = sub.replace([np.inf, -np.inf], np.nan).fillna(0.0)891 if mode == "min":892 s = sub.min(axis=1)893 elif mode == "mean":894 s = sub.mean(axis=1)895 elif mode == "max":896 s = sub.max(axis=1)897 else:898 raise ValueError("agg_mode 仅支持 'min'/'mean'/'max'。")899 s.name = "agg_z"900 return s901 902 def _get_sep_value(self, cid_i: str, cid_j: str) -> float:903 try:904 if cid_i in self.sep_df.index and cid_j in self.sep_df.columns:905 return float(self.sep_df.loc[cid_i, cid_j])906 if cid_j in self.sep_df.index and cid_i in self.sep_df.columns:907 return float(self.sep_df.loc[cid_j, cid_i])908 except Exception:909 pass910 return float("nan")911 912 def _build_sep_lookup(self, candidate_ids: List[str]) -> Tuple[Dict[Tuple[str, str], float], Dict[Tuple[str, str], float]]:913 pos_sep: Dict[Tuple[str, str], float] = {}914 neg_sep: Dict[Tuple[str, str], float] = {}915 for i, cid_i in enumerate(candidate_ids):916 for j in range(i + 1, len(candidate_ids)):917 cid_j = candidate_ids[j]918 val = self._get_sep_value(cid_i, cid_j)919 if not np.isfinite(val):920 continue921 if val > 0:922 pos_sep[(cid_i, cid_j)] = float(val)923 elif val < 0:924 neg_sep[(cid_i, cid_j)] = float(-val) # store magnitude925 return pos_sep, neg_sep926 927 def _build_herb_gene_maps(self, candidate_ids: List[str]) -> Tuple[Dict[str, List[str]], Dict[str, List[str]]]:928 herb_to_hits: Dict[str, List[str]] = {}929 gene_to_herbs: Dict[str, List[str]] = {}930 module_nodes = self.module_nodes931 932 for cid in candidate_ids:933 hits = [g for g in self.ns_targets.get(cid, []) if (not module_nodes) or (g in module_nodes)]934 herb_to_hits[cid] = hits935 for g in hits:936 gene_to_herbs.setdefault(g, []).append(cid)937 938 for g, herbs in gene_to_herbs.items():939 gene_to_herbs[g] = sorted(set(herbs))940 941 return herb_to_hits, gene_to_herbs942 943 def _compute_readable_metrics(944 self,945 *,946 selected: List[str],947 agg_z: pd.Series,948 herb_to_hits: Dict[str, List[str]],949 module_size: int,950 weighted_coverage: bool,951 node_weights: Optional[Dict[str, float]],952 ) -> Tuple[float, float, float]:953 z_vals = agg_z.loc[selected].values if len(selected) > 0 else np.array([])954 mean_z = float(np.mean(z_vals)) if len(z_vals) > 0 else float("nan")955 956 sep_vals: List[float] = []957 for i in range(len(selected)):958 for j in range(i + 1, len(selected)):959 v = self._get_sep_value(selected[i], selected[j])960 if np.isfinite(v):961 sep_vals.append(float(v))962 mean_sep = float(np.mean(sep_vals)) if sep_vals else float("nan")963 964 hit_genes: Set[str] = set()965 for cid in selected:966 hit_genes.update(herb_to_hits.get(cid, []))967 968 if weighted_coverage and node_weights is not None and len(hit_genes) > 0:969 denom = float(sum(node_weights.get(g, 1.0) for g in self.module_nodes)) or 1.0970 num = float(sum(node_weights.get(g, 1.0) for g in hit_genes))971 cov = num / denom972 else:973 cov = (len(hit_genes) / float(module_size)) if module_size > 0 else float("nan")974 975 return mean_z, mean_sep, cov976 977 # ============================978 # GRASP (multi-solution) - unchanged from modules1979 # ============================980 981 def _solve_grasp_multi(982 self,983 *,984 candidate_ids: List[str],985 scheme: Optional[Dict[str, int]],986 k: Optional[int],987 n_solutions: int,988 quality_rel_gap: float,989 quality_abs_gap: Optional[float],990 min_hamming: int,991 lambda_prox: float,992 lambda_sep: float,993 lambda_cov: float,994 weighted_coverage: bool,995 node_weights: Optional[Dict[str, float]],996 module_weights: Optional[Dict[str, float]],997 herb_to_hits: Dict[str, List[str]],998 gene_to_herbs: Dict[str, List[str]],999 module_size: int,1000 agg_z: pd.Series,1001 pos_sep: Dict[Tuple[str, str], float],1002 source_groups: Optional[Dict[str, Set[str]]],1003 rng: np.random.Generator,1004 pool_factor: int,1005 rcl_size: int,1006 verbose: bool,1007 scheme_idx: int,1008 ) -> List[FormulaSolution]:1009 if scheme is None:1010 if k is None:1011 raise ValueError("GRASP 同构 size 需要 k。")1012 k_total = int(k)1013 if k_total <= 0:1014 return []1015 else:1016 k_total = int(sum(scheme.values()))1017 if source_groups is None:1018 raise ValueError("GRASP 异构 size 需要 source_groups。")1019 1020 pool_size = max(int(n_solutions * pool_factor), 50)1021 pool: List[FormulaSolution] = []1022 1023 for t in range(pool_size):1024 seed_tag = int(rng.integers(0, 2**31 - 1))1025 sol_ids = self._grasp_construct(1026 candidate_ids=candidate_ids,1027 k_total=k_total,1028 scheme=scheme,1029 source_groups=source_groups,1030 lambda_prox=lambda_prox,1031 lambda_sep=lambda_sep,1032 lambda_cov=lambda_cov,1033 weighted_coverage=weighted_coverage,1034 node_weights=node_weights,1035 module_weights=module_weights,1036 herb_to_hits=herb_to_hits,1037 gene_to_herbs=gene_to_herbs,1038 module_size=module_size,1039 agg_z=agg_z,1040 rng=np.random.default_rng(seed_tag),1041 rcl_size=rcl_size,1042 )1043 if not sol_ids:1044 continue1045 1046 obj, comps = self._eval_submodular_objective(1047 selected=sol_ids,1048 candidate_ids=candidate_ids,1049 lambda_prox=lambda_prox,1050 lambda_sep=lambda_sep,1051 lambda_cov=lambda_cov,1052 weighted_coverage=weighted_coverage,1053 node_weights=node_weights,1054 module_weights=module_weights,1055 herb_to_hits=herb_to_hits,1056 gene_to_herbs=gene_to_herbs,1057 module_size=module_size,1058 )1059 mean_z, mean_sep, cov = self._compute_readable_metrics(1060 selected=sol_ids,1061 agg_z=agg_z,1062 herb_to_hits=herb_to_hits,1063 module_size=module_size,1064 weighted_coverage=weighted_coverage,1065 node_weights=node_weights,1066 )1067 1068 pool.append(1069 FormulaSolution(1070 selected=sol_ids,1071 objective=float(obj),1072 size=len(sol_ids),1073 size_by_group=scheme,1074 mean_z=mean_z,1075 mean_separation=mean_sep,1076 coverage=cov,1077 obj_components=comps,1078 meta={"method": "grasp", "scheme_idx": scheme_idx, "seed": seed_tag},1079 )1080 )1081 1082 if not pool:1083 return []1084 1085 picked = self._select_diverse_subset(1086 pool=pool,1087 n=n_solutions,1088 min_hamming=min_hamming,1089 quality_rel_gap=quality_rel_gap,1090 quality_abs_gap=quality_abs_gap,1091 )1092 1093 if verbose:1094 print(f"[FormulaOptimizer][GRASP] scheme={scheme_idx} pool={len(pool)} picked={len(picked)}")1095 1096 return picked1097 1098 def _grasp_construct(1099 self,1100 *,1101 candidate_ids: List[str],1102 k_total: int,1103 scheme: Optional[Dict[str, int]],1104 source_groups: Optional[Dict[str, Set[str]]],1105 lambda_prox: float,1106 lambda_sep: float,1107 lambda_cov: float,1108 weighted_coverage: bool,1109 node_weights: Optional[Dict[str, float]],1110 module_weights: Optional[Dict[str, float]],1111 herb_to_hits: Dict[str, List[str]],1112 gene_to_herbs: Dict[str, List[str]],1113 module_size: int,1114 agg_z: pd.Series,1115 rng: np.random.Generator,1116 rcl_size: int,1117 ) -> List[str]:1118 if k_total <= 0:1119 return []1120 1121 use_group = scheme is not None1122 remaining: Optional[Dict[str, int]] = None1123 if use_group:1124 if source_groups is None:1125 raise ValueError("scheme 非空时必须提供 source_groups。")1126 remaining = {str(g): int(v) for g, v in scheme.items() if int(v) > 0}1127 if sum(remaining.values()) != k_total:1128 k_total = sum(remaining.values())1129 1130 def pick_group(cid: str) -> Optional[str]:1131 if not use_group or remaining is None:1132 return None1133 gs = source_groups.get(cid, set()) if source_groups is not None else set()1134 feasible = [str(g) for g in gs if str(g) in remaining and remaining[str(g)] > 0]1135 if not feasible:1136 return None1137 feasible.sort(key=lambda g: remaining.get(g, 0), reverse=True)1138 return feasible[0]1139 1140 z_df = self.z_scores_by_module1141 cols = list(z_df.columns)1142 sub = z_df.loc[candidate_ids, cols].replace([np.inf, -np.inf], np.nan).fillna(0.0)1143 a_mat = np.maximum(0.0, -sub.to_numpy(dtype=float))1144 1145 if module_weights is not None:1146 w = np.array([float(module_weights.get(c, 1.0)) for c in cols], dtype=float)1147 else:1148 w = np.ones(len(cols), dtype=float)1149 if not np.isfinite(w).all() or w.sum() <= 0.0:1150 w = np.ones(len(cols), dtype=float)1151 w = w / w.sum()1152 1153 n = len(candidate_ids)1154 idx_of = {cid: i for i, cid in enumerate(candidate_ids)}1155 sim_mat = np.zeros((n, n), dtype=float)1156 for i, cid_i in enumerate(candidate_ids):1157 for j in range(i + 1, n):1158 s = self._get_sep_value(cid_i, candidate_ids[j])1159 sim = 0.0 if (np.isnan(s) or s >= 0) else -float(s)1160 sim_mat[i, j] = sim1161 sim_mat[j, i] = sim1162 1163 selected: List[str] = []1164 selected_set: Set[str] = set()1165 1166 sum_a = np.zeros((len(cols),), dtype=float)1167 best_sim = np.zeros((n,), dtype=float)1168 1169 covered: Set[str] = set()1170 if weighted_coverage and node_weights is not None:1171 cov_total = float(sum(node_weights.get(g, 1.0) for g in self.module_nodes)) or 1.01172 cov_curr = 0.01173 else:1174 cov_total = float(module_size) or 1.01175 cov_curr = 0.01176 1177 for step in range(k_total):1178 gains: List[Tuple[float, str, Optional[str], Any]] = []1179 1180 for cid in candidate_ids:1181 if cid in selected_set:1182 continue1183 1184 g_assigned = pick_group(cid)1185 if use_group and g_assigned is None:1186 continue1187 1188 i = idx_of[cid]1189 1190 if lambda_prox != 0.0:1191 new_sum_a = sum_a + a_mat[i, :]1192 delta_prox = float(np.sum(w * (np.log1p(new_sum_a) - np.log1p(sum_a))))1193 else:1194 new_sum_a = sum_a1195 delta_prox = 0.01196 1197 genes = herb_to_hits.get(cid, [])1198 new_genes = [g for g in genes if g not in covered]1199 delta_cov = 0.01200 new_cov_curr = cov_curr