///|
/// Bounds for deterministic uncertainty sampling, expressed in permille.
pub(all) struct UncertaintyConfig {
  solar_min_permille : Int
  solar_max_permille : Int
  base_load_min_permille : Int
  base_load_max_permille : Int
  tariff_min_permille : Int
  tariff_max_permille : Int
  outage_start_jitter_slots : Int
  outage_duration_jitter_slots : Int
  temporary_spike_probability_permille : Int
  temporary_spike_min_w : Int
  temporary_spike_max_w : Int
  correlated_signal_permille : Int
} derive(Debug, Eq, ToJson, FromJson)

///|
pub fn UncertaintyConfig::typical() -> UncertaintyConfig {
  {
    solar_min_permille: 700,
    solar_max_permille: 1150,
    base_load_min_permille: 900,
    base_load_max_permille: 1250,
    tariff_min_permille: 950,
    tariff_max_permille: 1100,
    outage_start_jitter_slots: 1,
    outage_duration_jitter_slots: 2,
    temporary_spike_probability_permille: 80,
    temporary_spike_min_w: 300,
    temporary_spike_max_w: 1600,
    correlated_signal_permille: 700,
  }
}

///|
pub fn UncertaintyConfig::conservative() -> UncertaintyConfig {
  {
    solar_min_permille: 450,
    solar_max_permille: 1000,
    base_load_min_permille: 950,
    base_load_max_permille: 1500,
    tariff_min_permille: 1000,
    tariff_max_permille: 1300,
    outage_start_jitter_slots: 2,
    outage_duration_jitter_slots: 4,
    temporary_spike_probability_permille: 160,
    temporary_spike_min_w: 500,
    temporary_spike_max_w: 2500,
    correlated_signal_permille: 800,
  }
}

///|
pub fn UncertaintyConfig::none() -> UncertaintyConfig {
  {
    solar_min_permille: 1000,
    solar_max_permille: 1000,
    base_load_min_permille: 1000,
    base_load_max_permille: 1000,
    tariff_min_permille: 1000,
    tariff_max_permille: 1000,
    outage_start_jitter_slots: 0,
    outage_duration_jitter_slots: 0,
    temporary_spike_probability_permille: 0,
    temporary_spike_min_w: 0,
    temporary_spike_max_w: 0,
    correlated_signal_permille: 1000,
  }
}

///|
priv struct DeterministicRng {
  mut state : UInt
}

///|
fn DeterministicRng::new(seed : UInt) -> DeterministicRng {
  { state: if seed == 0U { 0x9E3779B9U } else { seed } }
}

///|
fn DeterministicRng::next_u32(self : DeterministicRng) -> UInt {
  let mut x = self.state
  x = x ^ (x << 13)
  x = x ^ (x >> 17)
  x = x ^ (x << 5)
  self.state = x
  x
}

///|
fn DeterministicRng::range(
  self : DeterministicRng,
  lower : Int,
  upper : Int,
) -> Int {
  if upper <= lower {
    return lower
  }
  let width = (upper - lower + 1).reinterpret_as_uint()
  lower + (self.next_u32() % width).reinterpret_as_int()
}

///|
fn DeterministicRng::chance_permille(
  self : DeterministicRng,
  probability : Int,
) -> Bool {
  self.range(0, 999) < probability
}

///|
fn blend_factor(previous : Int, sampled : Int, correlation : Int) -> Int {
  previous * correlation / 1000 + sampled * (1000 - correlation) / 1000
}

///|
fn sampled_series(
  source : IntSeries,
  minimum_permille : Int,
  maximum_permille : Int,
  correlation_permille : Int,
  rng : DeterministicRng,
) -> (IntSeries, Int) {
  let values : Array[Int] = []
  let mut factor = rng.range(minimum_permille, maximum_permille)
  let mut factor_sum = 0
  for value in source.values {
    let sampled = rng.range(minimum_permille, maximum_permille)
    factor = blend_factor(factor, sampled, correlation_permille)
    factor = clamp(factor, minimum_permille, maximum_permille)
    values.push(maximum(0, value * factor / 1000))
    factor_sum = factor_sum + factor
  }
  let average_factor = if source.values.length() == 0 {
    1000
  } else {
    factor_sum / source.values.length()
  }
  ({ ..source, values, }, average_factor)
}

