///|
/// Arithmetic mean of an array.
pub fn mean(values : Array[Double]) -> Double? {
  if values.length() == 0 {
    return None
  }
  let mut total = 0.0
  for value in values {
    total = total + value
  }
  Some(total / values.length().to_double())
}

///|
/// Population standard deviation.
pub fn standard_deviation(values : Array[Double]) -> Double? {
  match mean(values) {
    None => None
    Some(average) => {
      let mut squared_total = 0.0
      for value in values {
        let difference = value - average
        squared_total = squared_total + difference * difference
      }
      Some((squared_total / values.length().to_double()).sqrt())
    }
  }
}

///|
/// Linear-interpolated percentile in the inclusive 0–100 range.
pub fn percentile(values : Array[Double], percent : Double) -> Double? {
  if values.length() == 0 {
    return None
  }
  let ordered = copy_array(values)
  ordered.sort()
  let bounded = clamp_double(percent, 0.0, 100.0)
  if ordered.length() == 1 {
    return Some(ordered[0])
  }
  let position = bounded / 100.0 * (ordered.length() - 1).to_double()
  let lower = position.floor().to_int()
  let upper = position.ceil().to_int()
  if lower == upper {
    Some(ordered[lower])
  } else {
    let fraction = position - lower.to_double()
    Some(ordered[lower] + (ordered[upper] - ordered[lower]) * fraction)
  }
}

///|
/// Median convenience wrapper.
pub fn median(values : Array[Double]) -> Double? {
  percentile(values, 50.0)
}

///|
/// Extract usable temperatures.
fn usable_temperatures(readings : Array[Reading]) -> Array[Double] {
  let values : Array[Double] = []
  for sample in readings {
    if sample.is_usable() {
      values.push(sample.temperature_c)
    }
  }
  values
}

///|
/// Minimum and maximum in one pass.
fn extrema(values : Array[Double]) -> (Double?, Double?) {
  if values.length() == 0 {
    return (None, None)
  }
  let mut minimum = values[0]
  let mut maximum = values[0]
  for value in values {
    if value < minimum {
      minimum = value
    }
    if value > maximum {
      maximum = value
    }
  }
  (Some(minimum), Some(maximum))
}

///|
/// Integrate a value over each interval using the trapezoid rule.
fn time_weighted_temperature(readings : Array[Reading]) -> (Double?, Int64) {
  if readings.length() < 2 {
    return (None, 0L)
  }
  let mut weighted_sum = 0.0
  let mut total_seconds = 0L
  for index = 1; index < readings.length(); index = index + 1 {
    let previous = readings[index - 1]
    let current = readings[index]
    if !previous.is_usable() || !current.is_usable() {
      continue
    }
    let duration = current.timestamp - previous.timestamp
    if duration <= 0L {
      continue
    }
    weighted_sum = weighted_sum +
      (previous.temperature_c + current.temperature_c) /
      2.0 *
      duration.to_double()
    total_seconds = total_seconds + duration
  }
  if total_seconds == 0L {
    (None, 0L)
  } else {
    (Some(weighted_sum / total_seconds.to_double()), total_seconds)
  }
}

///|
/// Mean kinetic temperature using a default activation energy of 83.144 kJ/mol.
/// This is an analytical indicator, not a regulatory disposition decision.
pub fn mean_kinetic_temperature(
  readings : Array[Reading],
  activation_energy_j_per_mol? : Double = 83144.0,
) -> Double? {
  if readings.length() < 2 {
    return None
  }
  let gas_constant = 8.314472
  let mut weighted_exponential = 0.0
  let mut total_seconds = 0L
  for index = 1; index < readings.length(); index = index + 1 {
    let previous = readings[index - 1]
    let current = readings[index]
    if !previous.is_usable() {
      continue
    }
    let duration = current.timestamp - previous.timestamp
    let kelvin = previous.temperature_c + 273.15
    if duration <= 0L || kelvin <= 0.0 {
      continue
    }
    weighted_exponential = weighted_exponential +
      @math.exp(-activation_energy_j_per_mol / (gas_constant * kelvin)) *
      duration.to_double()
    total_seconds = total_seconds + duration
  }
  if total_seconds == 0L || weighted_exponential <= 0.0 {
    return None
  }
  let average_exponential = weighted_exponential / total_seconds.to_double()
  let kelvin = -activation_energy_j_per_mol /
    (gas_constant * @math.ln(average_exponential))
  Some(kelvin - 273.15)
}

