From dd9bd638b8f74449740ff9635847b028648cc4d7 Mon Sep 17 00:00:00 2001 From: guzhujushi Date: Mon, 5 Oct 2026 23:11:04 +0800 Subject: [PATCH] =?UTF-8?q?feat(solve):=20=E5=AE=A1=E6=9F=A5=20markov=5Fsi?= =?UTF-8?q?mulate=EF=BC=9A=E5=B7=A5=E7=A8=8B=E4=BE=A7=E6=94=B6=E5=8F=A3=20?= =?UTF-8?q?+=209=20=E6=9D=A1=E6=B5=8B=E8=AF=95?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 用户实现的第一版 markov_simulate 算法本体经探针实测正确（终点频率 [0.657, 0.343] 对照平稳分布 [2/3, 1/3]，同种子可复现），本次只动工程侧，未改任何算法行： - ruff 8 条全清（I001 导入排序 + 5 处 W293 / 2 处 W291） - docstring 从一行补全为 Parameters / Returns / Raises / Notes / Examples - 文案：马尔科夫→马尔可夫（3 处）、initail→initial、半角逗号→全角、paths[:,0]→paths[:, 0] - arr, kind = as_array(...) 的 kind 全程未用 → arr, _ - tests/test_solve.py 新增 9 条 markov_simulate 测试（此前 0 条） - STATUS.md 新增决策 56（initial 语义待用户拍板，A/B/C 三选项）与 5 条待办 - CHANGELOG 与 docs/dev/personal-pitfalls.md 易错点 27（入口分流后两条分支口径矛盾） 校验：pytest 188 passed / 31 skipped（219 收集）；ruff 0；pyright 0 --- CHANGELOG.md | 21 +++++ STATUS.md | 78 ++++++++++++++---- docs/dev/personal-pitfalls.md | 15 ++++ src/mcm_toolkit/solve/simulate.py | 130 ++++++++++++++++++++++++++++-- tests/test_solve.py | 102 +++++++++++++++++++++++ 5 files changed, 323 insertions(+), 23 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8705822..fb286b6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -508,6 +508,27 @@ `solve` 段 42 → **43** 条，全量 **210 收集 / 179 passed / 31 skipped / 0 xfail**（ruff 0 / pyright 0）。 - `STATUS.md`：决策 55 改为"已选 B 并落地"；`random_walk` 的 7 条待办全部办结；验证计数更新为 179 / 31。 +### Added +- `solve.markov_simulate` 落地（**L1，S 档**，实现由用户本人写）：`arr[current_states]` 取当前状态那一行 → + `cumsum` 得累积分布 → 用 `[0, 1)` 均匀数"落格"（`(cum_probs < r[:, None]).sum(axis=1)`）→ 写进下一列； + `data` 形状 `(n_paths, steps + 1)`、值是整数状态编号。Codex 审查结论：**算法本体正确** + （探针实测 200 步后终点频率 `[0.657, 0.343]` 对照平稳分布 `[2/3, 1/3]`、同种子可复现）。 +- `tests/test_solve.py`：新增 9 条 `markov_simulate` 测试（此前 0 条）—— 形状 / 首列 / 状态取值 / 整数 dtype、 + 单位矩阵吸收态自锁、可复现（同种子相同、异种子不同）、终点频率收敛到 `[2/3, 1/3]`、metrics = 终点状态、 + 非法转移矩阵（非方阵 / 行和不为 1）、非法 `steps` / `n_paths`（0 / float / bool）、`np.int64` 可收、docstring Example 真能跑。 + +### Changed +- `markov_simulate` 的 docstring 从一行补到完整的 `Parameters` / `Returns` / `Raises` / `Notes` / `Examples`。 +- `markov_simulate` 的工程侧清理：`kind` 全程未用 → `arr, _`；`马尔科夫` → `马尔可夫`（3 处）、 + `initail` → `initial`、半角逗号 → 全角；`paths[:,0]` → `paths[:, 0]`；ruff 8 条全清（`I001` + `W291` / `W293`）。 +- `STATUS.md`：新增**决策 56**（`initial` 的语义待拍板，含 A / B / C 三个选项）与 5 条待办； + 验证计数更新为 **188 passed / 31 skipped（219 收集）**。 + +### Notes +- **`markov_simulate` 的 5 条口径留给用户**（决策 49）：`initial` 语义（标量按路径数 / 数组按状态数，两条口径互相矛盾）、 + 转移矩阵未查非负（行和恰为 1.0 的 `[[1.5, -0.5], [0, 1]]` 会静默给出全 0）、`paths[:, 0]` 赋值在 `try` 之外 + （长度不匹配时裸抛 `ValueError`）、状态越界没在入口拦、`std_terminal` 对分类状态的意义。 + ## [0.0.1] - 2026-10-01 ### Added diff --git a/STATUS.md b/STATUS.md index 8c24448..d460d0c 100644 --- a/STATUS.md +++ b/STATUS.md @@ -21,9 +21,12 @@ 当前 **L0 11/12**：`Result`、`get_rng`/`set_seed`、`load_table`/`save_table`、`clean_missing`、 `detect_outliers`、`normalize`/`standardize`、`describe`、`correlation_matrix`、`bootstrap_ci`、 `monte_carlo`（**用户写**，2026-10-04）；另**提前落地**了 L1 的 `confidence_interval` / `t_test` / -`monte_carlo_probability`、`monte_carlo_integral`、`random_walk`（**用户写**，2026-10-05）。剩余 1 项：`save_figure`（Codex 写）。 -`monte_carlo` 的口径已拍板（决策 48）；`monte_carlo_probability` / `random_walk` 的待办**均已办结**， -剩下 `ci_low`/`ci_high` 的语义、`monte_carlo_integral` 的 1 条口径（第一道 `isinstance` 分支已冗余），见"下一步"。 +`monte_carlo_probability`、`monte_carlo_integral`、`random_walk`、`markov_simulate` +（**均为用户写**，2026-10-05）。剩余 1 项：`save_figure`（Codex 写）。 +`monte_carlo` 的口径已拍板（决策 48）；`monte_carlo_probability` / `random_walk` 的待办**均已办结**； +`markov_simulate` 本轮审查通过（算法本体正确），留 **5 条入口口径待你拍板**（见"下一步"，决策 56 的 3 个选项）。 +另 3 个小口径仍悬着：`ci_low`/`ci_high` 的语义、`monte_carlo_integral` 的第一道 `isinstance` 分支可塌成两行、 +`random_state` 注解与 `std` 的 `ddof` 是否全项目对齐。 ## 任务清单 @@ -317,13 +320,28 @@ - **第三轮（2026-10-05）**：`start` 校验恢复、注解放宽为 `int | np.integer`；最后一条"方阵歧义"用户选**方案 B**， `step_sampler` 契约改为"整数 `size=n_paths * n_steps` + 必须返回一维 + 内部行主序 `reshape`"（决策 55，Codex 代改） +- [x] `solve/simulate.py` 的 `markov_simulate`（**L1，S 档，算法由用户本人写**，2026-10-05） + - 核心：`arr[current_states]` 取当前状态那一行 → `cumsum` 得累积分布 → 用 `[0, 1)` 均匀数"落格" + （`(cum_probs < r[:, None]).sum(axis=1)`）→ 写进下一列。`data` 形状 `(n_paths, steps + 1)`，值是整数状态编号 + - **Codex 审查结论：算法本体正确**（探针实测 `(2000, 201)`、状态 ∈ {{0,1}}、200 步后终点频率 + `[0.657, 0.343]` 对照平稳分布 `[2/3, 1/3]`、同种子可复现） + - **Codex 只碰工程侧**：① ruff 8 条全清（`I001` 导入排序 + 5 处 `W293` / 2 处 `W291` 空白）； + ② `arr, kind = as_array(...)` 的 `kind` 全程未用 → 改 `arr, _`； + ③ docstring 从一行补到完整的 `Parameters` / `Returns` / `Raises` / `Notes` / `Examples`； + ④ 文案：`马尔科夫` → `马尔可夫`（3 处）、`initail` → `initial`、半角逗号 → 全角；`paths[:,0]` → `paths[:, 0]` + - **测试：新增 9 条**（此前 `markov_simulate` 一条都没有）—— 形状 / 首列 / 状态取值 / 整数 dtype、 + 单位矩阵吸收态自锁、可复现（同种子相同、异种子不同）、终点频率收敛到 `[2/3, 1/3]`、 + metrics = 终点状态、非法转移矩阵（非方阵 / 行和不为 1）、非法 `steps` / `n_paths`（0 / float / bool）、 + `np.int64` 可收、docstring Example 真能跑 + - ⏳ **留给你的 5 条**见"下一步 → markov_simulate 待你拍板"（全是入口口径与少考虑情况，决策 49） + - [ ] 剩余 1 项：`save_figure`（A 档 = 工程，**Codex 直接写**） ### 验证结果 | 检查 | venv | conda `mcm-toolkit` | |---|---|---| -| `pytest` | 见下（单列） | **179 passed, 31 skipped（0 xfail）** | +| `pytest` | 见下（单列） | **188 passed, 31 skipped（0 xfail）** | | `import mcm_toolkit` | 0.0.1 / 16 项 | 0.0.1 / 16 项 | | `ruff check` | All checks passed | **All checks passed**（0 error） | | `scripts/doctor.ps1` | 跑通 | — | @@ -333,14 +351,15 @@ （修配置之前是 17 errors，全是"无法解析导入 numpy / pandas"，见决策 29 与 `docs/04-environment.md` 坑 7） - base 亦已复核：numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6（与两套环境一致） -- **210 条用例：179 passed / 31 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 - `test_core.py` 23、`test_solve.py` 43（本轮 +1：采样器契约"整数 `size` + 行主序"）、`test_viz.py` 8、`test_decision.py` 7、 +- **219 条用例：188 passed / 31 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 + `test_core.py` 23、`test_solve.py` 51（本轮 +9：`markov_simulate` 从 0 条起步）、`test_viz.py` 8、`test_decision.py` 7、 `test_stats.py` 53、`test_models.py` 5 - **`xfail` 再次清零**：`correlation_matrix` 的契约冲突已按方案 A 解决（决策 41），那条用例转正 - 上一轮那 7 条 xfail 已全部转正，决策 30 那套信号（XPASS → 摘 mark）跑完过一个完整闭环 - `examples/01_clean_and_describe.py` 实测跑通（exit 0）；`00_quickstart.py` 仍等 `stats` / `viz` -- 31 条 skip 是等实现的算法（`train_test_split`、`random_walk`、`markov_*`、`anova`、`fit_distribution` 等）； - `monte_carlo_probability` / `monte_carlo_integral` 已落地，分别摘掉 2 条与 1 条 skip； +- 31 条 skip 是等实现的算法（`train_test_split`、`markov_steady_state`、`mm1_metrics`、`anova`、`fit_distribution`、 + `knapsack_dp`、`linear_program`、`simulated_annealing` 等）； + `monte_carlo_probability` / `monte_carlo_integral` 已落地，分别摘掉 2 条与 1 条 skip；`markov_simulate` 新增的 9 条全部实跑； `TestValidate` 的 4 条仍挂着 —— `check_*` 还是 `NotImplementedError`，但它们实现时应当直接复用 `as_array` - **本轮回归（已修）**：`transform.py` 上一版有 3 处真 bug —— 修前 **7 failed / 75 passed**， 且 pyright 在同一处报 `reportArgumentType`。**不是静态误报**：pyright CLI 复现 + 测试真红。 @@ -499,22 +518,47 @@ docstring 与 `test_describe_reports_expected_keys` 早就写死了这个方向 歧义原样回来。**代价（已写进 docstring）**：`step_sampler` 不能再返回二维，步与路径的对应关系是**行主序**。 探针实测：默认采样器 `(5, 101)` 正常；后四种（转置二维 / 方阵二维 / 长度 +1 / `(N, 1)`）一律判失败。 +| 56 | `markov_simulate` 的 `initial` 语义**待用户拍板**（Codex 只报选项，未代改，决策 49） | +现状用 `isinstance(initial, int)` 区分"标量状态 / 状态数组"，两条分支的口径互相打架：标量分支按**路径数**铺满 +（`np.full(n_paths, initial)`），数组分支却按**状态数**校验长度（`len(initial) != s` 直接报错）—— +于是数组分支只有 `n_paths == s` 时才可能跑通；而 `paths[:, 0] = current_states` 又写在 `try` **外面**， +长度不匹配时**裸抛** `ValueError`（本该 `success=False`，同易错点 26 一族）。另有三个口子： +`np.int64(0)` 掉进数组分支、被 `as_array` 拒（与 `steps`/`n_paths` 收 `np.integer` 不一致，易错点 1 复发）； +`True` 因 `isinstance(True, int)` 为真被当状态 1 **静默放行**；越界状态（如 `initial=5`、`s=2`）只在 +第一轮索引时才炸、被 `except` 吞成 `success=False`。可选口径见下一步的 A / B / C。 + ## 下一步 -### ⭐ 当前进度与待办（2026-10-05 第四轮） +### ⭐ 当前进度与待办（2026-10-05 第五轮） -**本轮：`random_walk` 第二版（用户按 4 条待办改定）已审查通过**。Codex 又抓到一处**真 bug**： -`call_params` 定义在 `try` 内、却在 `except` 里被引用，**所有失败路径**都抛 `UnboundLocalError`（易错点 26，见决策 54）； -另补 docstring、+3 条测试（共 9 条）。 +**本轮：`markov_simulate`（你写的第一版）已审查通过** —— 算法本体正确，没动你任何算法行； +Codex 只清工程侧（ruff 8 条、docstring、3 处 `马尔科夫`、`initail` 拼写、半角逗号、`kind` 未用）， +并补 **9 条测试**（此前 0 条）。三连：`pytest` **188 passed / 31 skipped**、`ruff` 0、`pyright` 0。 **接下来三件事** -1. `random_walk` 的 7 条待办**全部办结**（4 条 + 本轮 3 条），下面一节留档。 -2. `ci_low` / `ci_high` 的语义确认一句（也是你的活），确认后 Codex 补 `docs/algorithms/monte_carlo.md` 卡片。 +1. **`markov_simulate` 的 5 条待你拍板**（下节，全是入口口径，决策 49/56）。 +2. `ci_low` / `ci_high` 的语义确认一句（A / B / C 三个选项在下面挂着），确认后 Codex 补 `docs/algorithms/monte_carlo.md` 卡片。 3. L0 最后一项 `save_figure`（A 档 = 工程，`viz/export.py`）—— **Codex 直接实现**，你只读。 -**你下一个要写的（S 档，仍在 `solve/simulate.py`）**：`markov_simulate` → `markov_steady_state`。 -两者是同一件事的两面：前者按转移矩阵**逐代抽状态**（`rng.choice` 的 `p=` 版），后者用**迭代**求平稳分布。 -`markov_steady_state` 的 docstring 已写明口径：用迭代而不用特征分解，因为"收敛快不快"本身就是论文里要写的性质。 +### ⏳ markov_simulate 待你拍板（5 条，Codex 只报不改） + +**① `initial` 的语义（最大的一条，三个选项）** +现在标量分支按**路径数**铺满、数组分支按**状态数**校验长度，两条口径互相矛盾（详见决策 56）。 +- **A**：只收标量状态，另加 `initial_dist` 参数显式传"初始分布"（分布再抽一次，语义最干净，签名变长） +- **B**：只收长度 `n_paths` 的数组（标量走 `np.full` 铺满），`data` 首列就是给定的初始状态序列（与 `paths[:, 0]` 赋值天然一致） +- **C**：沿用现状，但把类型收成 `int | np.integer` + 先判 `bool`，把 `paths[:, 0] = ...` 挪进 `try`，docstring 写清"数组按路径数" + +**② 转移矩阵没查非负**：`[[1.5, -0.5], [0, 1]]` 行和恰为 1.0 → 通过校验，但 `cumsum` 非单调 → +落格结果**全是 0**（静默错）。建议入口加一条 `np.any(arr < 0)` → `ValueError`。 +**③ `paths[:, 0]` 赋值在 `try` 之外**：长度不匹配时裸抛 `ValueError`，与"出错统一 `success=False`"（决策 48）不一致。 +**④ 状态越界没在入口拦**：`initial=5`（`s=2`）第一轮才 `IndexError`，被吞成 `success=False`； +建议入口加 `0 <= initial < s`。 +**⑤ `std_terminal` 对**分类**状态的意义**：状态编号是标称量，均值和标准差只对"有序状态"（如排队人数）成立 —— +要不要把 `metrics` 换成状态频率（`terminal_freq`）？你定。 + +**你下一个要写的（S 档，仍在 `solve/simulate.py`）**：`markov_steady_state`。 +与 `markov_simulate` 是同一件事的两面：前者按转移矩阵逐代抽状态，后者用**迭代**求平稳分布。 +docstring 已写明口径：用迭代而不用特征分解，因为"收敛快不快"本身就是论文里要写的性质。 ### ✅ random_walk 的 4 条待办（2026-10-05 已办结） diff --git a/docs/dev/personal-pitfalls.md b/docs/dev/personal-pitfalls.md index daa8832..94198fa 100644 --- a/docs/dev/personal-pitfalls.md +++ b/docs/dev/personal-pitfalls.md @@ -453,3 +453,18 @@ ② 每写一个 `try`，必须配一个**注定失败**的用例（`assert r.success is False`），专门喂 `except` 分支； ③ 包装层的失败结果也要带 `params`，出问题才追得回去。 - 首次：2026-10-05｜复发：0 + +### 27. 入口按类型分成两条路，两条路的校验口径互相矛盾 + +- 现象：`markov_simulate` 的 `initial` 用 `isinstance(initial, int)` 分流。标量那条按**路径数**铺满 + （`np.full(n_paths, initial)`），数组那条却按**状态数**校验（`len(initial) != s` 报错）， + 而真正写进去的是 `paths[:, 0] = current_states`，长度要等于**路径数**。于是"数组 + `n_paths != s`" + 这条组合永远跑不通；`n_paths == s` 时又能跑通 —— **同一个入口，两个维度各管一半**。 + 实测 `initial=[0.5, 0.5], n_paths=5` 直接在 `try` 外抛 `ValueError`（本该 `success=False`）。 +- 为什么会错：写校验时是"顺着刚写完的那一行"补的 —— 标量分支当年是为了铺满路径，数组分支是照着"转移矩阵的状态数" + 顺手写的，两边各自看着都对。**没有回到参数本身的语义上问一句"它到底是按谁计数"**。 + 同族：`np.int64(0)` 掉进数组分支（`isinstance` 只写了 `int`，易错点 1）、`True` 被当状态 1 静默放行。 +- 以后怎么防：① 一个参数有多条类型分支时，**把每条分支的"计数维度"写在一张纸上**（按路径数？按状态数？按列数？）， + 对不上就说明口径没定清，**先定语义再写代码**；② 分流用的 `isinstance` 必须把 `np.integer` / `bool` 一起想到（易错点 1）； + ③ 任何"写进结果数组"的赋值，问一遍"它的长度约束和上面校验用的是同一个下标维度吗"。 +- 首次：2026-10-05｜复发：0 diff --git a/src/mcm_toolkit/solve/simulate.py b/src/mcm_toolkit/solve/simulate.py index 1d98feb..08e2f52 100644 --- a/src/mcm_toolkit/solve/simulate.py +++ b/src/mcm_toolkit/solve/simulate.py @@ -13,7 +13,7 @@ from typing import Any import numpy as np from scipy.stats import norm -from mcm_toolkit.core import Result, get_rng +from mcm_toolkit.core import Result, as_array, get_rng __all__ = [ "markov_simulate", @@ -463,13 +463,131 @@ def random_walk( def markov_simulate( transition: Any, initial: Any, - steps: int = 100, - n_paths: int = 1, + steps: int | np.integer = 100, + n_paths: int | np.integer = 1, random_state: int | None = None, ) -> Any: - """马尔可夫链模拟，返回状态序列。""" - # TODO(L1) - raise NotImplementedError + """马尔可夫链模拟：从 `initial` 出发，按 `transition` 逐行走 `steps` 步。 + + Parameters + ---------- + transition: + 行随机矩阵（每行和为 1）的二维方阵，状态数就是它的阶数 `s`。 + initial: + 初始状态。**标量整数**时所有路径都从同一状态出发； + 否则按一维数组处理，长度必须等于状态数 `s`，再写进每条路径的第 0 列 + —— 只有 `n_paths == s` 时两者才对得上，见 Notes。 + steps: + 每条路径的步数，正整数。`data` 的列数是 `steps + 1`（多出的第 0 列是初始状态）。 + n_paths: + 路径条数，正整数，默认 1 条。 + random_state: + 随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。 + + Returns + ------- + `Result`：`data` 是形状 `(n_paths, steps + 1)` 的整数状态序列，第 0 列是初始状态； + `metrics={"mean_terminal", "std_terminal"}`（终点状态的均值与标准差，`std` 用 `ddof=0`）。 + + Raises + ------ + TypeError: + `transition` 不是二维数组；`initial` 数组长度不等于状态数； + `steps` / `n_paths` 不是整数。 + ValueError: + `transition` 不是方阵，或某一行和不等于 1.0；`steps` / `n_paths` 不是正整数。 + + Notes + ----- + 抽样用逆变换法：`cumsum` 得到累积分布，再拿 `[0, 1)` 的均匀数去"落格"。 + 一次抽完 `n_paths` 个状态，路径之间相互独立；步与步之间只能串行，因为下一步依赖当前状态。 + `cum_probs[:, -1] = 1.0` 是为了压掉浮点累加误差，避免随机数落在 `[1 - 1e-16, 1)` 时索引越界。 + 终点分布随步数收敛到平稳分布：这既是本函数的自检手段，也是论文里值得写一句的性质。 + + 遗留（决策 56 待定）：`initial` 现在靠 `isinstance(initial, int)` 区分"标量状态 / 状态数组"， + `np.int64`、`bool` 会掉进数组分支，"按状态数校验长度"与"按路径数写入"也互相矛盾。 + 口径定下来后，上面 `initial` 这一节与本段随之改写。 + + Examples + -------- + >>> import numpy as np + >>> P = np.array([[0.9, 0.1], [0.2, 0.8]]) + >>> r = markov_simulate(P, initial=0, steps=50, n_paths=4, random_state=0) + >>> r.data.shape + (4, 51) + >>> bool(np.all(r.data[:, 0] == 0)) + True + """ + arr, _ = as_array(transition) + if arr.ndim != 2 or arr.shape[0] != arr.shape[1]: + raise ValueError(f"transition 必须是二维方阵，你传入的形状是 {arr.shape}") + if not np.allclose(arr.sum(axis=1), 1.0): + raise ValueError("transition 每一行的和必须是 1.0") + s = arr.shape[0] + + if isinstance(n_paths, bool) or not isinstance(n_paths, (int, np.integer)): + raise TypeError(f"n_paths 必须是整数，你传入了 {type(n_paths).__name__}") + if n_paths <= 0: + raise ValueError(f"n_paths 必须是正整数，你传入了 {n_paths}") + if isinstance(steps, bool) or not isinstance(steps, (int, np.integer)): + raise TypeError(f"steps 必须是整数，你传入了 {type(steps).__name__}") + if steps <= 0: + raise ValueError(f"steps 必须是正整数，你传入了 {steps}") + + rng = get_rng(random_state) + + if isinstance(initial, int): + current_states = np.full(n_paths, initial) + else: + current_states, _ = as_array(initial, ndim=1) + if len(current_states) != s: + raise TypeError(f"initial 的长度必须是 {s}，而你传入的是 {len(current_states)}") + + paths = np.zeros((n_paths, steps + 1), dtype=int) + paths[:, 0] = current_states + + call_params = { + "n_steps": steps, + "n_paths": n_paths, + "initial": initial, + "random_state": random_state, + } + + try: + for t in range(steps): + probs = arr[current_states] + + cum_probs = np.cumsum(probs, axis=1) + cum_probs[:, -1] = 1.0 + + r = rng.random(n_paths) + next_states = (cum_probs < r[:, None]).sum(axis=1) + + paths[:, t + 1] = next_states + current_states = next_states + + terminal_positions = paths[:, -1] + metrics = { + "mean_terminal": float(np.mean(terminal_positions)), + "std_terminal": float(np.std(terminal_positions)), + } + return Result( + name="markov_simulate", + params=call_params, + metrics=metrics, + data=paths, + success=True, + message=f"马尔可夫链模拟成功，生成了 {n_paths} 条路径，每条路径 {steps} 步", + ) + except Exception as e: + return Result( + name="markov_simulate", + params=call_params, + metrics={}, + data=None, + success=False, + message=f"马尔可夫链模拟失败：{type(e).__name__} - {str(e)}", + ) def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e-10) -> Any: diff --git a/tests/test_solve.py b/tests/test_solve.py index 7e68f82..5549176 100644 --- a/tests/test_solve.py +++ b/tests/test_solve.py @@ -407,6 +407,108 @@ def test_random_walk_docstring_example_runs(): assert r.data.shape == (5, 101) assert np.all(np.abs(np.diff(r.data, axis=1)) == 1) +def test_markov_simulate_shape_start_and_states(): + """`data` 形状 (n_paths, steps + 1)，第 0 列是初始状态，取值只能是合法状态编号。""" + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_simulate(P, initial=0, steps=100, n_paths=5, random_state=0) + assert r.success + assert r.data.shape == (5, 101) + assert np.all(r.data[:, 0] == 0) + assert set(np.unique(r.data)) <= {0, 1} + assert np.issubdtype(r.data.dtype, np.integer) + + +def test_markov_simulate_absorbing_states_never_move(): + """单位转移矩阵：状态自锁，整条路径恒等于初始状态（确定性自检）。""" + from mcm_toolkit.solve import markov_simulate + + r = markov_simulate(np.eye(2), initial=1, steps=30, n_paths=4, random_state=0) + assert r.success + assert np.all(r.data == 1) + assert r.metrics["mean_terminal"] == 1.0 + assert r.metrics["std_terminal"] == 0.0 + + +def test_markov_simulate_is_reproducible(): + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.7, 0.3], [0.4, 0.6]]) + a = markov_simulate(P, initial=0, steps=40, n_paths=6, random_state=11) + b = markov_simulate(P, initial=0, steps=40, n_paths=6, random_state=11) + c = markov_simulate(P, initial=0, steps=40, n_paths=6, random_state=12) + assert np.array_equal(a.data, b.data) + assert not np.array_equal(a.data, c.data) + + +def test_markov_simulate_terminal_frequency_matches_stationary(): + """终点频率收敛到平稳分布 [2/3, 1/3] —— 马尔可夫链模拟的核心性质。""" + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_simulate(P, initial=0, steps=200, n_paths=20000, random_state=0) + freq = np.bincount(r.data[:, -1], minlength=2) / r.data.shape[0] + assert np.allclose(freq, [2 / 3, 1 / 3], atol=0.02) + + +def test_markov_simulate_metrics_are_terminal_states(): + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_simulate(P, initial=0, steps=50, n_paths=500, random_state=3) + assert r.metrics["mean_terminal"] == float(np.mean(r.data[:, -1])) + assert r.metrics["std_terminal"] == float(np.std(r.data[:, -1])) + + +def test_markov_simulate_rejects_bad_transition(): + """非方阵、行和不为 1 都要在入口直接报错，不能悄悄模拟。""" + from mcm_toolkit.solve import markov_simulate + + with pytest.raises(ValueError): + markov_simulate(np.array([[0.5, 0.5, 0.0], [0.2, 0.8, 0.0]]), initial=0) + with pytest.raises(ValueError): + markov_simulate(np.array([[0.5, 0.4], [0.2, 0.8]]), initial=0) + + +def test_markov_simulate_rejects_bad_arguments(): + """§4.7：步数 / 路径数的类型与正负都要卡住（bool 先判，易错点 1）。""" + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + with pytest.raises(ValueError): + markov_simulate(P, initial=0, steps=0) + with pytest.raises(ValueError): + markov_simulate(P, initial=0, n_paths=0) + with pytest.raises(TypeError): + markov_simulate(P, initial=0, steps=1.5) # type: ignore[arg-type] # 就是要传错 + with pytest.raises(TypeError): + markov_simulate(P, initial=0, n_paths=1.5) # type: ignore[arg-type] # 就是要传错 + with pytest.raises(TypeError): + markov_simulate(P, initial=0, steps=True) # type: ignore[arg-type] # 就是要传错 + with pytest.raises(TypeError): + markov_simulate(P, initial=0, n_paths=True) # type: ignore[arg-type] # 就是要传错 + + +def test_markov_simulate_accepts_numpy_integer(): + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_simulate(P, initial=0, steps=np.int64(10), n_paths=np.int64(2), random_state=0) + assert r.success + assert r.data.shape == (2, 11) + + +def test_markov_simulate_docstring_example_runs(): + from mcm_toolkit.solve import markov_simulate + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_simulate(P, initial=0, steps=50, n_paths=4, random_state=0) + assert r.data.shape == (4, 51) + assert np.all(r.data[:, 0] == 0) + + + @pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_markov_steady_state_two_state_chain():