///|
fn sampled_base_load(
  source : IntSeries,
  config : UncertaintyConfig,
  rng : DeterministicRng,
) -> (IntSeries, Int, Int, Int) {
  let (scaled, average_factor) = sampled_series(
    source,
    config.base_load_min_permille,
    config.base_load_max_permille,
    config.correlated_signal_permille,
    rng,
  )
  let values = scaled.copy_values()
  let mut spike_slot = -1
  let mut spike_power_w = 0
  for slot = 0; slot < values.length(); slot = slot + 1 {
    if rng.chance_permille(config.temporary_spike_probability_permille) {
      let spike = rng.range(
        config.temporary_spike_min_w,
        config.temporary_spike_max_w,
      )
      values[slot] = values[slot] + spike
      if spike > spike_power_w {
        spike_slot = slot
        spike_power_w = spike
      }
    }
  }
  ({ ..scaled, values, }, average_factor, spike_slot, spike_power_w)
}

///|
fn sampled_outage(
  outage : OutageEvent,
  horizon : Int,
  config : UncertaintyConfig,
  rng : DeterministicRng,
) -> OutageEvent {
  let start_jitter = rng.range(
    -config.outage_start_jitter_slots,
    config.outage_start_jitter_slots,
  )
  let duration_jitter = rng.range(
    -config.outage_duration_jitter_slots,
    config.outage_duration_jitter_slots,
  )
  let original_duration = outage.duration_slots()
  let duration = maximum(1, original_duration + duration_jitter)
  let start = clamp(
    outage.start_slot + start_jitter,
    0,
    maximum(0, horizon - 1),
  )
  let end = minimum(horizon, start + duration)
  { ..outage, start_slot: start, end_slot: end }
}

///|
pub(all) struct ScenarioSample {
  index : Int
  seed : UInt
  solar_factor_permille : Int
  base_load_factor_permille : Int
  tariff_factor_permille : Int
  largest_spike_slot : Int
  largest_spike_w : Int
  outage_start_slots : Array[Int]
  outage_duration_slots : Array[Int]
  input : PlanningInput
} derive(Debug, Eq, ToJson, FromJson)

///|
/// Generate a single reproducible uncertainty sample.
pub fn sample_scenario(
  input : PlanningInput,
  config : UncertaintyConfig,
  index : Int,
) -> ScenarioSample {
  let seed = input.random_seed ^
    ((index + 1).reinterpret_as_uint() * 0x9E3779B9U)
  let rng = DeterministicRng::new(seed)
  let (solar, solar_factor) = sampled_series(
    input.solar_w,
    config.solar_min_permille,
    config.solar_max_permille,
    config.correlated_signal_permille,
    rng,
  )
  let (base_load, base_factor, spike_slot, spike_w) = sampled_base_load(
    input.base_load_w,
    config,
    rng,
  )
  let (tariff, tariff_factor) = sampled_series(
    input.tariff_micro_per_kwh,
    config.tariff_min_permille,
    config.tariff_max_permille,
    config.correlated_signal_permille,
    rng,
  )
  let outages = input.outages.map(outage => {
    sampled_outage(outage, input.horizon_slots, config, rng)
  })
  let outage_starts = outages.map(outage => outage.start_slot)
  let outage_durations = outages.map(outage => outage.duration_slots())
  let sampled_input = {
    ..input,
    title: input.title + " sample " + index.to_string(),
    solar_w: solar,
    base_load_w: base_load,
    tariff_micro_per_kwh: tariff,
    outages,
    random_seed: seed,
  }
  {
    index,
    seed,
    solar_factor_permille: solar_factor,
    base_load_factor_permille: base_factor,
    tariff_factor_permille: tariff_factor,
    largest_spike_slot: spike_slot,
    largest_spike_w: spike_w,
    outage_start_slots: outage_starts,
    outage_duration_slots: outage_durations,
    input: sampled_input,
  }
}

