MCM-ICM_Toolkit/tests/test_solve.py

"""solve 包测试。"""

from __future__ import annotations

import numpy as np
import pytest

_NOT_IMPLEMENTED = "未实现：见 docs/03-roadmap.md L0/L1"


def _normal_batch(rng, size, **kwargs):
    """`monte_carlo` 的 sampler 契约：收下实现传来的 `size`，一次抽够。"""
    return rng.normal(0, 1, size)


def test_monte_carlo_integral_converges_to_analytic_value():
    """∫₀¹ x² dx = 1/3。n=2e5 时误差应 < 0.01。"""
    from mcm_toolkit.solve import monte_carlo_integral

    r = monte_carlo_integral(lambda X: X[:, 0] ** 2, bounds=[(0, 1)], n=200_000, random_state=0)
    assert r.success
    assert abs(r.metrics["mean"] - 1 / 3) < 0.01   # 体积必须是正的
    assert len(r.data) == 200_000                 # data 是 n 个样本各自的贡献


def test_monte_carlo_integral_odd_dimension_sign():
    """∫∫∫ x dV over [0,1]³ = 1/2 —— 奇数维最容易暴露"体积符号写反"。"""
    from mcm_toolkit.solve import monte_carlo_integral

    r = monte_carlo_integral(lambda X: X[:, 0], bounds=[(0, 1)] * 3, n=200_000, random_state=0)
    assert abs(r.metrics["mean"] - 0.5) < 0.01


def test_monte_carlo_integral_signed_bounds():
    """∫₀² x dx = 2 —— 边界带宽度时，体积算错和符号写反都会现形。"""
    from mcm_toolkit.solve import monte_carlo_integral

    r = monte_carlo_integral(lambda X: X[:, 0], bounds=[(0, 2)], n=200_000, random_state=0)
    assert abs(r.metrics["mean"] - 2.0) < 0.02


def test_monte_carlo_integral_2d_and_params():
    from mcm_toolkit.solve import monte_carlo_integral

    bounds = [(0, 1), (0, 1)]
    r = monte_carlo_integral(lambda X: X[:, 0] * X[:, 1], bounds=bounds, n=200_000, random_state=0)
    assert abs(r.metrics["mean"] - 0.25) < 0.01
    assert r.name == "monte_carlo_integral"
    assert r.params["dimension"] == 2
    assert r.params["bounds"] == bounds


def test_monte_carlo_integral_rejects_bad_bounds():
    """四类非法 bounds 都要在入口抛 ValueError：区间反了 / 漏括号 / 空 / 参差不齐。"""
    from mcm_toolkit.solve import monte_carlo_integral

    for bad in ([(1, 0)], [0, 1], [], [(0, 1), (0, 2, 3)]):
        with pytest.raises(ValueError):
            monte_carlo_integral(lambda X: X[:, 0], bounds=bad, n=100, random_state=0)


def test_monte_carlo_integral_failure_message_stays_failure():
    """失败时 message 不能被包装层改写成"成功"（易错点 22）。"""
    from mcm_toolkit.solve import monte_carlo_integral

    r = monte_carlo_integral(lambda X: 1.0, bounds=[(0, 1)], n=100, random_state=0)
    assert r.success is False
    assert "成功" not in r.message
    assert "失败" in r.message
    assert r.params["dimension"] == 1          # 失败也要留下追溯信息


def test_monte_carlo_integral_accepts_numpy_integer_n():
    """§4.7：`np.int64` 不是 `int` 的子类，注解与运行期都要收（易错点 1）。"""
    from mcm_toolkit.solve import monte_carlo_integral

    r = monte_carlo_integral(lambda X: X[:, 0], bounds=[(0, 1)], n=np.int64(1000), random_state=0)
    assert r.success
    assert len(r.data) == 1000


def test_monte_carlo_probability_of_standard_normal():
    """P(Z > 0) = 0.5，n=5000 时误差应 < 0.05。"""
    from mcm_toolkit.solve import monte_carlo_probability

    r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=5000, random_state=0)
    assert r.success
    assert abs(r.metrics["probability"] - 0.5) < 0.05


