aps-agent/server/aps_domain/sensitivity.py

# ============================================================
# 主控参数敏感性分析（moduleId: domain-sensitivity, SC-06，可重生 ✅）
# 沙盒 one-at-a-time 扰动 + 固定种子蒙特卡洛；永不写主干
# ============================================================
from __future__ import annotations

import copy
import uuid
from collections.abc import Callable
from datetime import date
from typing import Any

from server.aps_domain.params import get_schedule_params
from server.aps_domain.robustness import run_monte_carlo
from server.contracts import UIBlock
from server.engines import get_engine
from server.engines.base import EngineParams
from server.timeutil import add_minutes, fmt_date, today0

World = dict[str, Any]

SOBOL_METRICS = frozenset({"tardiness", "conflicts", "utilization", "changeoverMin"})
_SOBOL_FACTORS: tuple[dict[str, Any], ...] = (
    {"id": "planningHorizonDays", "label": "展望期(天)", "low": 3.0, "high": 60.0},
    {"id": "vipWeight", "label": "VIP 等级权重", "low": 1.0, "high": 10.0},
    {"id": "deliveryBufferRatio", "label": "交期缓冲比", "low": 0.85, "high": 1.0},
    {"id": "freezeWindowHours", "label": "冻结窗口(时)", "low": 0.0, "high": 48.0},
    {"id": "lineEfficiency", "label": "产线效率系数", "low": 0.8, "high": 1.2},
)


def _sandbox_counter():
    counters: dict[str, int] = {}

    def next_id(kind: str) -> int:
        counters[kind] = counters.get(kind, 0) + 2000000
        return counters[kind]

    return next_id


def _kpi_from_result(result, sandbox: World) -> dict[str, float]:
    ver = sandbox["scheduleVersions"][-1] if sandbox.get("scheduleVersions") else {}
    return {
        "tardiness": round(float(result.totalTardiness), 2),
        "conflicts": float(result.conflictCount),
        "utilization": round(float(result.avgUtilization), 4),
        "changeoverMin": float(ver.get("totalChangeoverMin") or 0),
        "poCount": float(result.poCount),
        "woCount": float(result.woCount),
    }


def _run_sandbox(
    world: World,
    *,
    strategy: str,
    horizon: int,
    start: str,
    vip_weight: float | None = None,
    delivery_buffer: float | None = None,
    freeze_hours: float | None = None,
    efficiency_scale: float | None = None,
) -> dict[str, float]:
    sandbox = copy.deepcopy(world)
    sandbox["scheduleVersions"] = []
    sandbox["productionOrders"] = [p for p in sandbox.get("productionOrders", []) if p.get("status") == "PUBLISHED"]
    sandbox["workOrders"] = [w for w in sandbox.get("workOrders", []) if w.get("status") not in ("PENDING", "DRAFT", None)]
    # 清空草稿工单/PO，避免干扰；简化：只保留非本沙盒相关——种子通常无已发布，直接清空排产结果
    sandbox["productionOrders"] = []
    sandbox["workOrders"] = []
    sandbox["conflicts"] = []

    sp = sandbox.setdefault("scheduleParams", {})
    if vip_weight is not None:
        levels = dict(sp.get("customerLevelWeights") or {})
        levels["VIP"] = float(vip_weight)
        sp["customerLevelWeights"] = levels
    if efficiency_scale is not None and efficiency_scale > 0:
        for ln in sandbox.get("lines") or []:
            base = float(ln.get("efficiencyFactor") or 1.0)
            ln["efficiencyFactor"] = round(base * float(efficiency_scale), 4)

    params = EngineParams(
        orderIds=[],
        engineType="RULE",
        strategyTemplate=strategy,
        planningHorizonDays=horizon,
        startDate=start,
        deliveryBufferRatio=delivery_buffer,
        freezeWindowHours=freeze_hours,
        name=f"sensitivity-{uuid.uuid4().hex[:6]}",
        constraints={
            "materialKit": False, "equipment": True, "personnel": True, "changeover": True,
            "capacity": True, "dueDate": True, "tooling": True, "freeze": True,
        },
    )
    result = get_engine("RULE").solve(sandbox, params, _sandbox_counter())
    return _kpi_from_result(result, sandbox)


