///|
/// Distribution-free interval estimate used for robust uncertainty reporting.
pub struct ProbabilityInterval {
  estimate : Double
  lower : Double
  upper : Double
  confidence : Double
  samples : Int
}

///|
pub struct ProbabilityMass {
  values : Array[Double]
  probabilities : Array[Double]
  entropy : Double
  effective_count : Double
}

///|
pub struct ProbabilityCalibration {
  nominal : Double
  empirical : Double
  error : Double
  sample_count : Int
}

///|
pub fn probability_make_interval(
  estimate : Double,
  lower : Double,
  upper : Double,
  confidence : Double,
  samples : Int,
) -> ProbabilityInterval {
  {
    estimate,
    lower: if lower > upper {
      upper
    } else {
      lower
    },
    upper: if upper < lower {
      lower
    } else {
      upper
    },
    confidence: transform_clip(confidence, 0.0, 1.0),
    samples: if samples < 0 {
      0
    } else {
      samples
    },
  }
}

///|
pub fn probability_interval_vector(
  interval : ProbabilityInterval,
) -> Array[Double] {
  [
    interval.estimate,
    interval.lower,
    interval.upper,
    interval.confidence,
    interval.samples.to_double(),
    interval.upper - interval.lower,
  ]
}

///|
pub fn probability_interval_lines(
  interval : ProbabilityInterval,
) -> Array[String] {
  [
    "estimate=" + interval.estimate.to_string(),
    "lower=" + interval.lower.to_string(),
    "upper=" + interval.upper.to_string(),
    "confidence=" + interval.confidence.to_string(),
    "samples=" + interval.samples.to_string(),
    "width=" + (interval.upper - interval.lower).to_string(),
  ]
}

///|
pub fn probability_interval_string(interval : ProbabilityInterval) -> String {
  probability_interval_lines(interval).join("\n")
}

///|
pub fn probability_bootstrap_mean(
  data : Array[Double],
  replicates : Int,
  seed : Int,
  confidence : Double,
) -> ProbabilityInterval {
  let samples = bootstrap_replicates(
    data,
    if replicates < 0 {
      0
    } else {
      replicates
    },
    sample_size=data.length(),
    seed~,
  )
  let estimates = []
  for sample in samples {
    estimates.push(mean(sample))
  }
  if estimates.length() == 0 {
    return probability_make_interval(
      mean(data),
      mean(data),
      mean(data),
      confidence,
      0,
    )
  }
  let alpha = (1.0 - transform_clip(confidence, 0.0, 1.0)) / 2.0
  probability_make_interval(
    mean(data),
    quantile(estimates, alpha),
    quantile(estimates, 1.0 - alpha),
    confidence,
    estimates.length(),
  )
}

///|
pub fn probability_bootstrap_median(
  data : Array[Double],
  replicates : Int,
  seed : Int,
  confidence : Double,
) -> ProbabilityInterval {
  let samples = bootstrap_replicates(
    data,
    if replicates < 0 {
      0
    } else {
      replicates
    },
    sample_size=data.length(),
    seed~,
  )
  let estimates = []
  for sample in samples {
    estimates.push(median(sample))
  }
  if estimates.length() == 0 {
    return probability_make_interval(
      median(data),
      median(data),
      median(data),
      confidence,
      0,
    )
  }
  let alpha = (1.0 - transform_clip(confidence, 0.0, 1.0)) / 2.0
  probability_make_interval(
    median(data),
    quantile(estimates, alpha),
    quantile(estimates, 1.0 - alpha),
    confidence,
    estimates.length(),
  )
}

///|
pub fn probability_jackknife_mean(
  data : Array[Double],
  confidence : Double,
) -> ProbabilityInterval {
  let samples = sampling_jackknife(data)
  let estimates = []
  for sample in samples {
    estimates.push(mean(sample))
  }
  let center = mean(data)
  if estimates.length() < 2 {
    return probability_make_interval(
      center,
      center,
      center,
      confidence,
      estimates.length(),
    )
  }
  let scale = sample_stddev(estimates)
  let width = 1.96 * scale
  probability_make_interval(
    center,
    center - width,
    center + width,
    confidence,
    estimates.length(),
  )
}

///|
pub fn probability_empirical(data : Array[Double], value : Double) -> Double {
  empirical_cdf(data, value)
}

///|
pub fn probability_tail(data : Array[Double], value : Double) -> Double {
  empirical_survival(data, value)
}

///|
pub fn probability_two_sided(data : Array[Double], value : Double) -> Double {
  let left = probability_empirical(data, value)
  let right = probability_tail(data, value)
  if left < right {
    2.0 * left
  } else {
    2.0 * right
  }
}

///|
pub fn probability_mass(
  values : Array[Double],
  weights : Array[Double],
) -> ProbabilityMass {
  let probabilities = []
  let mut total = 0.0
  let count = if values.length() < weights.length() {
    values.length()
  } else {
    weights.length()
  }
  for index = 0; index < count; index = index + 1 {
    total += if weights[index] < 0.0 { 0.0 } else { weights[index] }
  }
  for index = 0; index < count; index = index + 1 {
    let weight = if weights[index] < 0.0 { 0.0 } else { weights[index] }
    probabilities.push(if total == 0.0 { 0.0 } else { weight / total })
  }
  let mut entropy = 0.0
  let mut square_total = 0.0
  for probability in probabilities {
    if probability > 0.0 {
      entropy -= probability * drift_log(probability)
    }
    square_total += probability * probability
  }
  {
    values: values.copy(),
    probabilities,
    entropy,
    effective_count: if square_total == 0.0 {
      0.0
    } else {
      1.0 / square_total
    },
  }
}