def test_monte_carlo_is_reproducible():
    from mcm_toolkit.solve import monte_carlo_probability

    a = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=1000, random_state=9)
    b = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=1000, random_state=9)
    assert a.metrics["probability"] == b.metrics["probability"]


def test_monte_carlo_returns_n_results_and_documented_metrics():
    from mcm_toolkit.solve import monte_carlo

    r = monte_carlo(_normal_batch, lambda x: x**2, n=2000, random_state=0)
    assert r.success
    assert r.data.shape == (2000,)
    assert set(r.metrics) == {"mean", "std", "median", "p05", "p95", "ci_low", "ci_high"}


def test_monte_carlo_calls_simulator_only_once():
    """`simulator` 只被调用一次 —— 这是"必须向量化"的可执行版本。"""
    from mcm_toolkit.solve import monte_carlo

    counter = {"calls": 0}

    def simulator(params):
        counter["calls"] += 1
        return params

    r = monte_carlo(_normal_batch, simulator, n=1000, random_state=0)
    assert r.success
    assert counter["calls"] == 1


def test_monte_carlo_data_is_reproducible_with_same_seed():
    from mcm_toolkit.solve import monte_carlo

    a = monte_carlo(_normal_batch, lambda x: x, n=2000, random_state=11)
    b = monte_carlo(_normal_batch, lambda x: x, n=2000, random_state=11)
    assert np.array_equal(a.data, b.data)


def test_monte_carlo_follows_set_seed_when_random_state_omitted():
    from mcm_toolkit.core import set_seed
    from mcm_toolkit.solve import monte_carlo

    set_seed(5)
    a = monte_carlo(_normal_batch, lambda x: x, n=500)
    set_seed(5)
    b = monte_carlo(_normal_batch, lambda x: x, n=500)
    assert np.array_equal(a.data, b.data)


def test_monte_carlo_estimates_known_expectation():
    """E[X²] = 1（X ~ N(0, 1)）；n = 2e5 时误差应 < 0.02。"""
    from mcm_toolkit.solve import monte_carlo

    r = monte_carlo(_normal_batch, lambda x: x**2, n=200_000, random_state=0)
    assert abs(r.metrics["mean"] - 1.0) < 0.02


def test_monte_carlo_docstring_example_runs():
    """docstring 的 Examples 必须真的能跑通（此前那一版是坏的）。"""
    from mcm_toolkit.solve import monte_carlo

    draw = lambda rng, size, **k: rng.normal(0, 1, size)  # noqa: E731
    r = monte_carlo(draw, lambda x: x**2, n=10000, random_state=0)
    assert r.success
    assert abs(r.metrics["mean"] - 1.0) < 0.1


def test_monte_carlo_accepts_numpy_integer_n():
    """§4.7：`np.int64` 不是 `int` 的子类，不能漏。"""
    from mcm_toolkit.solve import monte_carlo

    r = monte_carlo(_normal_batch, lambda x: x, n=np.int64(100), random_state=0)
    assert r.success
    assert r.data.shape == (100,)


def test_monte_carlo_counts_boolean_outcomes():
    """布尔结果按 0/1 统计 —— 事件概率是 MC 最常见的用法。"""
    from mcm_toolkit.solve import monte_carlo

    r = monte_carlo(
        lambda rng, size, **k: rng.normal(0, 1, size) > 0,
        lambda x: x,
        n=5000,
        random_state=0,
    )
    assert r.success
    assert r.data.dtype == float
    assert abs(r.metrics["mean"] - 0.5) < 0.05


def test_monte_carlo_rejects_nan_results():
    """结果里出现 NaN 就判失败，不静默给出一列 nan 指标。"""
    from mcm_toolkit.solve import monte_carlo

    r = monte_carlo(
        lambda rng, size, **k: np.r_[np.zeros(size - 1), np.nan],
        lambda x: x,
        n=10,
        random_state=0,
    )
    assert not r.success
    assert "NaN" in r.message


