///|
/// Resampling diagnostics for estimating estimator variability and bias.
pub struct ResamplingDiagnostic {
  estimate : Double
  resample_mean : Double
  bias : Double
  standard_error : Double
  lower : Double
  upper : Double
  replicates : Int
}

///|
pub fn resampling_diagnostic(
  estimate : Double,
  samples : Array[Double],
  confidence : Double,
) -> ResamplingDiagnostic {
  let alpha = (1.0 - transform_clip(confidence, 0.0, 1.0)) / 2.0
  {
    estimate,
    resample_mean: mean(samples),
    bias: mean(samples) - estimate,
    standard_error: sample_stddev(samples),
    lower: quantile(samples, alpha),
    upper: quantile(samples, 1.0 - alpha),
    replicates: samples.length(),
  }
}

///|
pub fn resampling_mean_diagnostic(
  data : Array[Double],
  replicates : Int,
  seed : Int,
  confidence : Double,
) -> ResamplingDiagnostic {
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(mean(sample))
  }
  resampling_diagnostic(mean(data), estimates, confidence)
}

///|
pub fn resampling_median_diagnostic(
  data : Array[Double],
  replicates : Int,
  seed : Int,
  confidence : Double,
) -> ResamplingDiagnostic {
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(median(sample))
  }
  resampling_diagnostic(median(data), estimates, confidence)
}

///|
pub fn resampling_diagnostic_vector(
  diagnostic : ResamplingDiagnostic,
) -> Array[Double] {
  [
    diagnostic.estimate,
    diagnostic.resample_mean,
    diagnostic.bias,
    diagnostic.standard_error,
    diagnostic.lower,
    diagnostic.upper,
    diagnostic.replicates.to_double(),
  ]
}

///|
pub fn resampling_diagnostic_lines(
  diagnostic : ResamplingDiagnostic,
) -> Array[String] {
  [
    "estimate=" + diagnostic.estimate.to_string(),
    "resample_mean=" + diagnostic.resample_mean.to_string(),
    "bias=" + diagnostic.bias.to_string(),
    "standard_error=" + diagnostic.standard_error.to_string(),
    "lower=" + diagnostic.lower.to_string(),
    "upper=" + diagnostic.upper.to_string(),
    "replicates=" + diagnostic.replicates.to_string(),
  ]
}

///|
pub fn resampling_diagnostic_string(
  diagnostic : ResamplingDiagnostic,
) -> String {
  resampling_diagnostic_lines(diagnostic).join("\n")
}

///|
pub fn resampling_interval_width(diagnostic : ResamplingDiagnostic) -> Double {
  diagnostic.upper - diagnostic.lower
}

///|
pub fn resampling_coverage(
  diagnostic : ResamplingDiagnostic,
  observations : Array[Double],
) -> Double {
  coverage_of_interval(observations, [diagnostic.lower, diagnostic.upper])
}

///|
pub fn resampling_relative_bias(diagnostic : ResamplingDiagnostic) -> Double {
  if diagnostic.estimate == 0.0 {
    diagnostic.bias
  } else {
    diagnostic.bias / diagnostic.estimate
  }
}

///|
pub fn resampling_stable(
  left : ResamplingDiagnostic,
  right : ResamplingDiagnostic,
  tolerance : Double,
) -> Bool {
  abs_double(left.estimate - right.estimate) <= tolerance &&
  abs_double(left.standard_error - right.standard_error) <= tolerance
}

///|
pub fn resampling_compare(
  left : ResamplingDiagnostic,
  right : ResamplingDiagnostic,
) -> Array[Double] {
  [
    left.standard_error,
    right.standard_error,
    resampling_interval_width(left),
    resampling_interval_width(right),
    left.bias,
    right.bias,
  ]
}

///|
pub fn resampling_jackknife_diagnostic(
  data : Array[Double],
  confidence : Double,
) -> ResamplingDiagnostic {
  let samples = sampling_jackknife(data)
  let estimates = []
  for sample in samples {
    estimates.push(mean(sample))
  }
  resampling_diagnostic(mean(data), estimates, confidence)
}

///|
pub fn resampling_permutation_diagnostic(
  left : Array[Double],
  right : Array[Double],
  replicates : Int,
  seed : Int,
) -> ResamplingDiagnostic {
  let observed = mean(left) - mean(right)
  let samples = sampling_permutation_difference(left, right, replicates, seed)
  resampling_diagnostic(observed, samples, 0.95)
}

///|
pub fn resampling_batch(
  data_sets : Array[Array[Double]],
  replicates : Int,
  seed : Int,
) -> Array[Double] {
  let result = []
  for index = 0; index < data_sets.length(); index = index + 1 {
    result.push(
      resampling_mean_diagnostic(
        data_sets[index],
        replicates,
        seed + index,
        0.95,
      ).standard_error,
    )
  }
  result
}

///|
pub fn resampling_quality(
  diagnostic : ResamplingDiagnostic,
  tolerance : Double,
) -> Double {
  let stability = if diagnostic.standard_error <= tolerance {
    1.0
  } else {
    tolerance / diagnostic.standard_error
  }
  transform_clip(stability, 0.0, 1.0)
}