MCM-ICM_Toolkit/src/mcm_toolkit/solve/simulate.py

"""随机模拟：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

import numpy as np
from scipy.stats import norm

from mcm_toolkit.core import Result, as_array, get_rng

__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 | np.integer = 10000,
    random_state: int | None = None,
    **kwargs: Any,
) -> Any:
    """通用蒙特卡洛：采样 → 仿真 → 汇总。

    Parameters
    ----------
    sampler:
        `f(rng, size, **kwargs) -> 一批参数`。**一次就要抽够 `n` 组**，
        例如 `lambda rng, size: rng.normal(0, 1, size)`。
        批量大小由本函数以 `size=n` 传进去：`sampler` 必须能接住 `size`
        （显式声明，或带 `**kwargs`）。**别再往 `**kwargs` 里塞 `size`**，否则重复传参。
    simulator:
        `f(params) -> 一批结果`。收到 `sampler` 的那一整批参数，
        返回长度同为 `n` 的数组（标量 / 长度对不上都算契约错误）。
        **内部不要写 Python 循环**。
    n:
        模拟次数，正整数。与 `simulator` 返回的长度必须一致。
    random_state:
        随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。
    **kwargs:
        原样转交给 `sampler`，例如分布参数 `mu` / `sigma`（`size` 由本函数自己传，见上）。
        会一并记进 `Result.params`，便于复现。

    Returns
    -------
    `Result`：`data` 是 `n` 个结果组成的数组（一律转成 float，布尔结果按 0/1 统计），
    `metrics={"mean", "std", "median", "p05", "p95", "ci_low", "ci_high"}`。

    Raises
    ------
    TypeError:
        `sampler` / `simulator` 不可调用，`n` 不是整数，或 `random_state` 类型非法。
    ValueError:
        `n` 不是正整数。

    Notes
    -----
    比赛里 Monte Carlo 的收敛速度是 `O(1/√n)` —— 想少一半误差就要 4 倍样本。
    所以论文里要给"区间随 n 收敛"的图，证明 `n` 已经够大
    （这也是 `viz.plot_convergence` 的典型用法）。

    结果里只要出现 **NaN** 就判失败（`success=False`）：MC 的 NaN 通常来自除零、
    取负数的根号或模型发散，静默统计只会得到一列 nan 指标。

    Examples
    --------
    >>> import numpy as np
    >>> draw = lambda rng, size, **k: rng.normal(0, 1, size)
    >>> r = monte_carlo(draw, lambda x: x ** 2, n=10000, random_state=0)
    >>> bool(abs(r.metrics["mean"] - 1.0) < 0.1)  # E[X²] = 1
    True
    """
    if not callable(sampler) or not callable(simulator):
        raise TypeError("sampler 和 simulator 必须是可调用函数")
    if isinstance(n, bool) or not isinstance(n, (int, np.integer)):
        raise TypeError(f"n 必须是整数，你传入了 {type(n).__name__}")
    if int(n) <= 0:
        raise ValueError(f"n 必须是正整数，你传入了 {n}")
    n = int(n)

    rng = get_rng(random_state)
    call_params = {"n": n, "random_state": random_state, **kwargs}

    try:
        params = sampler(rng, size=n, **kwargs)
        data = np.asarray(simulator(params), dtype=float).ravel()
        if np.isnan(data).any():
            raise ValueError("模拟结果中包含 NaN，请检查 sampler 或 simulator 的传入参数")
        if len(data) != n:
            raise ValueError(f"simulator 返回的结果数量({len(data)})必须与模拟次数 n({n})一致")

        metrics = {
            "mean": float(np.mean(data)),
            "std": float(np.std(data, ddof=0)),
            "median": float(np.median(data)),
            "p05": float(np.percentile(data, 5)),
            "p95": float(np.percentile(data, 95)),
            "ci_low": float(np.percentile(data, 2.5)),
            "ci_high": float(np.percentile(data, 97.5)),
        }

        return Result(
            name="monte_carlo",
            params=call_params,
            metrics=metrics,
            data=data,
            success=True,
            message=f"蒙特卡洛模拟成功：已完成 {n} 次模拟",
        )
    except Exception as e:
        return Result(
            name="monte_carlo",
            params=call_params,
            metrics={},
            data=None,
            success=False,
            message=f"蒙特卡洛模拟失败：{type(e).__name__} - {e}",
        )


def monte_carlo_probability(
    event: Callable[..., Any],
    sampler: Callable[..., Any],
    n: int | np.integer = 10000,
    confidence: float = 0.95,
    random_state: int | None = None,
) -> Any:
    """估计事件发生概率，并给出区间。

    Parameters
    ----------
    event:
        `f(params) -> 一批 bool`。吃到 `sampler` 的**整批**参数，
        返回长度为 `n` 的布尔数组（写成 `z > 0` 这种逐元素判断即可，不要写 Python 循环）。
    sampler:
        `f(rng, size, **kwargs) -> 一批参数`，契约同 `monte_carlo`：
        批量大小由本函数以 `size=n` 传入，`sampler` 必须能接住 `size`。
    n:
        模拟次数，正整数。与 `event` 返回的长度必须一致。
    confidence:
        区间置信水平，必须落在开区间 `(0, 1)` 内。
    random_state:
        随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。

    Returns
    -------
    `Result`：`data` 是 `n` 个 0/1 结果组成的数组，
    `metrics={"probability", "ci_low", "ci_high", "std_error"}`。
    区间是 Wilson score 区间（`p` 贴近 0 / 1 时比 Wald 稳），并夹在 `[0, 1]` 内（比例量纲，见决策 42）。

    Raises
    ------
    TypeError:
        `n` 不是整数，或 `random_state` 类型非法。
    ValueError:
        `n` 不是正整数，或 `confidence` 不在 `(0, 1)` 内。

    Notes
    -----
    概率估计的相对误差 ≈ `1/√(n·p)`。**小概率事件（p < 1e-3）用朴素 MC 是不可行的**，
    需要重要性采样或方差缩减——遇到这类题要提前意识到，别硬跑。
    本函数当前的口径：`probability <= 1e-3` 时**直接判 `success=False`**（只给提示，不给区间）。

    Examples
    --------
    >>> import numpy as np
    >>> draw = lambda rng, size, **k: rng.normal(0, 1, size)
    >>> r = monte_carlo_probability(lambda z: z > 0, draw, n=2000, random_state=0)
    >>> abs(r.metrics["probability"] - 0.5) < 0.1
    True
    """
    if not callable(sampler) or not callable(event):
        raise TypeError("sampler 和 event 必须是可调用函数")
    if isinstance(n, bool) or not isinstance(n, (int, np.integer)):
        raise TypeError(f"n 必须是整数，你传入了 {type(n).__name__}")
    if int(n) <= 0:
        raise ValueError(f"n 必须是正整数，你传入了 {n}")
    n = int(n)
    if confidence <= 0.0 or confidence >= 1:
        raise ValueError(f"confidence必须在0和1之间,而你输入了{confidence}")

    rng = get_rng(random_state)
    call_params = {"n": n, "random_state": random_state, "confidence": confidence}

    try:
        params = sampler(rng, size=n)
        data = np.asarray(event(params), dtype=float).ravel()
        if np.isnan(data).any():
            raise ValueError("模拟结果中包含 NaN，请检查 event 或 sampler 的传入参数")
        if len(data) != n:
            raise ValueError(f"event 返回的结果数量({len(data)})必须与模拟次数 n({n})一致")
        p_value = float(np.mean(data))
        std_error = float(np.sqrt(p_value * (1.0 - p_value) / n))

        z_score = norm.ppf(1 - (1 - confidence) / 2)
        center = (p_value + z_score**2 / (2 * n)) / (1 + z_score ** 2 / n)
        half_width = (z_score / (1 + z_score ** 2 / n)) * np.sqrt(p_value * (1 - p_value) / n + z_score ** 2 / (4 * n ** 2))
        ci_low = max(0.0, center - half_width)
        ci_high = min(1.0, center + half_width)

        warning_message = ""
        if p_value <= 1e-3 or p_value >= 1 - 1e-3:
            warning_message = "小概率事件 P < 1e-3 或者高概率事件 P > 1 - 1e-3,不建议使用monte_carlo_probability"

        metrics = {
            "probability": p_value,
            "ci_low": float(ci_low),
            "ci_high": float(ci_high),
            "std_error": std_error,
        }
        return Result(
            name="monte_carlo_probability",
            params=call_params,
            metrics=metrics,
            data=data,
            success=True,
            message=(
                f"蒙特卡洛模拟成功：已完成 {n} 次模拟"
                + (f"；{warning_message}" if warning_message else "")
            ),
        )
    except Exception as e:
        return Result(
            name="monte_carlo_probability",
            params=call_params,
            metrics={},
            data=None,
            success=False,
            message=f"蒙特卡洛模拟失败：{type(e).__name__} - {e}",
        )


def monte_carlo_integral(
    func: Callable[..., Any],
    bounds: Any,
    n: int | np.integer = 10000,
    random_state: int | None = None,
) -> Any:
    """用蒙特卡洛求多重积分。

    Parameters
    ----------
    func:
        `f(X) -> 一批函数值`。`X` 是形状 `(n, d)` 的样本矩阵（一行一个样本），
        返回长度 `n` 的数组（写成 `X[:, 0] ** 2` 这种逐行写法，不要写 Python 循环）。
    bounds:
        `[(lo, hi), ...]`（非空 list / tuple），维度即积分重数；要求每项都是二元组且 `lo < hi`。
    n:
        模拟次数，正整数。
    random_state:
        随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。

    Returns
    -------
    `Result`：**`metrics["mean"]` 才是积分估计值** ≈ `∫f dV`；
    `data` 是 `n` 个样本各自的贡献（`func(X_i) × 体积`），可拿去画随 `n` 的收敛图。
    `params` 里额外记了 `bounds` 与 `dimension`。

    Raises
    ------
    ValueError:
        `bounds` 不是非空 list / tuple、形状不是 `(d, 2)`（例如漏括号写成 `[0, 1]`），或某个 `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.metrics["mean"] - 1 / 3) < 0.01   # ∫₀¹ x² dx = 1/3
    True
    """
    if not isinstance(bounds, np.ndarray):
        if not isinstance(bounds, (list, tuple)) or len(bounds) == 0:
            raise ValueError("bounds 必须是非空的列表或元组")

        bounds_arr = np.asarray(bounds)
    else:
        bounds_arr = bounds.copy()

    if bounds_arr.ndim != 2 or bounds_arr.shape[1] != 2:
        raise ValueError(f"bounds 必须形如 [(lo, hi), ...]，当前形状为 {bounds_arr.shape}，请检查是否漏了括号")

    lows = bounds_arr[:, 0]
    highs = bounds_arr[:, 1]
    d = len(lows)
    if np.any(lows >= highs):
        raise ValueError("积分边界非法：下界必须严格小于上界")

    widths = highs - lows

    volume = np.prod(widths)

    def sampler(rng, size=n):
        return rng.uniform(low=lows, high=highs, size=(size, d))

    def simulator(X):
        vals = func(X)
        if len(vals) != n:
            raise ValueError(f"func 返回的结果数量({len(vals)})必须与模拟次数 n({n})一致")
        return np.asarray(vals, dtype=float).ravel() * volume

    result = monte_carlo(
        sampler=sampler,
        simulator=simulator,
        n=n,
        random_state=random_state,
    )

    result.name = "monte_carlo_integral"
    result.params["bounds"] = bounds
    result.params["dimension"] = d
    if result.success:
        result.message = f"蒙特卡洛多重积分成功：维度 {d}，体积 {volume:.4f}，模拟 {n} 次。"
    else:
        result.message = f"蒙特卡洛多重积分失败：维度 {d}，体积 {volume:.4f}，模拟 {n} 次。{result.message}"

    return result


def random_walk(
    n_steps: int | np.integer,
    n_paths: int | np.integer = 1,
    step_sampler: Callable[..., Any] | None = None,
    start: float = 0.0,
    random_state: int | None = None,
) -> Any:
    """随机游走轨迹：一次抽完整批增量，再 `cumsum` 成路径。

    Parameters
    ----------
    n_steps:
        每条路径的步数，正整数。`data` 的列数是 `n_steps + 1`（多出的那一列是起点）。
    n_paths:
        路径条数，正整数，默认 1 条。
    step_sampler:
        `f(rng, size) -> 一批增量`。**`size` 是整数** `n_paths * n_steps`（与 `monte_carlo` 的 `sampler` 同一契约），
        必须返回**一维**长度 `n_paths * n_steps` 的数组；本函数负责 `reshape(n_paths, n_steps)`（行主序）。
        返回二维一律判失败 —— 否则"转置的采样器"在 `n_paths == n_steps` 时会静默给出转置结果。
        默认 ±1 等概率（对称随机游走）。
    start:
        起点，写进每条路径的第 0 列。
    random_state:
        随机种子。`None` 时跟随 `core.set_seed` 的全局流，见 `core.random.get_rng`。

    Returns
    -------
    `Result`：`data` 形状 `(n_paths, n_steps + 1)`，第 0 列恒为 `start`，之后是逐列累加出来的位置；
    `metrics={"mean_terminal", "std_terminal"}`（终点的均值与标准差，`std` 用 `ddof=0`）。

    Raises
    ------
    TypeError:
        `n_steps` / `n_paths` 不是整数。
        `start` 不是数值。
    ValueError:
        `n_steps` / `n_paths` 不是正整数。

    Notes
    -----
    随机游走是"平方根扩散"的原型：位置的标准差随步数按 `√t` 增长，所以典型图形是 `±√t` 的喇叭口。
    本函数只负责 `cumsum`；**建模差异全在 `step_sampler` 里** —— 换成 `rng.normal(0, 1, size)` 就是高斯游走，
    换成带正漂移的增量就是有趋势的游走。
    `step_sampler` 返回的不是**一维**、长度不等于 `n_paths * n_steps`，或增量里出现 NaN，都判 `success=False`
    （NaN 口径与 `monte_carlo` 一致）；其余异常同样统一走失败分支。
    步与路径的对应关系是**行主序**：一维数组的前 `n_steps` 个元素属于第 1 条路径。

    Examples
    --------
    >>> import numpy as np
    >>> r = random_walk(n_steps=100, n_paths=5, random_state=0)
    >>> r.data.shape
    (5, 101)
    >>> bool(np.all(np.abs(np.diff(r.data, axis=1)) == 1))   # 默认步长恒为 1
    True
    """
    if isinstance(n_steps, bool) or not isinstance(n_steps, (int, np.integer)):
        raise TypeError(f"n_steps 必须是整数，你传入了 {type(n_steps).__name__}")
    if n_steps <= 0:
        raise ValueError(f"n_steps 必须是正整数，你传入了 {n_steps}")

    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 not isinstance(start, (int, float, np.number)):
        raise TypeError(f"start 必须是数值，你传入了 {type(start).__name__}")

    rng = get_rng(random_state)

    if step_sampler is None:
        def _default_step_sampler(rng, size):
            return rng.choice([-1.0, 1.0], size=size)

        step_sampler = _default_step_sampler

    call_params = {
        "n_steps": n_steps,
        "n_paths": n_paths,
        "start": start,
        "random_state": random_state,
    }

    try:
        steps = np.asarray(step_sampler(rng, size=n_paths * n_steps), dtype=float)
        if steps.shape != (n_paths * n_steps,):
            raise ValueError(
                f"step_sampler 返回长度必须为 {n_paths * n_steps}，实际得到了 {steps.shape}"
            )
        if np.isnan(steps).any():
            raise ValueError("step_sampler 生成的增量数据中包含 NaN，请检查采样逻辑")

        steps = steps.reshape(n_paths, n_steps)   # 行主序：前 n_steps 个属于第 1 条路径

        cumsum = np.cumsum(steps, axis=1)

        paths = np.zeros((n_paths, n_steps + 1), dtype=float)
        paths[:, 0] = start
        paths[:, 1:] = start + cumsum

        terminal_positions = paths[:, -1]
        metrics = {
            "mean_terminal": float(np.mean(terminal_positions)),
            "std_terminal": float(np.std(terminal_positions))
        }

        return Result(
            name="random_walk",
            params=call_params,
            metrics=metrics,
            data=paths,
            success=True,
            message=f"随机游走模拟成功：生成了 {n_paths} 条路径，每条路径 {n_steps} 步。"
        )
    except Exception as e:
        return Result(
            name="random_walk",
            params=call_params,
            metrics={},
            data=None,
            success=False,
            message=f"随机游走模拟失败：{type(e).__name__} - {str(e)}"
        )


def markov_simulate(
    transition: Any,
    initial: Any,
    steps: int | np.integer = 100,
    n_paths: int | np.integer = 1,
    random_state: int | None = None,
) -> Any:
    """马尔可夫链模拟：从 `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", "terminal_distribution"}` —— 前者是终点状态的均值，
    后者是 `{状态: 出现频率}` 的 dict（按终点统计，只含**实际出现过**的状态，没出现的状态没有 key）。

    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.integer))` 区分"标量状态 / 状态数组"，
    `np.int64` 已收，但 `bool` 仍会被当成状态 1（`isinstance(True, int)` 为真）；
    且数组分支"按**状态数**校验长度"与真正写进去时"要按**路径数**"仍互相矛盾 —— 只有 `n_paths == s` 才跑得通。
    状态越界（如 `initial=5`、`s=2`）也要等第一轮索引才炸。口径定了这段随之改写。

    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 np.any(arr < -1e-12):
        raise ValueError("transition 矩阵的元素必须是非负的")
    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,np.integer)):
        current_states = np.full(n_paths, initial)
        initial_param = int(initial)
    else:
        current_states, _ = as_array(initial, ndim=1)
        if len(current_states) != s:
            raise TypeError(f"initial 的长度必须是 {s}，而你传入的是 {len(current_states)}")
        initial_param = current_states.tolist()

    call_params = {
        "steps": steps,
        "n_paths": n_paths,
        "initial": initial_param,
        "random_state": random_state
    }

    try:
        paths = np.zeros((n_paths, steps + 1), dtype=int)
        paths[:, 0] = current_states
        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]
        unique, counts = np.unique(terminal_positions, return_counts=True)
        terminal_dist = {int(state): float(count) / n_paths for state, count in zip(unique, counts, strict=True)}

        metrics = {
            "mean_terminal": float(np.mean(terminal_positions)),
            "terminal_distribution": terminal_dist
        }
        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 | np.integer = 10000, tol: float = 1e-10) -> Any:
    """迭代法求马尔可夫链的平稳分布（`v @ P` 反复推进，直到相邻两次几乎不变）。

    Parameters
    ----------
    transition:
        行随机矩阵（每行和为 1）的二维方阵。
    max_iter:
        迭代上限，正整数，默认 10000。
    tol:
        收敛阈值，正数，默认 1e-10。判据是**相邻两次迭代的最大分量差**
        `max|v_k+1 - v_k| < tol`，不是"到不动点的距离"，见 Notes。

    Returns
    -------
    `Result`：`data` 是一维平稳分布向量（和为 1）；
    `metrics={"converged", "iterations"}`（是否在 `max_iter` 内向容、实际迭代次数）。

    Raises
    ------
    TypeError:
        `transition` 不是二维方阵；`max_iter` 不是整数；`tol` 不是数值。
    ValueError:
        `transition` 含负元素、或某一行和不等于 1.0；
        `max_iter` 不是正整数；`tol` 不是正数（NaN 也算非法）。

    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
    --------
    >>> 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
    """
    arr, _ = as_array(transition)
    if arr.ndim != 2 or arr.shape[0] != arr.shape[1]:
        raise TypeError(f"transition 必须是二维方阵，你传入的形状是 {arr.shape}")
    if np.any(arr < -1e-12):
        raise ValueError("transition 矩阵的元素必须是非负的（概率不能为负）")
    if not np.allclose(arr.sum(axis=1), 1.0):
        raise ValueError("transition 每一行的和必须是 1.0")

    if isinstance(max_iter, bool) or not isinstance(max_iter, (int, np.integer)):
        raise TypeError(f"max_iter 必须是正整数，你传入了 {type(max_iter).__name__}")
    if max_iter <= 0:
        raise ValueError(f"max_iter 必须是正整数，你传入了 {max_iter}")
    if not isinstance(tol, (int, float)) or isinstance(tol, bool):
        raise TypeError(f"tol 必须是浮点数，你传入了 {type(tol).__name__}")
    if tol <= 0 or np.isnan(tol):
        raise ValueError(f"tol 必须是正数，你传入了 {tol}")

    s = arr.shape[0]
    v = np.full(s, 1.0 / s)

    iterations = 0
    converged = False
    call_params = {
        "transition": arr.tolist(),
        "max_iter": max_iter,
        "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=call_params,
            metrics=metrics,
            data=v,
            success=converged,
            message=f"平稳分布求解{'成功' if converged else '失败：未在最大迭代次数内收敛'}，迭代 {iterations} 次。"
        )
    except Exception as e:
        return Result(
            name="markov_steady_state",
            params=call_params,
            metrics={},
            data=None,
            success=False,
            message=f"平稳分布求解异常：{type(e).__name__} - {str(e)}"
        )


def mm1_metrics(arrival_rate: float, service_rate: float) -> Any:
    """M/M/1 排队系统的闭式指标。

    Parameters
    ----------
    arrival_rate, service_rate:
        到达率 λ 与服务率 μ，都收 `int` / `float` / numpy 标量。
        `λ` 允许为 0（空系统：`L = Lq = Wq = 0`，`W` 取极限值 `1/μ`）；
        `λ >= μ` 时系统不稳定、队列长度无界 —— 这种情况返回 `success=False` 并说明原因，而不是返回 `inf`。

    Returns
    -------
    `Result`：`metrics={"rho", "L", "Lq", "W", "Wq"}`
    （系统利用率、系统内平均顾客数、队列长度、平均逗留时间、平均等待时间）；
    `data` 与 `metrics` 同值（闭式解没有样本数组），但**不是同一个对象**。
    不稳定时 `success=False`、`metrics={{"rho"}}`（只报利用率）、`data=None`。

    Raises
    ------
    TypeError:
        `arrival_rate` / `service_rate` 不是数值（`bool` 也判非法）。
    ValueError:
        `arrival_rate` 是负数；`service_rate` 不是正数。

    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
    """
    if not isinstance(arrival_rate, (int, float, np.number)) or isinstance(arrival_rate, bool):
        raise TypeError(f"arrival_rate 必须是浮点数，你传入了 {type(arrival_rate).__name__}")
    if arrival_rate < 0:
        raise ValueError(f"arrival_rate 必须是非负数，你传入了 {arrival_rate}")
    if not isinstance(service_rate, (int, float, np.number)) or isinstance(service_rate, bool):
        raise TypeError(f"service_rate 必须是浮点数，你传入了 {type(service_rate).__name__}")
    if service_rate <= 0:
        raise ValueError(f"service_rate 必须是正数，你传入了 {service_rate}")

    call_params = {"arrival_rate": float(arrival_rate), "service_rate": float(service_rate)}

    if arrival_rate >= service_rate:
        rho = arrival_rate / service_rate
        return Result(
            name="mm1_metrics",
            params=call_params,
            metrics={"rho": float(rho)},
            data=None,
            success=False,
            message=f"系统不稳定：λ >= μ，队列无界。当前利用率 rho={rho:.4f}"
        )
    try:
        rho = arrival_rate / service_rate
        L = rho / (1.0 - rho)                   # 系统内平均顾客数
        Lq = rho ** 2 / (1.0 - rho)             # 队列中平均等待顾客数
        W = 1.0 / (service_rate - arrival_rate) # 平均逗留时间
        Wq = W - 1.0 / service_rate             # 平均等待时间

        metrics = {
            "rho": float(rho),
            "L": float(L),
            "Lq": float(Lq),
            "W": float(W),
            "Wq": float(Wq)
        }

        return Result(
            name="mm1_metrics",
            params=call_params,
            metrics=metrics,
            data=metrics.copy(),
            success=True,
            message="M/M/1 排队系统指标计算成功。"
        )

    except Exception as e:
        return Result(
            name="mm1_metrics",
            params=call_params,
            metrics={},
            data=None,
            success=False,
            message=f"计算失败：{type(e).__name__} - {str(e)}"
        )