diff --git a/CHANGELOG.md b/CHANGELOG.md index fb286b6..8c315ff 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -529,6 +529,32 @@ 转移矩阵未查非负（行和恰为 1.0 的 `[[1.5, -0.5], [0, 1]]` 会静默给出全 0）、`paths[:, 0]` 赋值在 `try` 之外 （长度不匹配时裸抛 `ValueError`）、状态越界没在入口拦、`std_terminal` 对分类状态的意义。 +### Added +- `solve.markov_steady_state` 落地（**L1，S 档**，实现由用户本人写）：`v` 从均匀分布出发，`v_next = v @ arr` + 反复推进（行随机矩阵按行向量右乘），直到 `max|v_next - v| < tol` 或耗尽 `max_iter`； + `data` 是平稳分布向量、`metrics={"converged", "iterations"}`。Codex 审查结论：**算法本体正确** —— + 3 状态与 `np.linalg.eig(P.T)` 特征值 1 的特征向量逐位一致、200×200 一致到 `1.9e-12`（6 次迭代 / 0.10s）、 + 与 `markov_simulate` 的终点频率互验 `[0.6702, 0.3298]` vs `[0.6667, 0.3333]`。 +- `tests/test_solve.py`：摘掉 `test_markov_steady_state_two_state_chain` 的 skip，另新增 7 条 —— + 与特征分解对拍、`data` 是分布（和为 1 / 非负 / `float64`）、慢混合链迭代次数 > 100、`np.int64` 的 `max_iter`、 + 非法转移矩阵、`max_iter=0` / `tol=0`、docstring Example 真能跑。 + +### Changed +- `markov_steady_state` 的 docstring 补全：`max_iter` / `tol` 两个参数、`Raises` 一节， + 并写清 `tol` 是**相邻差口径**（实测慢混合链 `tol=1e-10` 的实际误差 `4.6e-9`，被 `1/(1-|λ2|)` 放大约 50 倍）。 +- `markov_steady_state` 的 Notes 修正"周期链会不收敛"的说法：可约链的实情是**平稳分布不唯一**， + 迭代法给出的是"离初值最近的那个"（均匀初值对双随机矩阵本身就是不动点，所以 `iterations=1` 不代表链唯一平稳）。 +- `markov_steady_state` 的工程侧清理：ruff 5 条 `W293`；`mm1_metrics` 前 4 个空行 → 2；失败文案 `失败: ` → `失败：`。 +- `STATUS.md`：新增**决策 57**（6 条待拍板口径）与「下一步」第六轮；验证计数更新为 **196 passed / 30 skipped（226 收集）**。 +- `docs/dev/personal-pitfalls.md`：新增易错点 28（数值阈值参数的 `bool` / NaN 漏网）。 + +### Notes +- **`markov_steady_state` 的 6 条口径留给用户**（决策 49）：转移矩阵未查非负（行和恰为 1.0 的 `[[1.5, -0.5], [0, 1]]` + 会溢出成 `nan`，而 `data=[nan, nan]` 却 `success=True`，违反决策 48 的 NaN 口径；附带的 `RuntimeWarning` + 在 pyproject 里被设成 error，**测试会直接报错**）、`tol=True` 被当 `tol=1` 静默放行（`isinstance(True, int)`）、 + `tol=np.nan` 穿过校验白跑迭代、`converged=False` 仍 `success=True`、 + 类型错也抛 `ValueError`（与决策 54 分工不一致）、`params` 没记转移矩阵。 + ## [0.0.1] - 2026-10-01 ### Added diff --git a/STATUS.md b/STATUS.md index d460d0c..b8f8984 100644 --- a/STATUS.md +++ b/STATUS.md @@ -21,10 +21,10 @@ 当前 **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`、`markov_simulate` -（**均为用户写**，2026-10-05）。剩余 1 项：`save_figure`（Codex 写）。 +`monte_carlo_probability`、`monte_carlo_integral`、`random_walk`、`markov_simulate`、 +`markov_steady_state`（**均为用户写**）。剩余 1 项：`save_figure`（Codex 写）。 `monte_carlo` 的口径已拍板（决策 48）；`monte_carlo_probability` / `random_walk` 的待办**均已办结**； -`markov_simulate` 本轮审查通过（算法本体正确），留 **5 条入口口径待你拍板**（见"下一步"，决策 56 的 3 个选项）。 +`markov_simulate`（决策 56，5 条）与 `markov_steady_state`（决策 57，6 条）的**入口口径待你拍板**。 另 3 个小口径仍悬着：`ci_low`/`ci_high` 的语义、`monte_carlo_integral` 的第一道 `isinstance` 分支可塌成两行、 `random_state` 注解与 `std` 的 `ddof` 是否全项目对齐。 @@ -335,13 +335,29 @@ `np.int64` 可收、docstring Example 真能跑 - ⏳ **留给你的 5 条**见"下一步 → markov_simulate 待你拍板"（全是入口口径与少考虑情况，决策 49） +- [x] `solve/simulate.py` 的 `markov_steady_state`（**L1，S 档，算法由用户本人写**，2026-10-09） + - 核心：`v` 从均匀分布出发，`v_next = v @ arr` 反复推进（行随机矩阵按**行向量右乘**）， + 直到 `max|v_next - v| < tol` 或耗尽 `max_iter`；`data` 是平稳分布向量，`metrics={"converged", "iterations"}` + - **Codex 审查结论：算法本体正确**（探针实测）—— 2 状态 `[2/3, 1/3]`（58 次迭代）、3 状态与 `np.linalg.eig(P.T)` + 特征值 1 对应的特征向量逐位一致（39 次）、随机 6 状态一致（13 次）、200×200 一致到 `1.9e-12`（6 次、0.10s）、 + 与 `markov_simulate` 的终点频率互验 `[0.6702, 0.3298]` vs `[0.6667, 0.3333]`、`sum(data) = 1.0` + - **Codex 只碰工程侧**：① ruff 5 条 `W293` 全清，`mm1_metrics` 前 4 个空行 → 2； + ② docstring 补 `max_iter` / `tol` 两个参数（原来 Parameters 里只有 `transition`）、补 `Raises`， + 写清 `tol` 是**相邻差口径**（实测慢混合链 `tol=1e-10` → 实际误差 `4.6e-9`，被 `1/(1-|λ2|)` 放大）； + ③ 修正 Notes 里"周期链会不收敛"的说法（可约链的实情是"平稳分布不唯一、迭代法给离初值最近的那个"）； + ④ 失败文案 `失败: ` → `失败：`（与文件内其它函数统一） + - **测试：+8 条、摘掉 1 条 skip**（此前该函数只有 1 条挂着的 skip）—— 2 状态 `[2/3, 1/3]` + `converged` / `iterations`、 + 3 状态与特征分解对拍、`data` 是分布（和为 1 / 非负 / `float64`）、慢混合链迭代次数 > 100、 + `np.int64` 的 `max_iter`、非法转移矩阵、`max_iter=0` / `tol=0`、docstring Example 真能跑 + - ⏳ **留给你的 6 条口径**见"下一步 → markov_steady_state 待你拍板"（决策 57） + - [ ] 剩余 1 项：`save_figure`（A 档 = 工程，**Codex 直接写**） ### 验证结果 | 检查 | venv | conda `mcm-toolkit` | |---|---|---| -| `pytest` | 见下（单列） | **188 passed, 31 skipped（0 xfail）** | +| `pytest` | 见下（单列） | **196 passed, 30 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` | 跑通 | — | @@ -351,15 +367,15 @@ （修配置之前是 17 errors，全是"无法解析导入 numpy / pandas"，见决策 29 与 `docs/04-environment.md` 坑 7） - base 亦已复核：numpy 2.4.6 / scipy 1.18.1 / pandas 3.0.6（与两套环境一致） -- **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、 +- **226 条用例：196 passed / 30 skipped / 0 xfail**。收集数按文件：`test_data.py` 52、`test_api_surface.py` 20、 + `test_core.py` 23、`test_solve.py` 58（本轮 +8、摘掉 1 条 skip）、`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`、`markov_steady_state`、`mm1_metrics`、`anova`、`fit_distribution`、 +- 30 条 skip 是等实现的算法（`train_test_split`、`mm1_metrics`、`anova`、`fit_distribution`、 `knapsack_dp`、`linear_program`、`simulated_annealing` 等）； - `monte_carlo_probability` / `monte_carlo_integral` 已落地，分别摘掉 2 条与 1 条 skip；`markov_simulate` 新增的 9 条全部实跑； + `monte_carlo_probability` / `monte_carlo_integral` / `markov_steady_state` 已落地，分别摘掉 2 / 1 / 1 条 skip；`markov_simulate` 新增的 9 条全部实跑； `TestValidate` 的 4 条仍挂着 —— `check_*` 还是 `NotImplementedError`，但它们实现时应当直接复用 `as_array` - **本轮回归（已修）**：`transform.py` 上一版有 3 处真 bug —— 修前 **7 failed / 75 passed**， 且 pyright 在同一处报 `reportArgumentType`。**不是静态误报**：pyright CLI 复现 + 测试真红。 @@ -527,17 +543,29 @@ docstring 与 `test_describe_reports_expected_keys` 早就写死了这个方向 `True` 因 `isinstance(True, int)` 为真被当状态 1 **静默放行**；越界状态（如 `initial=5`、`s=2`）只在 第一轮索引时才炸、被 `except` 吞成 `success=False`。可选口径见下一步的 A / B / C。 +| 57 | `markov_steady_state` 的 6 条口径**待用户拍板**（Codex 只报选项，未代改，决策 49） | +探针实测出的口子：① **含负元素**的"行和为 1"矩阵（`[[1.5, -0.5], [0, 1]]`）会溢出成 `nan`，而 +`data=[nan, nan]` 却 `success=True` —— 违反决策 48 的 NaN 口径（根因：没查 `arr >= 0`）；附带的 +`RuntimeWarning: overflow encountered in matmul` 因为 pyproject 设了 `filterwarnings = ["error::RuntimeWarning"]`， +**能让测试直接报错**，所以不能只靠测试兜；② `tol=True` 因 `isinstance(True, int)` 通过校验 → `tol=1` → +1 次就"收敛"，给出 `[0.55, 0.45]` 的假答案（`max_iter` 查了 `bool`、`tol` 没查，易错点 1 复发）； +③ `tol=np.nan` 因 `nan <= 0` 为假穿过校验 → 白跑 10000 次迭代；④ `converged=False` 仍 `success=True` +且 `data` 是未收敛的迭代值（`max_iter=5` → `[0.6387, 0.3613]`）；⑤ 类型错也抛 `ValueError` +（`max_iter=100.0`），与决策 54 的"类型错 `TypeError`"分工不一致；⑥ `params` 只记 `max_iter` / `tol`， +没记转移矩阵（`monte_carlo_integral` 记了 `bounds`），事后复现不了。可选口径见"下一步"。 + ## 下一步 -### ⭐ 当前进度与待办（2026-10-05 第五轮） +### ⭐ 当前进度与待办（2026-10-09 第六轮） -**本轮：`markov_simulate`（你写的第一版）已审查通过** —— 算法本体正确，没动你任何算法行； -Codex 只清工程侧（ruff 8 条、docstring、3 处 `马尔科夫`、`initail` 拼写、半角逗号、`kind` 未用）， -并补 **9 条测试**（此前 0 条）。三连：`pytest` **188 passed / 31 skipped**、`ruff` 0、`pyright` 0。 +**本轮：`markov_steady_state`（你写的）已审查通过** —— 算法本体正确，一行没动（3 状态与特征分解逐位一致、 +200×200 只用 6 次迭代 / 0.10s）；Codex 只清工程侧（ruff 5 条、docstring 补 `max_iter`/`tol`/`Raises`、 +修正"周期链"表述、全角冒号），并 **+8 条测试、摘掉 1 条 skip**。 +三连：`pytest` **196 passed / 30 skipped**、`ruff` 0、`pyright` 0。 **接下来三件事** -1. **`markov_simulate` 的 5 条待你拍板**（下节，全是入口口径，决策 49/56）。 -2. `ci_low` / `ci_high` 的语义确认一句（A / B / C 三个选项在下面挂着），确认后 Codex 补 `docs/algorithms/monte_carlo.md` 卡片。 +1. **`markov_steady_state` 的 6 条口径待你拍板**（下节，决策 57）。 +2. **`markov_simulate` 的 5 条口径**仍在挂着（决策 56 的 A / B / C），与 `solve` 的收尾一起定。 3. L0 最后一项 `save_figure`（A 档 = 工程，`viz/export.py`）—— **Codex 直接实现**，你只读。 ### ⏳ markov_simulate 待你拍板（5 条，Codex 只报不改） @@ -556,9 +584,31 @@ Codex 只清工程侧（ruff 8 条、docstring、3 处 `马尔科夫`、`initail **⑤ `std_terminal` 对**分类**状态的意义**：状态编号是标称量，均值和标准差只对"有序状态"（如排队人数）成立 —— 要不要把 `metrics` 换成状态频率（`terminal_freq`）？你定。 -**你下一个要写的（S 档，仍在 `solve/simulate.py`）**：`markov_steady_state`。 -与 `markov_simulate` 是同一件事的两面：前者按转移矩阵逐代抽状态，后者用**迭代**求平稳分布。 -docstring 已写明口径：用迭代而不用特征分解，因为"收敛快不快"本身就是论文里要写的性质。 +### ⏳ markov_steady_state 待你拍板（6 条，Codex 只报不改） + +**① 转移矩阵没查非负**（与 `markov_simulate` 的待办 ② 是同一块）：`[[1.5, -0.5], [0, 1]]` 行和恰为 1.0 → +通过校验 → `v` 指数发散 → `data=[nan, nan]` **但 `success=True`**（违反决策 48 的 NaN 口径）。 +实测还会先抛 `RuntimeWarning: overflow encountered in matmul`，而 pyproject 设了 +`filterwarnings = ["error::RuntimeWarning"]` —— 这条**能让测试直接报错**，所以得在校验层拦，不能只靠测试兜。 + +**② `tol` 没查 `bool`**：`tol=True` → `isinstance(True, int)` 为真 → `tol=1` → +1 次就"收敛"，给出 `[0.55, 0.45]` 的**假答案**（`success=True, converged=True, iterations=1`）。 +`max_iter` 查了 `bool`，`tol` 漏了（易错点 1 复发）。 + +**③ `tol` 没查 NaN**：`np.nan <= 0` 为假 → 穿过校验 → 白跑 10000 次迭代后返回（结果本身是对的，纯浪费）。 + +**④ `converged=False` 仍 `success=True`**，且 `data` 是未收敛的迭代值（`max_iter=5` → `[0.6387, 0.3613]`）。 +选项：**A** 保持现状（`converged` 自己看）；**B** 不收敛判 `success=False`；**C** 不收敛时 `data=None`、只留 `metrics`。 + +**⑤ 类型错也抛 `ValueError`**（`max_iter=100.0`），与文件内其它函数（决策 54：类型错 `TypeError`、取值错 `ValueError`） +不一致 —— 要不要对齐？`tol` 同理。 + +**⑥ `params` 没记转移矩阵**：`{"max_iter": 10000, "tol": 1e-10}` 复现不了"算的是哪个矩阵" +（`monte_carlo_integral` 记了 `bounds`）—— 要不要补？ + +**你下一个要写的（S 档，仍在 `solve/simulate.py`）**：`mm1_metrics`（M/M/1 闭式指标，10 行）。 +docstring 已写明口径：`λ >= μ` 时系统不稳定、队列长度无界 —— 这种情况返回 `success=False` 并说明原因， +而不是返回 `inf`。写完 `solve` 这一层就齐了，L1 只剩 `decision` / `models` 那几条。 ### ✅ random_walk 的 4 条待办（2026-10-05 已办结） diff --git a/docs/dev/personal-pitfalls.md b/docs/dev/personal-pitfalls.md index 94198fa..6cf04ed 100644 --- a/docs/dev/personal-pitfalls.md +++ b/docs/dev/personal-pitfalls.md @@ -468,3 +468,22 @@ 对不上就说明口径没定清，**先定语义再写代码**；② 分流用的 `isinstance` 必须把 `np.integer` / `bool` 一起想到（易错点 1）； ③ 任何"写进结果数组"的赋值，问一遍"它的长度约束和上面校验用的是同一个下标维度吗"。 - 首次：2026-10-05｜复发：0 + +### 28. 数值阈值参数的两个漏网：`bool` 和 NaN + +- 现象：`markov_steady_state` 的 `tol` 只写了 `if not isinstance(tol, (int, float)) or tol <= 0`。 + `tol=True` → `isinstance(True, int)` 为真、`True > 0` 也为真 → 阈值变成 `1`， + **第 1 次迭代就"收敛"**，返回 `[0.55, 0.45]` 这种假答案，而且 `success=True, converged=True`。 + `tol=np.nan` → `nan <= 0` 为假 → 穿过校验 → 白跑满 `max_iter` 次迭代。 + 同一个函数里 `max_iter` 查了 `bool`、`tol` 没查 —— 说明是**逐条手动补的**，不是照着清单过的。 +- 为什么会错：脑子里"阈值 = 一个正数"，于是校验只写"是数 & 大于 0"；而 `bool` 是 `int` 的子类、 + NaN 与任何数比较都是 `False`，这两条都恰好从"是数 & 大于 0"里漏出去。**与易错点 1 同根**， + 只是这次漏的是"比较"而不是"类型判断"。 +- 同族（同一轮发现的）：只校验"转移矩阵行和 = 1.0"就以为矩阵合法 —— `[[1.5, -0.5], [0, 1]]` 行和也是 1.0， + 但含负元素，迭代时 `v` 指数发散、最后 `data=[nan, nan]` 却 `success=True`。 + **"行和 = 1"不等于"是概率矩阵"**，还要 `arr >= 0`。另外这类溢出会抛 `RuntimeWarning`， + 而本项目 `filterwarnings = ["error::RuntimeWarning"]` —— **测试会直接报错**，别指望测试能兜住静默错。 +- 以后怎么防：① 校验数值阈值时，脑子里过一遍"三件套"：**类型（含 `bool` / `np.integer` 例外）、 + 有限性（`np.isfinite`，挡 NaN / inf）、取值区间**；② 同一函数里同类参数（`max_iter` / `tol` / `n` / `tol`） + 用**同一套模板**写，别一个查 `bool` 一个不查；③ 校验"概率/权重"这类输入时，除了和为 1，还要问一句"有没有可能是负数/NaN"。 +- 首次：2026-10-09｜复发：0 diff --git a/src/mcm_toolkit/solve/simulate.py b/src/mcm_toolkit/solve/simulate.py index 08e2f52..37ee2b8 100644 --- a/src/mcm_toolkit/solve/simulate.py +++ b/src/mcm_toolkit/solve/simulate.py @@ -590,23 +590,42 @@ def markov_simulate( ) -def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e-10) -> Any: - """平稳分布。 +def markov_steady_state(transition: Any, max_iter: int | np.integer = 10000, tol: float = 1e-10) -> Any: + """迭代法求马尔可夫链的平稳分布（`v @ P` 反复推进，直到相邻两次几乎不变）。 Parameters ---------- transition: - 行随机矩阵（每行和为 1）。 + 行随机矩阵（每行和为 1）的二维方阵。 + max_iter: + 迭代上限，正整数，默认 10000。 + tol: + 收敛阈值，正数，默认 1e-10。判据是**相邻两次迭代的最大分量差** + `max|v_k+1 - v_k| < tol`，不是"到不动点的距离"，见 Notes。 Returns ------- - `Result`：`data` 是平稳分布向量，`metrics={"converged", "iterations"}`。 + `Result`：`data` 是一维平稳分布向量（和为 1）； + `metrics={"converged", "iterations"}`（是否在 `max_iter` 内向容、实际迭代次数）。 + + Raises + ------ + ValueError: + `transition` 不是二维方阵、或某一行和不等于 1.0； + `max_iter` 不是正整数；`tol` 不是正数。 Notes ----- 用迭代法而不是特征分解：迭代法能顺便告诉你"收敛快不快"， 而收敛速度本身就是论文里值得写一句的性质。 - **不可约 + 非周期**才有唯一平稳分布，周期链会不收敛，要检查。 + 初值是**均匀分布** `1/s`；**不可约 + 非周期**才有唯一平稳分布。 + 可约链（如单位阵 `eye(2)`）的平稳分布不唯一，迭代法给出的是"离初值最近的那个" + —— 均匀初值对双随机矩阵本身就是不动点，所以 `converged=True, iterations=1` 只说明"初值恰是不动点"， + 不代表这条链有唯一平稳分布；周期链从一般初值出发会一直振荡、`converged=False`。 + + `tol` 是**相邻差**口径，实际误差约被 `1 / (1 - |λ2|)` 放大：实测双状态慢混合链 + `[[0.999, 0.001], [0.02, 0.98]]`（`|λ2| = 0.979`）用默认 `tol=1e-10` 停在第 867 次， + 与特征向量相差 `4.6e-9`（约 50 倍）。慢混合链要更准，把 `tol` 压小即可。 Examples -------- @@ -616,8 +635,57 @@ def markov_steady_state(transition: Any, max_iter: int = 10000, tol: float = 1e- >>> bool(np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6)) True """ - # TODO(L1) - raise NotImplementedError + 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] + v = np.full(s, 1.0 / s) + + iterations = 0 + converged = False + + if isinstance(max_iter, bool) or not isinstance(max_iter, (int, np.integer)) or max_iter <= 0: + raise ValueError(f"max_iter 必须是正整数，你传入了 {max_iter}") + if not isinstance(tol, (int, float)) or tol <= 0: + raise ValueError(f"tol 必须是正数，你传入了 {tol}") + + try: + for _ in range(max_iter): + v_next = v @ arr + error = np.max(np.abs(v_next - v)) + v = v_next + iterations += 1 + + if error < tol: + converged = True + break + + metrics = { + "converged": converged, + "iterations": iterations + } + + return Result( + name="markov_steady_state", + params={"max_iter": max_iter, "tol": tol}, + metrics=metrics, + data=v, + success=True, + message=f"平稳分布求解{'成功' if converged else '未完全收敛'}，迭代 {iterations} 次。" + ) + + except Exception as e: + return Result( + name="markov_steady_state", + params={"max_iter": max_iter, "tol": tol}, + metrics={}, + data=None, + success=False, + message=f"平稳分布求解失败：{type(e).__name__} - {str(e)}" + ) def mm1_metrics(arrival_rate: float, service_rate: float) -> Any: diff --git a/tests/test_solve.py b/tests/test_solve.py index 5549176..84d1df3 100644 --- a/tests/test_solve.py +++ b/tests/test_solve.py @@ -508,15 +508,88 @@ def test_markov_simulate_docstring_example_runs(): assert np.all(r.data[:, 0] == 0) - - -@pytest.mark.skip(reason=_NOT_IMPLEMENTED) def test_markov_steady_state_two_state_chain(): """P=[[0.9,0.1],[0.2,0.8]] 的平稳分布是 [2/3, 1/3]。""" from mcm_toolkit.solve import markov_steady_state r = markov_steady_state(np.array([[0.9, 0.1], [0.2, 0.8]])) + assert r.success assert np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6) + assert r.metrics["converged"] is True + assert r.metrics["iterations"] > 0 + + +def test_markov_steady_state_matches_eigenvector(): + """与特征分解对拍：取特征值 1 对应的左特征向量（只是检验手段，实现仍走迭代）。""" + from mcm_toolkit.solve import markov_steady_state + + P = np.array([[0.5, 0.25, 0.25], [0.1, 0.8, 0.1], [0.2, 0.2, 0.6]]) + eigenvalues, eigenvectors = np.linalg.eig(P.T) + ref = np.real(eigenvectors[:, int(np.argmin(np.abs(eigenvalues - 1.0)))]) + ref = ref / ref.sum() + + r = markov_steady_state(P) + assert np.allclose(r.data, ref, atol=1e-6) + + +def test_markov_steady_state_is_a_distribution(): + from mcm_toolkit.solve import markov_steady_state + + P = np.array([[0.7, 0.2, 0.1], [0.1, 0.6, 0.3], [0.05, 0.15, 0.8]]) + r = markov_steady_state(P) + assert r.data.shape == (3,) + assert np.isclose(float(np.sum(r.data)), 1.0) + assert np.all(r.data >= 0.0) + assert r.data.dtype == np.float64 + + +def test_markov_steady_state_slow_mixing_needs_more_iterations(): + """慢混合链：迭代次数显著变多，但结果仍要对得上特征向量（tol 是相邻差口径）。""" + from mcm_toolkit.solve import markov_steady_state + + P = np.array([[0.999, 0.001], [0.02, 0.98]]) + r = markov_steady_state(P) + assert r.metrics["iterations"] > 100 + assert np.allclose(r.data, [0.9523809524, 0.0476190476], atol=1e-6) + + +def test_markov_steady_state_accepts_numpy_integer(): + """§4.7：`np.int64` 不是 `int` 的子类，注解与运行期都要收（易错点 1）。""" + from mcm_toolkit.solve import markov_steady_state + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_steady_state(P, max_iter=np.int64(10000)) + assert r.success + assert r.metrics["converged"] is True + + +def test_markov_steady_state_rejects_bad_transition(): + from mcm_toolkit.solve import markov_steady_state + + with pytest.raises(ValueError): + markov_steady_state(np.ones((2, 3)) / 3) + with pytest.raises(ValueError): + markov_steady_state(np.array([[0.5, 0.4], [0.2, 0.8]])) + + +def test_markov_steady_state_rejects_bad_arguments(): + """取值错一律 `ValueError`（类型错要不要拆 `TypeError` 待定，见 STATUS 决策 57）。""" + from mcm_toolkit.solve import markov_steady_state + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + with pytest.raises(ValueError): + markov_steady_state(P, max_iter=0) + with pytest.raises(ValueError): + markov_steady_state(P, tol=0) + + +def test_markov_steady_state_docstring_example_runs(): + from mcm_toolkit.solve import markov_steady_state + + P = np.array([[0.9, 0.1], [0.2, 0.8]]) + r = markov_steady_state(P) + assert bool(np.allclose(r.data, [2 / 3, 1 / 3], atol=1e-6)) + @pytest.mark.skip(reason=_NOT_IMPLEMENTED)