def _resolve_sobol_start_date(world: World, start_date: str | None) -> tuple[str, str]:
    source = "request" if start_date is not None else "world.businessDate"
    raw = start_date if start_date is not None else world.get("businessDate")
    if not isinstance(raw, str) or not raw.strip():
        raise ValueError("startDate 必填；仅可省略于 world.businessDate 已设置时")
    value = raw.strip()
    try:
        parsed = date.fromisoformat(value)
    except ValueError as exc:
        raise ValueError("startDate 须为 YYYY-MM-DD") from exc
    if parsed.isoformat() != value:
        raise ValueError("startDate 须为 YYYY-MM-DD")
    return value, source


def _validate_sobol_options(*, base_samples: int, seed: int, metric: str) -> None:
    if type(base_samples) is not int or not 8 <= base_samples <= 64 or base_samples & (base_samples - 1):
        raise ValueError("baseSamples 须为 8、16、32 或 64")
    if type(seed) is not int or not 0 <= seed <= 2**32 - 1:
        raise ValueError("seed 须为 0~4294967295 的整数")
    if metric not in SOBOL_METRICS:
        raise ValueError(f"metric 须为 {' / '.join(sorted(SOBOL_METRICS))}")


def _sobol_factor_values(unit_row: Any) -> dict[str, float]:
    values: dict[str, float] = {}
    for factor, unit_value in zip(_SOBOL_FACTORS, unit_row, strict=True):
        value = factor["low"] + (factor["high"] - factor["low"]) * float(unit_value)
        values[factor["id"]] = float(round(value, 8))
    values["planningHorizonDays"] = float(round(values["planningHorizonDays"]))
    return values


def run_sobol_sensitivity(
    world: World,
    *,
    strategy: str = "COMPREHENSIVE",
    base_samples: int = 32,
    seed: int = 20260818,
    metric: str = "tardiness",
    start_date: str | None = None,
    cancel_check: Callable[[], None] | None = None,
) -> dict[str, Any]:
    """Run scrambled Sobol sampling with Jansen first/total-order estimators."""
    _validate_sobol_options(base_samples=base_samples, seed=seed, metric=metric)
    resolved_start, start_source = _resolve_sobol_start_date(world, start_date)
    from server.aps_domain.explore_boundary import run_explore
    return run_explore(world, lambda sandbox: _sobol_impl(
        sandbox,
        strategy=strategy,
        base_samples=base_samples,
        seed=seed,
        metric=metric,
        start_date=resolved_start,
        start_date_source=start_source,
        cancel_check=cancel_check,
    ))