def test_monte_carlo_rejects_bad_arguments():
    from mcm_toolkit.solve import monte_carlo

    with pytest.raises(TypeError):
        monte_carlo(123, lambda x: x)  # type: ignore[arg-type]  # 就是要传错
    with pytest.raises(TypeError):
        monte_carlo(_normal_batch, 123)  # type: ignore[arg-type]  # 就是要传错
    with pytest.raises(TypeError):
        monte_carlo(_normal_batch, lambda x: x, n=True)
    with pytest.raises(TypeError):
        monte_carlo(_normal_batch, lambda x: x, n=10.0)  # type: ignore[arg-type]  # 就是要传错
    with pytest.raises(ValueError):
        monte_carlo(_normal_batch, lambda x: x, n=0)


def test_monte_carlo_probability_returns_documented_metrics():
    from mcm_toolkit.solve import monte_carlo_probability

    r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=20000, random_state=0)
    assert r.success
    assert set(r.metrics) == {"probability", "ci_low", "ci_high", "std_error"}
    assert r.data.shape == (20000,)
    assert 0.0 <= r.metrics["ci_low"] <= r.metrics["probability"]
    assert r.metrics["probability"] <= r.metrics["ci_high"] <= 1.0


def test_monte_carlo_probability_std_error_matches_binomial_formula():
    """`std_error` 就是比例的标准误 √(p(1-p)/n)。"""
    from mcm_toolkit.solve import monte_carlo_probability

    n = 20000
    r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=n, random_state=0)
    p = r.metrics["probability"]
    assert np.isclose(r.metrics["std_error"], np.sqrt(p * (1 - p) / n))


def test_monte_carlo_probability_ci_narrows_with_n():
    """区间宽度随 n 增大而收窄（O(1/√n)）。"""
    from mcm_toolkit.solve import monte_carlo_probability

    def width(n):
        r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=n, random_state=0)
        return r.metrics["ci_high"] - r.metrics["ci_low"]

    assert width(100_000) < width(1_000)


def test_monte_carlo_probability_higher_confidence_is_wider():
    from mcm_toolkit.solve import monte_carlo_probability

    def width(confidence):
        r = monte_carlo_probability(
            lambda z: z > 0, _normal_batch, n=5000, random_state=0, confidence=confidence
        )
        return r.metrics["ci_high"] - r.metrics["ci_low"]

    assert width(0.99) > width(0.90)


def test_monte_carlo_probability_ci_matches_statsmodels_wilson():
    """区间必须对得上独立实现（statsmodels 的 wilson）—— `^` 写成异或那次就靠它拦住。"""
    from statsmodels.stats.proportion import proportion_confint

    from mcm_toolkit.solve import monte_carlo_probability

    n = 5000
    r = monte_carlo_probability(lambda z: z > 0, _normal_batch, n=n, random_state=0)
    count = int(r.data.sum())
    lo, hi = proportion_confint(count, n, alpha=0.05, method="wilson")
    assert np.isclose(r.metrics["ci_low"], lo)
    assert np.isclose(r.metrics["ci_high"], hi)


def test_monte_carlo_probability_rejects_bad_confidence():
    from mcm_toolkit.solve import monte_carlo_probability

    for bad in (0.0, 1.0, -0.5, 1.5):
        with pytest.raises(ValueError):
            monte_carlo_probability(lambda z: z > 0, _normal_batch, n=100, confidence=bad)


def test_random_walk_shape_and_start():
    """`data` 形状 (n_paths, n_steps + 1)，第 0 列是起点。"""
    from mcm_toolkit.solve import random_walk

    r = random_walk(n_steps=100, n_paths=5, start=2.5, random_state=0)
    assert r.success
    assert r.data.shape == (5, 101)
    assert np.all(r.data[:, 0] == 2.5)