///|
pub(all) struct SimulationRun {
  index : Int
  seed : UInt
  solar_factor_permille : Int
  base_load_factor_permille : Int
  tariff_factor_permille : Int
  largest_spike_slot : Int
  largest_spike_w : Int
  outage_duration_slots : Array[Int]
  status : PlanStatus
  metrics : PlanMetrics
  scheduled_tasks : Int
  explanation_count : Int
} derive(Debug, Eq, ToJson, FromJson)

///|
pub fn SimulationRun::from_sample(
  sample : ScenarioSample,
  result : PlanResult,
) -> SimulationRun {
  {
    index: sample.index,
    seed: sample.seed,
    solar_factor_permille: sample.solar_factor_permille,
    base_load_factor_permille: sample.base_load_factor_permille,
    tariff_factor_permille: sample.tariff_factor_permille,
    largest_spike_slot: sample.largest_spike_slot,
    largest_spike_w: sample.largest_spike_w,
    outage_duration_slots: sample.outage_duration_slots,
    status: result.status,
    metrics: result.metrics,
    scheduled_tasks: result.schedule.length(),
    explanation_count: result.explanations.length(),
  }
}

///|
pub(all) enum RiskBand {
  LowRisk
  ModerateRisk
  HighRisk
  CriticalRisk
} derive(Debug, Eq, ToJson, FromJson)

///|
pub fn RiskBand::label(self : RiskBand) -> String {
  match self {
    LowRisk => "low"
    ModerateRisk => "moderate"
    HighRisk => "high"
    CriticalRisk => "critical"
  }
}

///|
pub(all) struct DistributionSummary {
  minimum : Int
  p10 : Int
  median : Int
  p90 : Int
  p95 : Int
  maximum : Int
  mean : Int
} derive(Debug, Eq, ToJson, FromJson)

///|
fn percentile_index(length : Int, percentile : Int) -> Int {
  if length <= 1 {
    0
  } else {
    clamp((length - 1) * percentile / 100, 0, length - 1)
  }
}

///|
pub fn summarize_distribution(values : Array[Int]) -> DistributionSummary {
  if values.length() == 0 {
    return {
      minimum: 0,
      p10: 0,
      median: 0,
      p90: 0,
      p95: 0,
      maximum: 0,
      mean: 0,
    }
  }
  let sorted = values.map(value => value)
  sorted.sort()
  let mut total : Int64 = 0L
  for value in sorted {
    total = total + value.to_int64()
  }
  {
    minimum: sorted[0],
    p10: sorted[percentile_index(sorted.length(), 10)],
    median: sorted[percentile_index(sorted.length(), 50)],
    p90: sorted[percentile_index(sorted.length(), 90)],
    p95: sorted[percentile_index(sorted.length(), 95)],
    maximum: sorted[sorted.length() - 1],
    mean: (total / sorted.length().to_int64()).to_int(),
  }
}

///|
pub(all) struct SimulationSummary {
  title : String
  requested_runs : Int
  completed_runs : Int
  feasible_runs : Int
  infeasible_runs : Int
  risk_band : RiskBand
  cost_micro : DistributionSummary
  carbon_g : DistributionSummary
  unserved_energy_wh : DistributionSummary
  critical_unserved_wh : DistributionSummary
  peak_grid_w : DistributionSummary
  resilience_permille : DistributionSummary
  worst_run_index : Int
  best_run_index : Int
  runs : Array[SimulationRun]
  recommendations : Array[String]
} derive(Debug, Eq, ToJson, FromJson)

///|
pub fn SimulationSummary::feasibility_permille(self : SimulationSummary) -> Int {
  if self.completed_runs == 0 {
    0
  } else {
    self.feasible_runs * 1000 / self.completed_runs
  }
}

///|
pub fn SimulationSummary::to_json_string(self : SimulationSummary) -> String {
  self.to_json().stringify(indent=2)
}