///|
/// Integrate time and degree-seconds by temperature band.
fn exposure_buckets(
  readings : Array[Reading],
  policy : TemperaturePolicy,
) -> (Int64, Int64, Int64, Int64, Double, Double) {
  let mut in_range_seconds = 0L
  let mut low_seconds = 0L
  let mut high_seconds = 0L
  let mut unknown_seconds = 0L
  let mut low_degree_seconds = 0.0
  let mut high_degree_seconds = 0.0
  for index = 1; index < readings.length(); index = index + 1 {
    let previous = readings[index - 1]
    let current = readings[index]
    let duration = current.timestamp - previous.timestamp
    if duration <= 0L {
      continue
    }
    if !previous.is_usable() || current.has_flag(GapBefore) {
      unknown_seconds = unknown_seconds + duration
      continue
    }
    if previous.temperature_c < policy.lower_c {
      low_seconds = low_seconds + duration
      low_degree_seconds = low_degree_seconds +
        (policy.lower_c - previous.temperature_c) * duration.to_double()
    } else if previous.temperature_c > policy.upper_c {
      high_seconds = high_seconds + duration
      high_degree_seconds = high_degree_seconds +
        (previous.temperature_c - policy.upper_c) * duration.to_double()
    } else {
      in_range_seconds = in_range_seconds + duration
    }
  }
  (
    in_range_seconds, low_seconds, high_seconds, unknown_seconds, low_degree_seconds,
    high_degree_seconds,
  )
}

///|
/// Compute all summary statistics for one sensor.
pub fn compute_sensor_statistics(
  sensor_id : String,
  input : Array[Reading],
  policy : TemperaturePolicy,
) -> SensorStatistics {
  let readings = sort_readings(readings_for_sensor(input, sensor_id))
  let values = usable_temperatures(readings)
  let (minimum, maximum) = extrema(values)
  let (weighted_mean, _) = time_weighted_temperature(readings)
  let (
    in_range_seconds,
    low_seconds,
    high_seconds,
    unknown_seconds,
    low_degree_seconds,
    high_degree_seconds,
  ) = exposure_buckets(readings, policy)
  {
    sensor_id,
    sample_count: readings.length(),
    usable_count: values.length(),
    first_timestamp: if readings.length() == 0 {
      None
    } else {
      Some(readings[0].timestamp)
    },
    last_timestamp: if readings.length() == 0 {
      None
    } else {
      Some(readings[readings.length() - 1].timestamp)
    },
    minimum_c: minimum,
    maximum_c: maximum,
    mean_c: mean(values),
    time_weighted_mean_c: weighted_mean,
    mean_kinetic_temperature_c: mean_kinetic_temperature(readings),
    standard_deviation_c: standard_deviation(values),
    median_c: median(values),
    p05_c: percentile(values, 5.0),
    p95_c: percentile(values, 95.0),
    in_range_seconds,
    low_seconds,
    high_seconds,
    unknown_seconds,
    low_degree_seconds,
    high_degree_seconds,
  }
}

///|
/// Compute summaries for every sensor.
pub fn compute_all_statistics(
  readings : Array[Reading],
  policy : TemperaturePolicy,
) -> Array[SensorStatistics] {
  let output : Array[SensorStatistics] = []
  for sensor_id in unique_sensor_ids(readings) {
    output.push(compute_sensor_statistics(sensor_id, readings, policy))
  }
  output
}