def test_random_walk_default_steps_are_plus_minus_one():
    """默认是对称随机游走：每一步都恰好是 ±1。"""
    from mcm_toolkit.solve import random_walk

    r = random_walk(n_steps=200, n_paths=3, random_state=0)
    assert set(np.unique(np.diff(r.data, axis=1))) == {-1.0, 1.0}


def test_random_walk_is_reproducible():
    from mcm_toolkit.solve import random_walk

    a = random_walk(n_steps=50, n_paths=4, random_state=7)
    b = random_walk(n_steps=50, n_paths=4, random_state=7)
    assert np.array_equal(a.data, b.data)


def test_random_walk_accepts_custom_step_sampler():
    """换成高斯增量就是带漂移的高斯游走 —— 建模差异全在 step_sampler 里。"""
    from mcm_toolkit.solve import random_walk

    r = random_walk(
        n_steps=100,
        n_paths=3,
        step_sampler=lambda rng, size: rng.normal(0.5, 1.0, size),
        random_state=0,
    )
    assert r.success
    assert r.data.shape == (3, 101)
    assert abs(r.data[:, -1].mean() - 50.0) < 20.0   # 100 步 × 漂移 0.5


def test_random_walk_rejects_bad_arguments():
    """§4.7 的分工：类型不对抛 TypeError，取值不对（<= 0）抛 ValueError。"""
    from mcm_toolkit.solve import random_walk

    with pytest.raises(TypeError):
        random_walk(n_steps=2.5)  # type: ignore[arg-type]  # 就是要传错
    with pytest.raises(TypeError):
        random_walk(n_steps=True)
    with pytest.raises(TypeError):
        random_walk(n_steps=10, n_paths=1.5)  # type: ignore[arg-type]  # 就是要传错
    with pytest.raises(ValueError):
        random_walk(n_steps=0)
    with pytest.raises(ValueError):
        random_walk(n_steps=10, n_paths=0)
    with pytest.raises(TypeError):
        random_walk(n_steps=10, start="a")  # type: ignore[arg-type]  # 就是要传错


def test_random_walk_rejects_2d_step_sampler():
    """`step_sampler` 必须返回**一维** —— 二维返回（含方阵）一律判失败（易错点 25 的根治）。"""
    from mcm_toolkit.solve import random_walk

    def transposed(rng, size):
        return rng.normal(0, 1, (4, 3))          # 无论要多少，一律给二维

    r = random_walk(n_steps=4, n_paths=3, step_sampler=transposed, random_state=0)
    assert r.success is False                    # 失败也要返回 Result，不能抛给调用方
    assert "长度" in r.message
    assert r.params["n_paths"] == 3              # 失败也要留下追溯信息
    assert r.data is None

    # 旧方案在 n_paths == n_steps 时对"转置"完全无感；现在同样拦住
    assert random_walk(n_steps=3, n_paths=3, step_sampler=transposed, random_state=0).success is False


def test_random_walk_step_sampler_gets_flat_size_and_is_row_major():
    """契约：`size` 是整数 `n_paths * n_steps`，一维增量按**行主序**切成 n_paths 行。"""
    from mcm_toolkit.solve import random_walk

    seen = {}

    def sampler(rng, size):
        seen["size"] = size
        return np.arange(1, size + 1, dtype=float)   # 1, 2, 3, 4, 5, 6

    r = random_walk(n_steps=3, n_paths=2, step_sampler=sampler, random_state=0)
    assert seen["size"] == 6                                  # 整数，不是 (n_paths, n_steps)
    assert r.data.tolist() == [[0, 1, 3, 6], [0, 4, 9, 15]]   # cumsum(1,2,3) / cumsum(4,5,6)


def test_random_walk_rejects_nan_steps():
    """增量含 NaN 判失败 —— 与 monte_carlo 的 NaN 口径（决策 48）一致。"""
    from mcm_toolkit.solve import random_walk

    r = random_walk(
        n_steps=4,
        n_paths=2,
        step_sampler=lambda rng, size: np.full(size, np.nan),
        random_state=0,
    )
    assert r.success is False
    assert "NaN" in r.message