///|
fn determine_risk_band(
  completed : Int,
  infeasible : Int,
  critical_p95_wh : Int,
  resilience_p10 : Int,
) -> RiskBand {
  if completed == 0 || infeasible * 4 >= completed || critical_p95_wh > 0 {
    CriticalRisk
  } else if infeasible * 10 >= completed || resilience_p10 < 900 {
    HighRisk
  } else if infeasible > 0 || resilience_p10 < 980 {
    ModerateRisk
  } else {
    LowRisk
  }
}

///|
fn simulation_recommendations(
  input : PlanningInput,
  risk : RiskBand,
  cost : DistributionSummary,
  unserved : DistributionSummary,
  resilience : DistributionSummary,
) -> Array[String] {
  let result : Array[String] = []
  match risk {
    LowRisk =>
      result.push(
        "The plan remains feasible under sampled uncertainty; keep monitoring forecasts.",
      )
    ModerateRisk =>
      result.push(
        "Keep additional battery reserve before the most uncertain hours.",
      )
    HighRisk =>
      result.push(
        "Move noncritical tasks outside outage windows or increase flexible capacity.",
      )
    CriticalRisk =>
      result.push(
        "Critical demand is at risk; add backup capacity or reduce declared critical load.",
      )
  }
  if cost.p95 > cost.median * 12 / 10 {
    result.push(
      "The upper cost tail is material; use tariff alerts or a stricter cost budget.",
    )
  }
  if unserved.p95 > 0 {
    result.push(
      "At least five percent of scenarios shed load; review task priorities and grid limit.",
    )
  }
  if resilience.p10 < 950 {
    result.push(
      "Ten-percent-tail resilience is below 95%; reserve more stored energy.",
    )
  }
  if input.battery is None {
    result.push(
      "No battery is configured; a small storage system can improve outage coverage.",
    )
  }
  result
}

///|
/// Run reproducible scenario analysis. Each sample is derived from the input seed.
pub fn simulate(
  input : PlanningInput,
  runs? : Int = 100,
  uncertainty? : UncertaintyConfig = UncertaintyConfig::typical(),
  solver_config? : SolverConfig = SolverConfig::default(),
) -> SimulationSummary {
  let count = clamp(runs, 1, 10000)
  let results : Array[SimulationRun] = []
  let costs : Array[Int] = []
  let carbon : Array[Int] = []
  let unserved : Array[Int] = []
  let critical_unserved : Array[Int] = []
  let peaks : Array[Int] = []
  let resilience : Array[Int] = []
  let mut feasible = 0
  let mut infeasible = 0
  let mut worst_index = 0
  let mut best_index = 0
  let mut worst_score : Int64 = -1L
  let mut best_score : Int64 = 0x7FFFFFFFFFFFFFFFL
  for index = 0; index < count; index = index + 1 {
    let sample = sample_scenario(input, uncertainty, index)
    let result = solve(sample.input, config=solver_config)
    let run = SimulationRun::from_sample(sample, result)
    results.push(run)
    costs.push(result.metrics.net_cost_micro())
    carbon.push(result.metrics.carbon_g)
    unserved.push(result.metrics.unserved_energy_wh)
    critical_unserved.push(result.metrics.critical_unserved_wh)
    peaks.push(result.metrics.peak_grid_w)
    resilience.push(result.metrics.resilience_permille)
    if result.status == Infeasible {
      infeasible = infeasible + 1
    } else {
      feasible = feasible + 1
    }
    let risk_score = result.metrics.unserved_energy_wh.to_int64() * 100000L +
      result.metrics.critical_unserved_wh.to_int64() * 1000000L +
      result.metrics.net_cost_micro().to_int64()
    if risk_score > worst_score {
      worst_score = risk_score
      worst_index = index
    }
    if risk_score < best_score {
      best_score = risk_score
      best_index = index
    }
  }
  let cost_summary = summarize_distribution(costs)
  let carbon_summary = summarize_distribution(carbon)
  let unserved_summary = summarize_distribution(unserved)
  let critical_summary = summarize_distribution(critical_unserved)
  let peak_summary = summarize_distribution(peaks)
  let resilience_summary = summarize_distribution(resilience)
  let risk = determine_risk_band(
    count,
    infeasible,
    critical_summary.p95,
    resilience_summary.p10,
  )
  {
    title: input.title + " uncertainty simulation",
    requested_runs: runs,
    completed_runs: count,
    feasible_runs: feasible,
    infeasible_runs: infeasible,
    risk_band: risk,
    cost_micro: cost_summary,
    carbon_g: carbon_summary,
    unserved_energy_wh: unserved_summary,
    critical_unserved_wh: critical_summary,
    peak_grid_w: peak_summary,
    resilience_permille: resilience_summary,
    worst_run_index: worst_index,
    best_run_index: best_index,
    runs: results,
    recommendations: simulation_recommendations(
      input, risk, cost_summary, unserved_summary, resilience_summary,
    ),
  }
}