///|
pub fn probability_mass_normalized(mass : ProbabilityMass) -> Bool {
  let mut total = 0.0
  for probability in mass.probabilities {
    total += probability
  }
  abs_double(total - 1.0) < 1.0e-9 || mass.probabilities.length() == 0
}

///|
pub fn probability_mass_mean(mass : ProbabilityMass) -> Double {
  let mut total = 0.0
  let count = if mass.values.length() < mass.probabilities.length() {
    mass.values.length()
  } else {
    mass.probabilities.length()
  }
  for index = 0; index < count; index = index + 1 {
    total += mass.values[index] * mass.probabilities[index]
  }
  total
}

///|
pub fn probability_mass_variance(mass : ProbabilityMass) -> Double {
  let center = probability_mass_mean(mass)
  let mut total = 0.0
  let count = if mass.values.length() < mass.probabilities.length() {
    mass.values.length()
  } else {
    mass.probabilities.length()
  }
  for index = 0; index < count; index = index + 1 {
    let delta = mass.values[index] - center
    total += delta * delta * mass.probabilities[index]
  }
  total
}

///|
pub fn probability_mass_quantile(
  mass : ProbabilityMass,
  probability : Double,
) -> Double {
  if mass.values.length() == 0 {
    return 0.0
  }
  let pairs = []
  let count = if mass.values.length() < mass.probabilities.length() {
    mass.values.length()
  } else {
    mass.probabilities.length()
  }
  for index = 0; index < count; index = index + 1 {
    pairs.push([mass.values[index], mass.probabilities[index]])
  }
  pairs.sort_by((left, right) => {
    if left[0] < right[0] {
      -1
    } else if left[0] > right[0] {
      1
    } else {
      0
    }
  })
  let target = transform_clip(probability, 0.0, 1.0)
  let mut cumulative = 0.0
  for pair in pairs {
    cumulative += pair[1]
    if cumulative >= target {
      return pair[0]
    }
  }
  mass.values[mass.values.length() - 1]
}

///|
pub fn probability_calibration(
  nominal : Array[Double],
  observed : Array[Bool],
  bins : Int,
) -> Array[ProbabilityCalibration] {
  let result = []
  let count = if nominal.length() < observed.length() {
    nominal.length()
  } else {
    observed.length()
  }
  let width = if bins < 1 { 1 } else { bins }
  for bin = 0; bin < width; bin = bin + 1 {
    let lower = bin.to_double() / width.to_double()
    let upper = (bin + 1).to_double() / width.to_double()
    let mut sum_probability = 0.0
    let mut positives = 0
    let mut samples = 0
    for index = 0; index < count; index = index + 1 {
      if nominal[index] >= lower && (nominal[index] < upper || bin == width - 1) {
        sum_probability += nominal[index]
        if observed[index] {
          positives += 1
        }
        samples += 1
      }
    }
    let empirical = if samples == 0 {
      0.0
    } else {
      positives.to_double() / samples.to_double()
    }
    let average = if samples == 0 {
      0.0
    } else {
      sum_probability / samples.to_double()
    }
    result.push({
      nominal: average,
      empirical,
      error: abs_double(average - empirical),
      sample_count: samples,
    })
  }
  result
}

///|
pub fn probability_expected_calibration_error(
  calibration : Array[ProbabilityCalibration],
) -> Double {
  let mut total = 0.0
  let mut count = 0
  for item in calibration {
    total += item.error * item.sample_count.to_double()
    count += item.sample_count
  }
  if count == 0 {
    0.0
  } else {
    total / count.to_double()
  }
}

///|
pub fn probability_max_calibration_error(
  calibration : Array[ProbabilityCalibration],
) -> Double {
  let mut result = 0.0
  for item in calibration {
    if item.error > result {
      result = item.error
    }
  }
  result
}

///|
pub fn probability_interval_coverage(
  interval : ProbabilityInterval,
  observations : Array[Double],
) -> Double {
  coverage_of_interval(observations, [interval.lower, interval.upper])
}

///|
pub fn probability_interval_score(
  interval : ProbabilityInterval,
  observations : Array[Double],
) -> Double {
  let coverage = probability_interval_coverage(interval, observations)
  let width = interval.upper - interval.lower
  coverage - width * (1.0 - interval.confidence)
}

///|
pub fn probability_reliability_score(
  calibration : Array[ProbabilityCalibration],
) -> Double {
  1.0 -
  transform_clip(probability_expected_calibration_error(calibration), 0.0, 1.0)
}

///|
pub fn probability_robust_interval(
  data : Array[Double],
  confidence : Double,
) -> ProbabilityInterval {
  let alpha = (1.0 - transform_clip(confidence, 0.0, 1.0)) / 2.0
  probability_make_interval(
    median(data),
    quantile(data, alpha),
    quantile(data, 1.0 - alpha),
    confidence,
    data.length(),
  )
}

///|
pub fn probability_summary(
  data : Array[Double],
  confidence : Double,
) -> Array[Double] {
  let interval = probability_robust_interval(data, confidence)
  [
    interval.estimate,
    interval.lower,
    interval.upper,
    interval.upper - interval.lower,
    probability_empirical(data, interval.estimate),
    probability_tail(data, interval.upper),
  ]
}