def test_random_walk_accepts_numpy_integer():
    """§4.7：`np.int64` 不是 `int` 的子类，注解与运行期都要收（易错点 1）。"""
    from mcm_toolkit.solve import random_walk

    r = random_walk(n_steps=np.int64(20), n_paths=np.int64(3), random_state=0)
    assert r.success
    assert r.data.shape == (3, 21)


def test_random_walk_metrics_are_terminal_positions():
    """终点位置的标准差按 √t 扩散（t = n_steps）—— 随机游走的核心性质。"""
    from mcm_toolkit.solve import random_walk

    n_steps = 100
    r = random_walk(n_steps=n_steps, n_paths=2000, start=5.0, random_state=0)
    assert abs(r.metrics["mean_terminal"] - 5.0) < 5.0              # 对称游走，期望还在起点附近
    assert abs(r.metrics["std_terminal"] - np.sqrt(n_steps)) < 5.0  # √t 扩散


def test_random_walk_docstring_example_runs():
    from mcm_toolkit.solve import random_walk

    r = random_walk(n_steps=100, n_paths=5, random_state=0)
    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["terminal_distribution"] == {1: 1.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():
    """`metrics` 只看终点：均值 + `{状态: 频率}`（频率之和为 1、与逐条统计一致）。"""
    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)
    terminal = r.data[:, -1]
    dist = r.metrics["terminal_distribution"]
    assert r.metrics["mean_terminal"] == float(np.mean(terminal))
    assert np.isclose(sum(dist.values()), 1.0)
    for state, freq in dist.items():
        assert np.isclose(freq, float(np.mean(terminal == state)))


def test_markov_simulate_terminal_distribution_only_lists_seen_states():
    """只含实际出现过的状态（没出现的状态没有 key）—— docstring 已写明这条口径。"""
    from mcm_toolkit.solve import markov_simulate

    r = markov_simulate(np.eye(3), initial=2, steps=5, n_paths=4, random_state=0)
    assert r.metrics["terminal_distribution"] == {2: 1.0}


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)
    with pytest.raises(ValueError):
        markov_simulate(np.array([[1.5, -0.5], [0.0, 1.0]]), initial=0)  # 行和是 1，但含负元素


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)


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(TypeError):
        markov_steady_state(np.ones((2, 3)) / 3)  # 形状不对 → TypeError（口径见决策 57）
    with pytest.raises(ValueError):
        markov_steady_state(np.array([[0.5, 0.4], [0.2, 0.8]]))
    with pytest.raises(ValueError):
        markov_steady_state(np.array([[1.5, -0.5], [0.0, 1.0]]))  # 行和是 1，但含负元素


def test_markov_steady_state_non_convergence_is_success_false():
    """`max_iter` 不够时按"未收敛"处理：`success=False` 且 `converged=False`（口径 B）。"""
    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=3)
    assert r.success is False
    assert r.metrics["converged"] is False
    assert r.metrics["iterations"] == 3
    assert "未在最大迭代次数内收敛" in r.message


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)
def test_knapsack_dp_known_optimum():
    from mcm_toolkit.solve import knapsack_dp

    r = knapsack_dp(weights=[1, 3, 4, 5], values=[1, 4, 5, 7], capacity=7)
    assert np.isclose(r.metrics["total_value"], 9.0)
    assert np.isclose(r.metrics["total_weight"], 7.0)


@pytest.mark.skip(reason=_NOT_IMPLEMENTED)
def test_linear_program_basic():
    from mcm_toolkit.solve import linear_program

    r = linear_program(c=[1, 2], A_ub=[[1, 1]], b_ub=[4], bounds=[(0, None)] * 2)
    assert r.success
    assert np.isclose(r.metrics["objective"], 4.0)


