chore: 初始化 MCM/ICM Toolkit 仓库骨架
审查并精简 ChatGPT 的建库建议后落地的第一版骨架，只定接口不写实现。 - 建议审查：14 个包 → 7 个；约 250 个函数 → 约 40 个；5 个 Result 类 → 1 个 （docs/01-suggestion-review.md） - docs/：架构、路线图（含永久不做清单）、环境、Git 工作流、函数卡片模板 - src/mcm_toolkit/：core / data / stats / models / solve / decision / viz 共 23 个模块，均为签名 + docstring + NotImplementedError - tests/：9 个文件，其中 test_api_surface.py 已通过（20 passed） - examples/ problems/ notebooks/ figures/ tables/ 骨架 - pyproject.toml（核心依赖 4 个 + 6 组 extras）、environment.yml - scripts/：setup-venv.ps1（ASCII + LF）、setup-conda.ps1、doctor.ps1 - 环境：.venv（system-site-packages）与 conda 环境 mcm-toolkit 均已建好， 数值栈版本与 base 完全对齐（numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6，MKL 后端）
This commit is contained in:
guzhujushi committed
commit
5319589dfd
76 files changed
+6812
No files matched your search
@@ -0,0 +1,58 @@
|
||||
"""solve —— 求解：优化 / 模拟 / 敏感性 / 图 / 动力学。
|
||||
|
||||
竞赛主力包。原建议把 `optimization` `simulation` `sensitivity` 拆成三个包，
|
||||
但它们本质是同一件事的三种形态——"把模型反复算"——共享
|
||||
`Result`、`rng`、并行与缓存约定。合成一个包，这些约定只写一遍。
|
||||
|
||||
`graph` 与 `dynamics` 是 v0.2 内容（低频，且 NetworkX / SymPy 已覆盖大部分需求）。
|
||||
"""
|
||||
|
||||
from .optimize import (
|
||||
dynamic_programming,
|
||||
integer_program,
|
||||
knapsack_dp,
|
||||
linear_program,
|
||||
minimize,
|
||||
pareto_front,
|
||||
simulated_annealing,
|
||||
)
|
||||
from .sensitivity import (
|
||||
local_sensitivity,
|
||||
morris_analyze,
|
||||
one_at_a_time,
|
||||
parameter_sweep,
|
||||
robustness_test,
|
||||
sobol_analyze,
|
||||
)
|
||||
from .simulate import (
|
||||
markov_simulate,
|
||||
markov_steady_state,
|
||||
mm1_metrics,
|
||||
monte_carlo,
|
||||
monte_carlo_integral,
|
||||
monte_carlo_probability,
|
||||
random_walk,
|
||||
)
|
||||
|
||||
__all__ = [
|
||||
"minimize",
|
||||
"linear_program",
|
||||
"integer_program",
|
||||
"dynamic_programming",
|
||||
"knapsack_dp",
|
||||
"simulated_annealing",
|
||||
"pareto_front",
|
||||
"monte_carlo",
|
||||
"monte_carlo_probability",
|
||||
"monte_carlo_integral",
|
||||
"random_walk",
|
||||
"markov_simulate",
|
||||
"markov_steady_state",
|
||||
"mm1_metrics",
|
||||
"local_sensitivity",
|
||||
"one_at_a_time",
|
||||
"parameter_sweep",
|
||||
"robustness_test",
|
||||
"sobol_analyze",
|
||||
"morris_analyze",
|
||||
]
|
||||
@@ -0,0 +1,79 @@
|
||||
"""微分方程与动力系统（v0.2）。
|
||||
|
||||
原建议把 `dynamics/` 做成独立包。考虑到：
|
||||
1. 只有 MCM A 用得上，一年一题；
|
||||
2. 底层 SymPy（解析）+ SciPy（数值）已经覆盖全部需求；
|
||||
|
||||
降级为一个模块。依赖：`pip install "mcm-toolkit[symbolic]"`（sympy）。
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from collections.abc import Callable
|
||||
from typing import Any
|
||||
|
||||
__all__ = [
|
||||
"find_equilibrium",
|
||||
"local_stability",
|
||||
"solve_ode",
|
||||
"solve_system",
|
||||
]
|
||||
|
||||
|
||||
def solve_ode(
|
||||
derivative: Callable[..., Any],
|
||||
t_span: tuple[float, float],
|
||||
y0: Any,
|
||||
t_eval: Any = None,
|
||||
method: str = "RK45",
|
||||
**kwargs: Any,
|
||||
) -> Any:
|
||||
"""解 ODE 初值问题（底层 `scipy.integrate.solve_ivp`）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
method:
|
||||
`"RK45"` 通用；`"LSODA"` 刚性系统自动切换（种群/化学反应常见）；
|
||||
`"Radau"` / `"BDF"` 强刚性。
|
||||
|
||||
Notes
|
||||
-----
|
||||
**刚性是 ODE 最常踩的坑**：非刚性方法遇到刚性问题会跑到天荒地老
|
||||
且以为是自己代码慢。判据：解里同时存在快慢两个时间尺度。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = solve_ode(lambda t, y: -y, t_span=(0, 5), y0=[1.0]) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def solve_system(equations: Any, variables: Any, **kwargs: Any) -> Any:
|
||||
"""解方程组（底层 `sympy.solve` / `scipy.optimize.root`）。
|
||||
|
||||
`equations` 是符号表达式列表时走 SymPy 解析解，是数值函数时走 root 数值解。
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def find_equilibrium(f: Callable[..., Any], x0: Any, n_starts: int = 10, random_state: int | None = None) -> Any:
|
||||
"""找平衡点。
|
||||
|
||||
**必须多初值搜索**：单个初值只能找到最近的平衡点，
|
||||
而"系统有几个平衡点、分别是什么类型"才是题目真正要问的。
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def local_stability(jacobian_matrix: Any, equilibrium: Any = None) -> Any:
|
||||
"""按 Jacobian 的特征值判定平衡点的局部稳定性。
|
||||
|
||||
判据（连续系统）：实部全负 → 稳定结点/焦点；有正实部 → 不稳定；
|
||||
纯虚根 → 中心（临界情形，线性化失效，需要非线性分析）。
|
||||
`equilibrium` 用于把特征值和解的位置一起报出来。
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
@@ -0,0 +1,86 @@
|
||||
"""图与网络（v0.2）。
|
||||
|
||||
原建议把 `network/` 做成一个独立包（5 个子模块、20+ 函数），
|
||||
但里面绝大多数是 NetworkX 的一行转发，且只有 ICM D 用得上。
|
||||
降级为一个模块，只留真正有封装价值的入口。
|
||||
|
||||
依赖：`pip install "mcm-toolkit[graph]"`（networkx）。
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from typing import Any
|
||||
|
||||
__all__ = [
|
||||
"betweenness_centrality",
|
||||
"build_graph",
|
||||
"community_detection",
|
||||
"max_flow",
|
||||
"min_spanning_tree",
|
||||
"pagerank",
|
||||
"shortest_path",
|
||||
]
|
||||
|
||||
|
||||
def build_graph(edges: Any, directed: bool = False, weighted: bool = True) -> Any:
|
||||
"""从边表建图。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
edges:
|
||||
`[(u, v), ...]` 或 `[(u, v, w), ...]`，也可以是 DataFrame（列名 `source/target/weight`）。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> G = build_graph([("A", "B"), ("B", "C")]) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def shortest_path(G: Any, source: Any, target: Any = None, weight: str | None = None) -> Any:
|
||||
"""最短路径。`target=None` 时返回从 source 出发到所有点的路径。"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def min_spanning_tree(G: Any, weight: str = "weight") -> Any:
|
||||
"""最小生成树（铺管网、修路这类"连通成本最小"题目）。"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def betweenness_centrality(G: Any, normalized: bool = True) -> Any:
|
||||
"""介数中心性。
|
||||
|
||||
节点的"关键程度"有三种主流度量（度/接近/介数），
|
||||
论文里选哪一个要说明理由：**介数**衡量"必经之路"，
|
||||
**度**衡量"连接多"，**特征向量**衡量"连到重要节点"。
|
||||
原建议把四个都列了，这里只留最常用的介数，其余直接调 networkx。
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def pagerank(G: Any, alpha: float = 0.85, max_iter: int = 100) -> Any:
|
||||
"""PageRank（有向图的节点重要性）。"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def max_flow(G: Any, source: Any, sink: Any, capacity: str = "capacity") -> Any:
|
||||
"""最大流 / 最小割（网络容量、瓶颈分析）。"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def community_detection(G: Any, method: str = "louvain", random_state: int | None = None) -> Any:
|
||||
"""社区发现。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是 `{节点: 社区编号}`，`metrics={"modularity"}`。
|
||||
**必须报模块度 Q**，否则社区划分好不好无从判断。
|
||||
"""
|
||||
# TODO(v0.2)
|
||||
raise NotImplementedError
|
||||
@@ -0,0 +1,234 @@
|
||||
"""优化。
|
||||
|
||||
分层（对应原建议的 A/B/C 三级）：
|
||||
|
||||
* **B 类（薄封装）**：`minimize` `linear_program` `integer_program` —— 底层 SciPy。
|
||||
* **A 类（自己实现，S 档）**：`dynamic_programming` `knapsack_dp` `simulated_annealing`。
|
||||
* **C 类（直接用库）**：NSGA-II 等，用 pymoo，本模块只提供 Pareto 前沿的**后处理**。
|
||||
|
||||
原建议里的 `maximize` 已删：它就是 `minimize` 取负，用 `sense="max"` 参数表达。
|
||||
`value_iteration` / `policy_iteration` 也已删：与 `dynamic_programming` 同族，
|
||||
第一版同时做三个属于过度设计。
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from collections.abc import Callable
|
||||
from typing import Any
|
||||
|
||||
__all__ = [
|
||||
"dynamic_programming",
|
||||
"integer_program",
|
||||
"knapsack_dp",
|
||||
"linear_program",
|
||||
"minimize",
|
||||
"pareto_front",
|
||||
"simulated_annealing",
|
||||
]
|
||||
|
||||
|
||||
def minimize(
|
||||
objective: Callable[..., float],
|
||||
x0: Any = None,
|
||||
bounds: Any = None,
|
||||
constraints: Any = None,
|
||||
method: str = "auto",
|
||||
sense: str = "min",
|
||||
random_state: int | None = None,
|
||||
**kwargs: Any,
|
||||
) -> Any:
|
||||
"""连续优化（底层 `scipy.optimize.minimize`）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
method:
|
||||
`"auto"` 时按"有无约束 / 有无梯度"自动选：
|
||||
无约束 → BFGS；有边界 → L-BFGS-B；有约束 → SLSQP。
|
||||
sense:
|
||||
`"min"` / `"max"`。取 max 时内部对目标函数取负并还原结果，
|
||||
调用方不用自己写负号（负号写错是这类题最常见的手滑）。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`params={"x": 最优解}`，`metrics={"objective": 目标值}`。
|
||||
|
||||
Notes
|
||||
-----
|
||||
**结果必须在论文里报"是否收敛"和"初值敏感性"**：
|
||||
用一个 `x0` 跑出来的局部最优，和换 10 个随机初值跑出来的全局最优，
|
||||
结论可能完全不同。返回值里带上 `success` 与 `message` 就是为了这个。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> f = lambda v: (v[0] - 3) ** 2 + (v[1] + 1) ** 2
|
||||
>>> r = minimize(f, x0=[0, 0])
|
||||
>>> round(r.metrics["objective"], 6)
|
||||
0.0
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def linear_program(
|
||||
c: Any,
|
||||
A_ub: Any = None,
|
||||
b_ub: Any = None,
|
||||
A_eq: Any = None,
|
||||
b_eq: Any = None,
|
||||
bounds: Any = None,
|
||||
sense: str = "min",
|
||||
) -> Any:
|
||||
"""线性规划（底层 `scipy.optimize.linprog`）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
c:
|
||||
目标函数系数。`sense="max"` 时内部取负。
|
||||
bounds:
|
||||
`None` 表示默认 `x >= 0`；要无下界显式写 `(None, None)`。
|
||||
**这是最容易踩的坑**：忘记声明变量非负，解出来的方案可能物理上不可行。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = linear_program(c=[1, 2], A_ub=[[1, 1]], b_ub=[4], bounds=[(0, None)] * 2)
|
||||
>>> r.metrics["objective"]
|
||||
4.0
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def integer_program(
|
||||
c: Any,
|
||||
integrality: Any,
|
||||
A_ub: Any = None,
|
||||
b_ub: Any = None,
|
||||
A_eq: Any = None,
|
||||
b_eq: Any = None,
|
||||
bounds: Any = None,
|
||||
sense: str = "min",
|
||||
) -> Any:
|
||||
"""整数规划（底层 `scipy.optimize.milp`）。
|
||||
|
||||
`integrality`：0 = 连续，1 = 整数，2 = 0-1 变量。
|
||||
|
||||
Notes
|
||||
-----
|
||||
整数规划比 LP 慢得多。变量数到几百以上时先想想能不能用
|
||||
线性松弛 + 取整 + 邻域搜索近似，别硬等求解器。
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def dynamic_programming(
|
||||
states: Any,
|
||||
actions: Any,
|
||||
transition: Callable[..., Any],
|
||||
reward: Callable[..., float],
|
||||
horizon: int,
|
||||
discount: float = 1.0,
|
||||
) -> Any:
|
||||
"""通用有限期动态规划（A 类，自己实现）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
transition:
|
||||
`f(state, action, rng) -> next_state`。
|
||||
reward:
|
||||
`f(state, action, next_state) -> float`。
|
||||
horizon:
|
||||
期数。`discount < 1` 时按 `γ^t` 折算。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是"最优值函数表 + 最优策略表"，`extras["policy"]` 可直接查。
|
||||
|
||||
Notes
|
||||
-----
|
||||
复杂度 `O(horizon × |S| × |A|)`。**状态空间是唯一瓶颈**——
|
||||
如果 `|S| × holiday` 超过内存，要么精简状态离散化粒度，要么改成近似 DP。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> # 背包问题用 knapsack_dp 更直接；本函数用于"逐期决策"类题
|
||||
>>> r = dynamic_programming(states, actions, trans, rew, horizon=10) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L1) ★ S 档：配完整函数卡片
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def knapsack_dp(weights: Any, values: Any, capacity: float) -> Any:
|
||||
"""0-1 背包（动态规划模板，A 类，自己实现）。
|
||||
|
||||
复杂度 `O(n × capacity)`。`capacity` 在题面里是小数时先做单位缩放取整，
|
||||
否则 DP 表会大到装不下。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`params={"selected": 选中的下标}`，`metrics={"total_value", "total_weight"}`。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = knapsack_dp(weights=[1, 3, 4, 5], values=[1, 4, 5, 7], capacity=7)
|
||||
>>> r.metrics["total_value"]
|
||||
9.0
|
||||
"""
|
||||
# TODO(L1) ★ S 档
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def simulated_annealing(
|
||||
objective: Callable[..., float],
|
||||
x0: Any,
|
||||
bounds: Any = None,
|
||||
temp0: float = 1.0,
|
||||
cooling: float = 0.95,
|
||||
n_iter: int = 1000,
|
||||
sense: str = "min",
|
||||
random_state: int | None = None,
|
||||
**kwargs: Any,
|
||||
) -> Any:
|
||||
"""模拟退火（A 类，自己实现）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cooling:
|
||||
降温系数，每步 `T *= cooling`。`0.95` 配 `n_iter=1000` 是保守默认值。
|
||||
temp0:
|
||||
初始温度。**经验做法：取"随机扰动的典型劣化幅度"同量级**，
|
||||
太大等于随机游走，太小等于爬山。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`extras["history"]` 是每步的最优值轨迹，可直接喂给 `viz.plot_convergence`。
|
||||
|
||||
Notes
|
||||
-----
|
||||
* 必须记录收敛曲线并画出来——论文里"算法收敛性分析"那一节就靠它。
|
||||
* 结果对随机种子敏感，**至少跑 10 个种子**，报最好值 + 分布，
|
||||
只报一个种子的结果是不可信的。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = simulated_annealing(lambda v: sum(x * x for x in v), x0=[3.0, 4.0], random_state=0) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L1) ★ S 档：配完整函数卡片
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def pareto_front(points: Any, sense: str = "min") -> Any:
|
||||
"""从解集里筛出 Pareto 非支配解。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`params={"mask": 布尔掩码, "rank": 每个解的支配层}`。
|
||||
|
||||
Notes
|
||||
-----
|
||||
多目标问题的**后处理**入口。真正的 NSGA-II 求解器用 pymoo
|
||||
（`pip install "mcm-toolkit[moo]"`），本函数只负责"给定一组解，挑出前沿"。
|
||||
它也是画 Pareto 前沿图的必要前置步骤。
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
@@ -0,0 +1,214 @@
|
||||
"""敏感性分析与鲁棒性检验。
|
||||
|
||||
★ 建议做成招牌模块。原建议这一节方向完全正确，只是函数名偏多，
|
||||
已合并同义项、删掉 `grid_sweep`（`parameter_sweep` 的特例）。
|
||||
|
||||
论文里的"敏感性分析"一节基本就是这四个函数的组合：
|
||||
局部扰动 → OAT 单因素扫描 → 多因素参数扫描 → 全局 Sobol/Morris。
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from collections.abc import Callable
|
||||
from typing import Any
|
||||
|
||||
__all__ = [
|
||||
"local_sensitivity",
|
||||
"one_at_a_time",
|
||||
"parameter_sweep",
|
||||
"robustness_test",
|
||||
]
|
||||
|
||||
|
||||
def local_sensitivity(
|
||||
model: Callable[..., Any],
|
||||
params: Any,
|
||||
step: float = 1e-4,
|
||||
relative: bool = True,
|
||||
**kwargs: Any,
|
||||
) -> Any:
|
||||
"""局部敏感性：各参数的偏导 / 弹性。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
relative:
|
||||
`True` 时算弹性 `(∂y/∂p)·(p/y)`——**无量纲，可以跨参数比较**，
|
||||
这是局部敏感性的正确做法；`False` 时是原始偏导（有量纲，只有同量纲参数之间可比）。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是参数 → 弹性 的表，按绝对值降序。
|
||||
|
||||
Notes
|
||||
-----
|
||||
有限差分的步长是精度与截断误差的权衡：太大截断误差大，
|
||||
太小被浮点误差淹没。`1e-4` 是稳健默认值；参数本身量级差很多时用相对步长。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> s = local_sensitivity(lambda p: p["a"] ** 2 * p["b"], {"a": 2.0, "b": 3.0}) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L1) ★ S 档
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def one_at_a_time(
|
||||
model: Callable[..., Any],
|
||||
params: Any,
|
||||
parameter: str | None = None,
|
||||
values: Any = None,
|
||||
n_points: int = 20,
|
||||
span: float = 0.5,
|
||||
) -> Any:
|
||||
"""OAT（一次只动一个参数）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
parameter:
|
||||
`None` 时对所有参数依次扫描。
|
||||
span:
|
||||
相对基准值的扫描范围，`0.5` 表示 ±50%。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是长表 `(parameter, value, output)`，配 `viz.plot` 直接出蛛网图/敏感性曲线。
|
||||
|
||||
Notes
|
||||
-----
|
||||
OAT 的**已知缺陷**：它看不见参数之间的交互作用。
|
||||
如果两个参数是相乘关系（很常见），OAT 会严重低估整体敏感性——
|
||||
论文里必须说明这一点，并补一个全局敏感性作为交叉验证。
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def parameter_sweep(
|
||||
model: Callable[..., Any],
|
||||
ranges: dict[str, Any],
|
||||
n_points: int = 20,
|
||||
method: str = "grid",
|
||||
n_samples: int = 1000,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""多参数扫描（合并了原建议的 `grid_sweep` 与 `random_parameter_sweep`）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
ranges:
|
||||
`{"alpha": (0.1, 0.9), "beta": (0, 5)}`。
|
||||
method:
|
||||
`"grid"`：笛卡尔积，点数 = `n_points^d`，**维度一高就爆**（d=5, n=20 是 320 万次）；
|
||||
`"lhs"`：拉丁超立方，维度高时用这个；
|
||||
`"random"`：纯随机。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是 `(n_points^d, d + 1)` 的结果表，
|
||||
`params["best"]` 是让目标最优的参数组合。
|
||||
|
||||
Notes
|
||||
-----
|
||||
`d >= 4` 时默认改走 `"lhs"`，否则 `n_points` 会被指数放大到跑不完。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = parameter_sweep(lambda a, b: -(a - 0.3) ** 2 - (b - 2) ** 2,
|
||||
... {"a": (0, 1), "b": (0, 4)}, n_points=11) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L0)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def robustness_test(
|
||||
model: Callable[..., Any],
|
||||
scenarios: Any = None,
|
||||
params: Any = None,
|
||||
perturbation: float = 0.1,
|
||||
n_samples: int = 1000,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""鲁棒性检验（合并了原建议的 `perturbation_test`）。
|
||||
|
||||
在参数上施加 ±`perturbation` 的随机扰动，看结论是否翻盘。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`metrics` 含 `"rank_stability"`（原结论在前 X% 的比例）、
|
||||
`"conclusion_flip_rate"`（结论反转频率）。**这两个数字才是审稿人想看的。**
|
||||
|
||||
Notes
|
||||
-----
|
||||
鲁棒性检验和敏感性分析不是一回事：
|
||||
敏感性问"参数变了输出变多少"，鲁棒性问"参数变了**结论**还成立吗"。
|
||||
论文里两节都要有，别混着写。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = robustness_test(model, params={"alpha": 0.5}, random_state=0) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def sobol_analyze(
|
||||
model: Callable[..., Any],
|
||||
bounds: Any,
|
||||
n_samples: int = 1024,
|
||||
calc_second_order: bool = False,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""Sobol 全局敏感性（底层 SALib，自己不要重写）。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
bounds:
|
||||
`{"alpha": (0, 1), "beta": (0, 5)}` 形式的参数上下界。
|
||||
n_samples:
|
||||
SALib 的基样本数 N；**实际模型调用次数是 `N·(2d + 2)`**
|
||||
（二阶时更多），d=5、N=1024 就是 12288 次调用。先估一下模型单次耗时。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 含 `S1`（一阶）、`ST`（总效应）、`S2`（二阶，可选）三列。
|
||||
|
||||
Notes
|
||||
-----
|
||||
* `ST - S1` 的差距就是**交互作用**的大小——这正是 OAT 看不见的东西。
|
||||
两者差得多，说明不能只看单因素。
|
||||
* 需要 `pip install "mcm-toolkit[sensitivity]"`；缺 SALib 时本函数给出
|
||||
可读的安装提示，而不是裸抛 `ModuleNotFoundError`。
|
||||
* 费用高：**先做 OAT 缩小参数范围，再对少数关键参数做 Sobol**，
|
||||
不要一上来就对全部参数做全局敏感性。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = sobol_analyze(model, bounds={"a": (0, 1), "b": (0, 5)}) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L1) ★ S 档：配完整函数卡片
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def morris_analyze(
|
||||
model: Callable[..., Any],
|
||||
bounds: Any,
|
||||
n_trajectories: int = 10,
|
||||
n_levels: int = 4,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""Morris 筛选法（底层 SALib）。
|
||||
|
||||
调用次数只与轨迹数成正比（`(d + 1) · n_trajectories`），比 Sobol 便宜一个量级，
|
||||
适合**参数多但只需要找出"哪几个重要"**的筛选阶段。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 含 `mu`、`mu_star`（推荐看的）、`sigma` 三列。
|
||||
`sigma` 大表示该参数的效应是非线性的或与其他参数有交互。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = morris_analyze(model, bounds={"a": (0, 1), "b": (0, 5)}) # doctest: +SKIP
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
@@ -0,0 +1,218 @@
|
||||
"""随机模拟：Monte Carlo / 马尔可夫 / 随机游走 / 排队。
|
||||
|
||||
★ 本模块的 `monte_carlo` 与 `monte_carlo_integral` 是 S 档，建议自己实现。
|
||||
原建议里的 `brownian_motion` `geometric_brownian_motion` 已删（niche），
|
||||
`simulate_queue` 也已删（等于自己写一个离散事件仿真引擎，成本远高于收益）。
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
from collections.abc import Callable
|
||||
from typing import Any
|
||||
|
||||
__all__ = [
|
||||
"markov_simulate",
|
||||
"markov_steady_state",
|
||||
"mm1_metrics",
|
||||
"monte_carlo",
|
||||
"monte_carlo_integral",
|
||||
"monte_carlo_probability",
|
||||
"random_walk",
|
||||
]
|
||||
|
||||
|
||||
def monte_carlo(
|
||||
sampler: Callable[..., Any],
|
||||
simulator: Callable[..., Any],
|
||||
n: int = 10000,
|
||||
random_state: int | None = None,
|
||||
**kwargs: Any,
|
||||
) -> Any:
|
||||
"""通用蒙特卡洛：采样 → 仿真 → 汇总。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
sampler:
|
||||
`f(rng, **kwargs) -> 一组参数`，例如从先验分布里抽参数。
|
||||
simulator:
|
||||
`f(params) -> float`，例如"给定参数下的年收益"。
|
||||
n:
|
||||
模拟次数。**必须向量化**：`simulator` 内部不要写 Python 循环，
|
||||
否则 10^4 次能跑几分钟，10^6 次直接放弃。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是 `n` 个结果组成的数组，
|
||||
`metrics={"mean", "std", "median", "p05", "p95", "ci_low", "ci_high"}`。
|
||||
|
||||
Notes
|
||||
-----
|
||||
比赛里 Monte Carlo 的收敛速度是 `O(1/√n)`——
|
||||
想少一半误差要 4 倍样本。所以论文里要给"区间随 n 收敛"的图，
|
||||
证明 `n` 已经够大（这也是 `viz.plot_convergence` 的典型用法）。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> r = monte_carlo(lambda rng, **k: rng.normal(0, 1), np.mean, n=1000, random_state=0)
|
||||
"""
|
||||
# TODO(L1) ★ S 档
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def monte_carlo_probability(
|
||||
event: Callable[..., bool],
|
||||
sampler: Callable[..., Any],
|
||||
n: int = 10000,
|
||||
confidence: float = 0.95,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""估计事件发生概率，并给出区间。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`metrics={"probability", "ci_low", "ci_high", "std_error"}`。
|
||||
|
||||
Notes
|
||||
-----
|
||||
概率估计的相对误差 ≈ `1/√(n·p)`。**小概率事件（p < 1e-3）用朴素 MC 是不可行的**，
|
||||
需要重要性采样或方差缩减——遇到这类题要提前意识到，别硬跑。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> r = monte_carlo_probability(lambda v: v > 0, lambda rng, **k: rng.normal(0, 1), n=2000, random_state=0)
|
||||
>>> abs(r.metrics["probability"] - 0.5) < 0.1
|
||||
True
|
||||
"""
|
||||
# TODO(L1) ★ S 档
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def monte_carlo_integral(
|
||||
func: Callable[..., Any],
|
||||
bounds: Any,
|
||||
n: int = 10000,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""用蒙特卡洛求多重积分。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
bounds:
|
||||
`[(lo, hi), ...]`，维度即积分重数。
|
||||
|
||||
Notes
|
||||
-----
|
||||
高维积分是 MC 相对网格法的**唯一绝对优势**：网格法误差 `O(h^{-d})` 随维度爆炸，
|
||||
MC 始终是 `O(1/√n)`。三维以上优先考虑 MC。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> r = monte_carlo_integral(lambda X: X[:, 0] ** 2, bounds=[(0, 1)], n=200000, random_state=0)
|
||||
>>> abs(r.data - 1 / 3) < 0.01
|
||||
True
|
||||
"""
|
||||
# TODO(L1) ★ S 档：配完整函数卡片
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def random_walk(
|
||||
n_steps: int,
|
||||
n_paths: int = 1,
|
||||
step_sampler: Callable[..., Any] | None = None,
|
||||
start: float = 0.0,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""随机游走轨迹。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
step_sampler:
|
||||
`f(rng, size) -> 增量`，默认为 ±1 等概率（对称随机游走）。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 形状 `(n_paths, n_steps + 1)`。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = random_walk(n_steps=100, n_paths=5, random_state=0)
|
||||
>>> r.data.shape
|
||||
(5, 101)
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def markov_simulate(
|
||||
transition: Any,
|
||||
initial: Any,
|
||||
steps: int = 100,
|
||||
n_paths: int = 1,
|
||||
random_state: int | None = None,
|
||||
) -> Any:
|
||||
"""马尔可夫链模拟，返回状态序列。"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e-10) -> Any:
|
||||
"""平稳分布。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
transition:
|
||||
行随机矩阵（每行和为 1）。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`data` 是平稳分布向量，`metrics={"converged", "iterations"}`。
|
||||
|
||||
Notes
|
||||
-----
|
||||
用迭代法而不是特征分解：迭代法能顺便告诉你"收敛快不快"，
|
||||
而收敛速度本身就是论文里值得写一句的性质。
|
||||
**不可约 + 非周期**才有唯一平稳分布，周期链会不收敛，要检查。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> import numpy as np
|
||||
>>> P = np.array([[0.9, 0.1], [0.2, 0.8]])
|
||||
>>> r = markov_steady_state(P)
|
||||
>>> bool(np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6))
|
||||
True
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
|
||||
|
||||
def mm1_metrics(arrival_rate: float, service_rate: float) -> Any:
|
||||
"""M/M/1 排队系统的闭式指标。
|
||||
|
||||
Parameters
|
||||
----------
|
||||
arrival_rate, service_rate:
|
||||
到达率 λ 与服务率 μ。`λ >= μ` 时系统不稳定，队列长度无界——
|
||||
这种情况下函数返回 `success=False` 并说明原因，而不是返回 `inf`。
|
||||
|
||||
Returns
|
||||
-------
|
||||
`Result`：`metrics={"rho", "L", "Lq", "W", "Wq"}`
|
||||
（系统利用率、系统内平均顾客数、队列长度、平均逗留时间、平均等待时间）。
|
||||
|
||||
Notes
|
||||
-----
|
||||
M/M/1 是唯一值得手写的排队模型：公式是闭式的，10 行搞定。
|
||||
多服务员（M/M/c）和一般服务时间分布（M/G/1）在美赛出现频率极低，
|
||||
真遇到时现场用仿真做，不做通用封装。
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> r = mm1_metrics(arrival_rate=2, service_rate=3)
|
||||
>>> round(r.metrics["rho"], 4)
|
||||
0.6667
|
||||
"""
|
||||
# TODO(L1)
|
||||
raise NotImplementedError
|
||||
Reference in new issue
Block a user