///|
pub(all) struct SensitivityCase {
  id : String
  label : String
  changed_parameter : String
  change_permille : Int
  status : PlanStatus
  metrics : PlanMetrics
  delta_cost_micro : Int
  delta_carbon_g : Int
  delta_unserved_wh : Int
  delta_resilience_permille : Int
} derive(Debug, Eq, ToJson, FromJson)

///|
fn sensitivity_case(
  baseline : PlanResult,
  id : String,
  label : String,
  parameter : String,
  change_permille : Int,
  input : PlanningInput,
) -> SensitivityCase {
  let result = solve(input)
  {
    id,
    label,
    changed_parameter: parameter,
    change_permille,
    status: result.status,
    metrics: result.metrics,
    delta_cost_micro: result.metrics.net_cost_micro() -
    baseline.metrics.net_cost_micro(),
    delta_carbon_g: result.metrics.carbon_g - baseline.metrics.carbon_g,
    delta_unserved_wh: result.metrics.unserved_energy_wh -
    baseline.metrics.unserved_energy_wh,
    delta_resilience_permille: result.metrics.resilience_permille -
    baseline.metrics.resilience_permille,
  }
}

///|
/// Evaluate transparent one-at-a-time changes for key planning assumptions.
pub fn sensitivity_analysis(input : PlanningInput) -> Array[SensitivityCase] {
  let baseline = solve(input)
  let result : Array[SensitivityCase] = []
  for factor in [700, 850, 1150, 1300] {
    result.push(
      sensitivity_case(
        baseline,
        "solar-" + factor.to_string(),
        "Solar forecast " + factor.to_string() + " permille",
        "solar_w",
        factor - 1000,
        { ..input, solar_w: input.solar_w.scale_permille(factor) },
      ),
    )
  }
  for factor in [800, 1200, 1500] {
    result.push(
      sensitivity_case(
        baseline,
        "load-" + factor.to_string(),
        "Base load " + factor.to_string() + " permille",
        "base_load_w",
        factor - 1000,
        { ..input, base_load_w: input.base_load_w.scale_permille(factor) },
      ),
    )
  }
  for factor in [800, 1200, 1500] {
    result.push(
      sensitivity_case(
        baseline,
        "tariff-" + factor.to_string(),
        "Tariff " + factor.to_string() + " permille",
        "tariff_micro_per_kwh",
        factor - 1000,
        {
          ..input,
          tariff_micro_per_kwh: input.tariff_micro_per_kwh.scale_permille(
            factor,
          ),
        },
      ),
    )
  }
  result
}

///|
pub fn most_influential_case(
  cases : Array[SensitivityCase],
) -> SensitivityCase? {
  if cases.length() == 0 {
    return None
  }
  let mut best = cases[0]
  let mut best_magnitude = absolute(best.delta_cost_micro / 1000) +
    absolute(best.delta_carbon_g) +
    absolute(best.delta_unserved_wh) * 100
  for item in cases {
    let magnitude = absolute(item.delta_cost_micro / 1000) +
      absolute(item.delta_carbon_g) +
      absolute(item.delta_unserved_wh) * 100
    if magnitude > best_magnitude {
      best = item
      best_magnitude = magnitude
    }
  }
  Some(best)
}