def _sobol_impl(
    world: World,
    *,
    strategy: str,
    base_samples: int,
    seed: int,
    metric: str,
    start_date: str,
    start_date_source: str,
    cancel_check: Callable[[], None] | None,
) -> dict[str, Any]:
    from scipy.stats import qmc

    dimension = len(_SOBOL_FACTORS)
    sampler = qmc.Sobol(d=dimension * 2, scramble=True, seed=seed)
    paired = sampler.random_base2(m=base_samples.bit_length() - 1)
    matrix_a = paired[:, :dimension]
    matrix_b = paired[:, dimension:]
    sp = get_schedule_params(world)

    def evaluate(unit_row: Any) -> float:
        if cancel_check is not None:
            cancel_check()
        values = _sobol_factor_values(unit_row)
        kpis = _run_sandbox(
            world,
            strategy=strategy,
            horizon=int(values["planningHorizonDays"]),
            start=start_date,
            vip_weight=values["vipWeight"],
            delivery_buffer=values["deliveryBufferRatio"],
            freeze_hours=values["freezeWindowHours"],
            efficiency_scale=values["lineEfficiency"],
        )
        return float(kpis[metric])

    values_a = [evaluate(row) for row in matrix_a]
    values_b = [evaluate(row) for row in matrix_b]
    hybrid_values: list[list[float]] = []
    for factor_index in range(dimension):
        hybrid = matrix_a.copy()
        hybrid[:, factor_index] = matrix_b[:, factor_index]
        hybrid_values.append([evaluate(row) for row in hybrid])

    all_base_values = values_a + values_b
    mean = sum(all_base_values) / len(all_base_values)
    variance = sum((value - mean) ** 2 for value in all_base_values) / len(all_base_values)
    degenerate = variance <= 1e-12
    factors: list[dict[str, Any]] = []
    for factor, hybrid in zip(_SOBOL_FACTORS, hybrid_values, strict=True):
        if degenerate:
            first_order = total_order = interaction = None
        else:
            total_order = sum(
                (a_value - hybrid_value) ** 2
                for a_value, hybrid_value in zip(values_a, hybrid, strict=True)
            ) / (2.0 * base_samples * variance)
            first_order = 1.0 - sum(
                (b_value - hybrid_value) ** 2
                for b_value, hybrid_value in zip(values_b, hybrid, strict=True)
            ) / (2.0 * base_samples * variance)
            interaction = total_order - first_order
        factors.append({
            "factorId": factor["id"],
            "label": factor["label"],
            "range": {"low": factor["low"], "high": factor["high"]},
            "firstOrder": round(first_order, 6) if first_order is not None else None,
            "totalOrder": round(total_order, 6) if total_order is not None else None,
            "interaction": round(interaction, 6) if interaction is not None else None,
            "rank": None,
        })
    if not degenerate:
        factors.sort(key=lambda row: (-row["totalOrder"], row["factorId"]))
        for rank, factor in enumerate(factors, start=1):
            factor["rank"] = rank

    return {
        "method": "sobol-jansen",
        "sampler": "scipy.stats.qmc.Sobol",
        "estimator": "Jansen first/total-order",
        "scrambled": True,
        "strategy": strategy,
        "seed": seed,
        "startDate": start_date,
        "startDateSource": start_date_source,
        "baseSamples": base_samples,
        "factorCount": dimension,
        "evaluations": base_samples * (dimension + 2),
        "metric": metric,
        "variance": round(variance, 12),
        "isDegenerate": degenerate,
        "scheduleParamsSnapshot": {
            "planningHorizonDays": int(
                sp["planningHorizonDays"] if sp.get("planningHorizonDays") is not None else 14
            ),
            "vipWeight": float(
                (sp.get("customerLevelWeights") or {}).get("VIP")
                if (sp.get("customerLevelWeights") or {}).get("VIP") is not None
                else 3.0
            ),
            "deliveryBufferRatio": float(
                sp["deliveryBufferRatio"] if sp.get("deliveryBufferRatio") is not None else 0.95
            ),
            "freezeWindowHours": float(
                sp["freezeWindowHours"] if sp.get("freezeWindowHours") is not None else 24
            ),
        },
        "factors": factors,
        "hint": (
            "指标无方差，Sobol 指数不可识别。"
            if degenerate
            else "Sobol 结论来自 Explore 沙盒；一阶指数衡量单因子贡献，总效应指数包含交互贡献。"
        ),
    }


def run_sensitivity(
    world: World,
    *,
    strategy: str = "COMPREHENSIVE",
) -> dict[str, Any]:
    """
    SC-06：Tornado 单因子排序 + 固定种子蒙特卡洛鲁棒性。

    经统一 Explore 通道执行（矩阵 55 行）：fn 只拿深拷贝沙盒，主干永不外泄写引用；
    _run_sandbox 内部再做子沙盒扰动。
    """
    from server.aps_domain.explore_boundary import run_explore
    return run_explore(world, lambda sandbox: _sensitivity_impl(
        sandbox, strategy=strategy,
    ))