///|
/// Fixed-width time windows anchored to the first sample.
pub fn compute_windows(
  sensor_id : String,
  input : Array[Reading],
  policy : TemperaturePolicy,
  width_seconds : Int64,
) -> Array[WindowStatistic] {
  let windows : Array[WindowStatistic] = []
  if width_seconds <= 0L {
    return windows
  }
  let readings = sort_readings(readings_for_sensor(input, sensor_id))
  if readings.length() == 0 {
    return windows
  }
  let anchor = readings[0].timestamp
  let groups : Map[Int64, Array[Reading]] = Map([])
  for sample in readings {
    let index = (sample.timestamp - anchor) / width_seconds
    match groups.get(index) {
      Some(group) => group.push(sample)
      None => groups[index] = [sample]
    }
  }
  let indices = groups.keys().to_array()
  indices.sort()
  for index in indices {
    let group = groups[index]
    let values = usable_temperatures(group)
    let (minimum, maximum) = extrema(values)
    let (_, _, _, _, low_degree_seconds, high_degree_seconds) = exposure_buckets(
      group, policy,
    )
    windows.push({
      sensor_id,
      started_at: anchor + index * width_seconds,
      ended_at: anchor + (index + 1L) * width_seconds,
      sample_count: group.length(),
      minimum_c: minimum,
      maximum_c: maximum,
      mean_c: mean(values),
      low_degree_seconds,
      high_degree_seconds,
    })
  }
  windows
}

///|
/// Convert a numeric score to a non-regulatory risk band.
pub fn risk_band_from_score(score : Int, confidence_percent : Int) -> RiskBand {
  if confidence_percent < 25 {
    return Indeterminate
  }
  if score < 10 {
    Minimal
  } else if score < 25 {
    Low
  } else if score < 50 {
    Moderate
  } else if score < 75 {
    High
  } else {
    Critical
  }
}

///|
/// Build an explainable risk assessment from events, gaps and data diagnostics.
pub fn assess_risk(
  events : Array[ExcursionEvent],
  gaps : Array[SamplingGap],
  diagnostics : Array[Diagnostic],
  health : Array[SensorHealth],
) -> RiskAssessment {
  let confirmed = confirmed_events(events)
  let mut temperature_component = 0
  let mut duration_component = 0
  let mut data_quality_component = 0
  let mut sensor_health_component = 0
  let reasons : Array[String] = []
  let peak = maximum_peak_deviation(confirmed)
  temperature_component = clamp_int((peak * 8.0).round().to_int(), 0, 30)
  if temperature_component > 0 {
    reasons.push(
      "peak temperature deviation contributed \{temperature_component} points",
    )
  }
  let duration = total_event_duration(confirmed, LowTemperature) +
    total_event_duration(confirmed, HighTemperature)
  duration_component = clamp_int((duration / 900L).to_int() * 3, 0, 30)
  if duration_component > 0 {
    reasons.push(
      "confirmed excursion duration contributed \{duration_component} points",
    )
  }
  data_quality_component = clamp_int(
    gaps.length() * 3 +
    count_diagnostics(diagnostics, Error) * 5 +
    count_diagnostics(diagnostics, Warning),
    0,
    20,
  )
  if data_quality_component > 0 {
    reasons.push(
      "data completeness and diagnostics contributed \{data_quality_component} points",
    )
  }
  if health.length() > 0 {
    let mut deficit_total = 0
    for item in health {
      deficit_total = deficit_total + (100 - clamp_int(item.score, 0, 100))
    }
    sensor_health_component = clamp_int(
      deficit_total / health.length() / 5,
      0,
      20,
    )
    if sensor_health_component > 0 {
      reasons.push(
        "sensor health contributed \{sensor_health_component} points",
      )
    }
  }
  let score = clamp_int(
    temperature_component +
    duration_component +
    data_quality_component +
    sensor_health_component,
    0,
    100,
  )
  let confidence = clamp_int(
    100 -
    gaps.length() * 8 -
    count_diagnostics(diagnostics, Error) * 15 -
    count_diagnostics(diagnostics, Warning) * 2,
    0,
    100,
  )
  if reasons.length() == 0 {
    reasons.push("no material risk factors were detected")
  }
  {
    score,
    band: risk_band_from_score(score, confidence),
    temperature_component,
    duration_component,
    data_quality_component,
    sensor_health_component,
    confidence_percent: confidence,
    reasons,
  }
}