#!/usr/bin/env python3
"""MSC-P-034 deterministic Monte Carlo reference. MIT; standard library only."""

from __future__ import annotations

import csv
import math
import random
import statistics
from pathlib import Path


DATA = Path(__file__).resolve().parents[1] / "datasets" / "msc-010-monte-carlo-inputs.csv"
SCENARIO = "base"
SEED = 20260822
DRAWS = 20_000


def load_inputs() -> dict[str, dict[str, str]]:
    rows = [row for row in csv.DictReader(DATA.open(encoding="utf-8", newline="")) if row["scenario"] == SCENARIO]
    values = {row["input"]: row for row in rows}
    required = {"baseline_price", "baseline_volume", "variable_cost", "elasticity", "price_change", "fixed_cost_change"}
    if set(values) != required:
        raise ValueError("base scenario does not contain the six required inputs")
    return values


def parameter(row: dict[str, str], index: int) -> float:
    return float(row[f"param_{index}"])


def quantile(sorted_values: list[float], probability: float) -> float:
    position = (len(sorted_values) - 1) * probability
    lower = math.floor(position)
    upper = math.ceil(position)
    if lower == upper:
        return sorted_values[lower]
    return sorted_values[lower] + (position - lower) * (sorted_values[upper] - sorted_values[lower])


def main() -> None:
    inputs = load_inputs()
    rng = random.Random(SEED)
    p0 = parameter(inputs["baseline_price"], 1)
    price_change = parameter(inputs["price_change"], 1)
    volume_mean, volume_sd = parameter(inputs["baseline_volume"], 1), parameter(inputs["baseline_volume"], 2)
    elasticity_mean, elasticity_sd = parameter(inputs["elasticity"], 1), parameter(inputs["elasticity"], 2)
    rho = float(inputs["baseline_volume"]["rho"])
    vc_low, vc_mode, vc_high = (parameter(inputs["variable_cost"], i) for i in (1, 2, 3))
    fc_low, fc_mode, fc_high = (parameter(inputs["fixed_cost_change"], i) for i in (1, 2, 3))
    outcomes: list[float] = []
    for _ in range(DRAWS):
        z_volume = rng.gauss(0, 1)
        z_independent = rng.gauss(0, 1)
        z_elasticity = rho * z_volume + math.sqrt(1 - rho * rho) * z_independent
        q0 = max(1.0, volume_mean + volume_sd * z_volume)
        elasticity = elasticity_mean + elasticity_sd * z_elasticity
        vc = rng.triangular(vc_low, vc_high, vc_mode)
        delta_fc = rng.triangular(fc_low, fc_high, fc_mode)
        p1 = p0 * (1 + price_change)
        q1 = q0 * (1 + price_change) ** elasticity
        outcomes.append((p1 - vc) * q1 - (p0 - vc) * q0 - delta_fc)

    outcomes.sort()
    negatives = sum(value < 0 for value in outcomes)
    risk = negatives / DRAWS
    mc_se = math.sqrt(risk * (1 - risk) / DRAWS)
    print(f"scenario={SCENARIO}")
    print(f"seed={SEED}")
    print(f"draws={DRAWS}")
    print(f"mean_incremental_profit={statistics.fmean(outcomes):.6f}")
    print(f"median_incremental_profit={statistics.median(outcomes):.6f}")
    print(f"p05={quantile(outcomes, 0.05):.6f}")
    print(f"p95={quantile(outcomes, 0.95):.6f}")
    print(f"probability_negative={risk:.6f}")
    print(f"mc_standard_error={mc_se:.6f}")


if __name__ == "__main__":
    main()