def _sensitivity_impl(world: World, *, strategy: str) -> dict[str, Any]:
    """Tornado/蒙特卡洛核心实现（在统一通道沙盒内运行）。"""
    sp = get_schedule_params(world)
    base_horizon = int(sp.get("planningHorizonDays") or 14)
    base_vip = float((sp.get("customerLevelWeights") or {}).get("VIP") or 3.0)
    base_buffer = float(sp.get("deliveryBufferRatio") or 0.95)
    base_freeze = float(sp.get("freezeWindowHours") or 24)
    start = fmt_date(add_minutes(today0(), 24 * 60))

    baseline = _run_sandbox(
        world, strategy=strategy, horizon=base_horizon, start=start,
        vip_weight=base_vip, delivery_buffer=base_buffer, freeze_hours=base_freeze,
        efficiency_scale=1.0,
    )

    # 每因子：低/高 两档（相对基线）
    factors: list[dict[str, Any]] = [
        {
            "id": "planningHorizonDays", "label": "展望期(天)",
            "baseline": base_horizon,
            "low": {"label": f"{max(3, base_horizon // 2)} 天", "kwargs": {"horizon": max(3, base_horizon // 2)}},
            "high": {"label": f"{min(60, base_horizon * 2)} 天", "kwargs": {"horizon": min(60, base_horizon * 2)}},
        },
        {
            "id": "vipWeight", "label": "VIP 等级权重",
            "baseline": base_vip,
            "low": {"label": "VIP=1", "kwargs": {"vip_weight": 1.0}},
            "high": {"label": "VIP=10", "kwargs": {"vip_weight": 10.0}},
        },
        {
            "id": "deliveryBufferRatio", "label": "交期缓冲比",
            "baseline": base_buffer,
            "low": {"label": "缓冲 0.85", "kwargs": {"delivery_buffer": 0.85}},
            "high": {"label": "缓冲 1.00", "kwargs": {"delivery_buffer": 1.0}},
        },
        {
            "id": "freezeWindowHours", "label": "冻结窗口(时)",
            "baseline": base_freeze,
            "low": {"label": "冻结 0h", "kwargs": {"freeze_hours": 0.0}},
            "high": {"label": "冻结 48h", "kwargs": {"freeze_hours": 48.0}},
        },
        {
            "id": "lineEfficiency", "label": "产线效率系数",
            "baseline": 1.0,
            "low": {"label": "效率×0.8", "kwargs": {"efficiency_scale": 0.8}},
            "high": {"label": "效率×1.2", "kwargs": {"efficiency_scale": 1.2}},
        },
    ]

    rows: list[dict[str, Any]] = []
    for fac in factors:
        common = {
            "strategy": strategy, "horizon": base_horizon, "start": start,
            "vip_weight": base_vip, "delivery_buffer": base_buffer,
            "freeze_hours": base_freeze, "efficiency_scale": 1.0,
        }
        low_kw = {**common, **fac["low"]["kwargs"]}
        high_kw = {**common, **fac["high"]["kwargs"]}
        # kwargs 用 horizon 键
        def _call(kw: dict) -> dict[str, float]:
            return _run_sandbox(
                world,
                strategy=kw["strategy"],
                horizon=int(kw["horizon"]),
                start=kw["start"],
                vip_weight=kw.get("vip_weight"),
                delivery_buffer=kw.get("delivery_buffer"),
                freeze_hours=kw.get("freeze_hours"),
                efficiency_scale=kw.get("efficiency_scale"),
            )

        low_kpi = _call(low_kw)
        high_kpi = _call(high_kw)
        d_low = round(low_kpi["tardiness"] - baseline["tardiness"], 2)
        d_high = round(high_kpi["tardiness"] - baseline["tardiness"], 2)
        rows.append({
            "factorId": fac["id"],
            "label": fac["label"],
            "baselineValue": fac["baseline"],
            "lowLabel": fac["low"]["label"],
            "highLabel": fac["high"]["label"],
            "low": {
                "tardiness": low_kpi["tardiness"],
                "conflicts": low_kpi["conflicts"],
                "utilization": low_kpi["utilization"],
                "deltaTardiness": d_low,
                "deltaConflicts": round(low_kpi["conflicts"] - baseline["conflicts"], 1),
            },
            "high": {
                "tardiness": high_kpi["tardiness"],
                "conflicts": high_kpi["conflicts"],
                "utilization": high_kpi["utilization"],
                "deltaTardiness": d_high,
                "deltaConflicts": round(high_kpi["conflicts"] - baseline["conflicts"], 1),
            },
            "swing": round(abs(d_low) + abs(d_high), 2),
            "impactMetric": "tardiness",
        })

    rows.sort(key=lambda r: r["swing"], reverse=True)

    monte_carlo = run_monte_carlo(
        world,
        strategy=strategy,
        engine_type="RULE",
        baseline_kpi=baseline,
        planning_horizon_days=base_horizon,
        start_date=start,
        delivery_buffer_ratio=base_buffer,
        freeze_window_hours=base_freeze,
        constraints={
            "materialKit": False, "equipment": True, "personnel": True, "changeover": True,
            "capacity": True, "dueDate": True, "tooling": True, "freeze": True,
        },
        reset_schedule_products=True,
    )

    md_lines = [
        "# 敏感性分析（Tornado）",
        "",
        f"- 策略基线：`{strategy}`",
        f"- 基线 KPI：延期 **{baseline['tardiness']}** h · 冲突 **{int(baseline['conflicts'])}** · "
        f"利用率 **{round(baseline['utilization'] * 100, 1)}%** · 换型 **{baseline['changeoverMin']}** 分",
        "- 说明：沙盒 one-at-a-time，不写主干；超时秒数对 RULE 无意义，未纳入。",
        "",
        "| 因子 | 低档 Δ延期 | 高档 Δ延期 | 摆幅 |",
        "| --- | ---: | ---: | ---: |",
    ]
    for r in rows:
        md_lines.append(
            f"| {r['label']}（{r['lowLabel']} / {r['highLabel']}） | "
            f"{r['low']['deltaTardiness']:+} | {r['high']['deltaTardiness']:+} | {r['swing']} |"
        )
    md_lines.append("")
    md_lines.append("摆幅 = |低档Δ延期| + |高档Δ延期|；越大表示该参数对延期越敏感。")
    md_lines.extend([
        "",
        "## 固定种子蒙特卡洛鲁棒性",
        "",
        f"- 种子：`{monte_carlo['seed']}` · 样本：**{monte_carlo['trials']}** 次",
        f"- 鲁棒性：**{round(monte_carlo['robustness'] * 100, 1)}%**",
        f"- 延期分布：P50 **{monte_carlo['distribution']['tardinessP50']}h** · "
        f"P90 **{monte_carlo['distribution']['tardinessP90']}h**",
    ])

    return {
        "strategy": strategy,
        "baseline": baseline,
        "rows": rows,
        "monteCarlo": monte_carlo,
        "markdown": "\n".join(md_lines),
        "hint": "敏感性只跑沙盒。若要固化参数，请到设置 → 排产参数 改权重/展望期并确认（P2）。",
    }


def sensitivity_as_block(report: dict[str, Any]) -> UIBlock:
    """产出 report UI 块（可预览/下载）。"""
    return UIBlock(
        blockId=f"sensitivity-{uuid.uuid4().hex[:8]}",
        type="report",
        props={
            "reportId": uuid.uuid4().hex[:10],
            "title": f"敏感性分析 Tornado（{report.get('strategy')}）",
            "markdown": report.get("markdown") or "",
            "reportType": "sensitivity-tornado",
            "baseline": report.get("baseline"),
            "rows": report.get("rows"),
            "monteCarlo": report.get("monteCarlo"),
        },
    )