def test_mm1_metrics_closed_form():
    """与闭式公式逐项对拍（ρ = λ/μ，L = ρ/(1-ρ)，Wq = ρ/(μ-λ)，W = Wq + 1/μ）。"""
    from mcm_toolkit.solve import mm1_metrics

    arrival_rate, service_rate = 2, 3
    rho = arrival_rate / service_rate
    r = mm1_metrics(arrival_rate=arrival_rate, service_rate=service_rate)

    assert r.success
    assert np.isclose(r.metrics["rho"], rho, rtol=1e-12)
    assert np.isclose(r.metrics["L"], rho / (1 - rho), rtol=1e-12)
    assert np.isclose(r.metrics["Lq"], rho**2 / (1 - rho), rtol=1e-12)
    assert np.isclose(r.metrics["W"], 1 / (service_rate - arrival_rate), rtol=1e-12)
    assert np.isclose(r.metrics["Wq"], rho / (service_rate - arrival_rate), rtol=1e-12)
    assert np.isclose(r.metrics["W"], r.metrics["Wq"] + 1 / service_rate, rtol=1e-12)


def test_mm1_rejects_unstable_system():
    """λ >= μ 时系统不稳定，应 success=False 而不是返回 inf（λ == μ 也要走这一支）。"""
    from mcm_toolkit.solve import mm1_metrics

    equal = mm1_metrics(arrival_rate=3, service_rate=3)
    assert equal.success is False
    assert "不稳定" in equal.message
    assert equal.metrics == {"rho": 1.0}  # 失败时也把 rho 报出来（决策 58 的 ④）
    assert equal.data is None

    assert mm1_metrics(arrival_rate=3, service_rate=2).success is False


def test_mm1_allows_zero_arrival_rate_as_empty_system():
    """λ = 0 是合法退化解：空系统（L = Lq = Wq = 0，W 取极限值 1/μ）。"""
    from mcm_toolkit.solve import mm1_metrics

    r = mm1_metrics(arrival_rate=0, service_rate=3)
    assert r.success
    assert r.metrics["rho"] == 0.0
    assert r.metrics["L"] == 0.0
    assert r.metrics["Lq"] == 0.0
    assert r.metrics["Wq"] == 0.0
    assert np.isclose(r.metrics["W"], 1 / 3)


def test_mm1_rejects_bad_arguments():
    """λ 非负、μ 必须为正；bool 与非数值都要在入口拦下。"""
    from mcm_toolkit.solve import mm1_metrics

    with pytest.raises(ValueError):
        mm1_metrics(arrival_rate=-1, service_rate=3)
    with pytest.raises(ValueError):
        mm1_metrics(arrival_rate=2, service_rate=0)
    with pytest.raises(TypeError):
        mm1_metrics(arrival_rate=True, service_rate=3)
    with pytest.raises(TypeError):
        mm1_metrics(arrival_rate=2, service_rate="3")  # type: ignore[arg-type]  # 就是要传错


def test_mm1_docstring_example_runs():
    from mcm_toolkit.solve import mm1_metrics

    r = mm1_metrics(arrival_rate=2, service_rate=3)
    assert round(r.metrics["rho"], 4) == 0.6667


def test_mm1_stability_boundary_is_success_not_exception():
    """λ == μ 不能靠 ZeroDivisionError 落到失败分支：必须是显式的"不稳定"口径。"""
    from mcm_toolkit.solve import mm1_metrics

    r = mm1_metrics(arrival_rate=5, service_rate=5)
    assert r.success is False
    assert "ZeroDivisionError" not in r.message


@pytest.mark.skip(reason=_NOT_IMPLEMENTED)
def test_simulated_annealing_finds_quadratic_minimum():
    from mcm_toolkit.solve import simulated_annealing

    f = lambda v: (v[0] - 3.0) ** 2 + (v[1] + 1.0) ** 2  # noqa: E731
    r = simulated_annealing(f, x0=[0.0, 0.0], n_iter=5000, random_state=0)
    assert np.allclose(r.params["x"], [3.0, -1.0], atol=0.1)