///|
/// Resampling and uncertainty summaries for calibration models.
pub(all) struct BootstrapFit {
  replicate : Int
  model : RegressionModel
  prediction_at_reference : Double
  residual_scale : Double
  accepted : Bool
} derive(Debug, Eq, ToJson)

///|
pub(all) struct FitUncertainty {
  reference_x : Double
  mean_prediction : Double
  standard_deviation : Double
  lower_prediction : Double
  upper_prediction : Double
  relative_uncertainty : Double
  replicates : Int
} derive(Debug, Eq, ToJson)

///|
pub(all) struct CoefficientSummary {
  index : Int
  mean : Double
  standard_deviation : Double
  minimum : Double
  maximum : Double
} derive(Debug, Eq, ToJson)

///|
pub fn residual_scale(
  model : RegressionModel,
  observations : ArrayView[Observation],
) -> Double {
  if observations.length() == 0 {
    return 0.0
  }
  let mut total = 0.0
  for observation in observations {
    let residual = model.residual(observation)
    total += residual * residual
  }
  (total / observations.length().to_double()).sqrt()
}

///|
pub fn bootstrap_sample(
  observations : ArrayView[Observation],
  replicate : Int,
) -> Array[Observation] raise @domain.ExtractionError {
  if observations.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let sample = []
  let seed = (replicate.abs() + 1) * 17
  for index in 0.. Array[BootstrapFit] raise @domain.ExtractionError {
  if replicates < 1 {
    raise @domain.ExtractionError::InvalidArgument(
      message="bootstrap replicates must be positive",
    )
  }
  let fits = []
  for replicate in 0.. Array[BootstrapFit] {
  let accepted = []
  for fit in fits {
    if fit.accepted {
      accepted.push(fit)
    }
  }
  accepted
}

///|
pub fn summarize_fit_uncertainty(
  fits : ArrayView[BootstrapFit],
  reference_x : Double,
  lower_probability : Double,
  upper_probability : Double,
) -> FitUncertainty raise @domain.ExtractionError {
  if fits.length() == 0 ||
    lower_probability < 0.0 ||
    upper_probability > 1.0 ||
    lower_probability > upper_probability {
    raise @domain.ExtractionError::InvalidArgument(
      message="fit uncertainty inputs are invalid",
    )
  }
  let predictions = [ for fit in fits => fit.prediction_at_reference ]
  predictions.sort()
  let mut mean = 0.0
  for prediction in predictions {
    mean += prediction
  }
  mean /= predictions.length().to_double()
  let mut variance = 0.0
  for prediction in predictions {
    let difference = prediction - mean
    variance += difference * difference
  }
  variance /= predictions.length().to_double()
  let lower_index = (lower_probability * (predictions.length() - 1).to_double()).to_int()
  let upper_index = (upper_probability * (predictions.length() - 1).to_double()).to_int()
  {
    reference_x,
    mean_prediction: mean,
    standard_deviation: variance.sqrt(),
    lower_prediction: predictions[lower_index],
    upper_prediction: predictions[upper_index],
    relative_uncertainty: if mean == 0.0 {
      0.0
    } else {
      variance.sqrt() / mean.abs()
    },
    replicates: predictions.length(),
  }
}

///|
pub fn summarize_coefficients(
  fits : ArrayView[BootstrapFit],
  coefficient_count : Int,
) -> Array[CoefficientSummary] {
  let summaries = []
  for index in 0.. values.push(value)
        None => ()
      }
    }
    if values.length() == 0 {
      summaries.push({
        index,
        mean: 0.0,
        standard_deviation: 0.0,
        minimum: 0.0,
        maximum: 0.0,
      })
    } else {
      let mut total = 0.0
      let mut minimum = values[0]
      let mut maximum = values[0]
      for value in values {
        total += value
        minimum = minimum.min(value)
        maximum = maximum.max(value)
      }
      let mean = total / values.length().to_double()
      let mut variance = 0.0
      for value in values {
        let difference = value - mean
        variance += difference * difference
      }
      summaries.push({
        index,
        mean,
        standard_deviation: (variance / values.length().to_double()).sqrt(),
        minimum,
        maximum,
      })
    }
  }
  summaries
}

///|
pub fn cross_validated_rmse(
  observations : ArrayView[Observation],
  degree : Int,
  folds : Int,
) -> Double raise @domain.ExtractionError {
  if observations.length() < 2 || folds < 2 || folds > observations.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="cross-validation inputs are invalid",
    )
  }
  let mut total_squared_error = 0.0
  let mut count = 0
  for fold in 0.. Array[Double] raise @domain.ExtractionError {
  if minimum_size < degree + 1 || minimum_size > observations.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="learning-curve minimum size is invalid",
    )
  }
  let errors = []
  let mut size = minimum_size
  while size <= observations.length() {
    let subset = []
    for index in 0.. Array[(Double, Double, Double)] {
  let scale = residual_scale(model, observations) * confidence_multiplier
  [
    for x in xs => {
      let prediction = model.predict(x)
      (prediction - scale, prediction, prediction + scale)
    }
  ]
}

///|
pub fn calibration_quality_label(
  rmse : Double,
  relative_uncertainty : Double,
) -> String {
  if rmse <= 1.0e-6 && relative_uncertainty <= 0.01 {
    "excellent"
  } else if rmse <= 1.0e-3 && relative_uncertainty <= 0.05 {
    "good"
  } else if rmse <= 1.0e-2 && relative_uncertainty <= 0.15 {
    "screening"
  } else {
    "needs-data"
  }
}

///|
pub fn fit_uncertainty_markdown(
  uncertainty : FitUncertainty,
  label : String,
) -> String {
  "| Reference x | Mean | Std dev | Lower | Upper | Relative uncertainty | Quality |\n" +
  "|---:|---:|---:|---:|---:|---:|---|\n" +
  "| " +
  uncertainty.reference_x.to_string() +
  " | " +
  uncertainty.mean_prediction.to_string() +
  " | " +
  uncertainty.standard_deviation.to_string() +
  " | " +
  uncertainty.lower_prediction.to_string() +
  " | " +
  uncertainty.upper_prediction.to_string() +
  " | " +
  uncertainty.relative_uncertainty.to_string() +
  " | " +
  label +
  " |"
}

///|
pub fn compare_fit_predictions(
  first : RegressionModel,
  second : RegressionModel,
  xs : ArrayView[Double],
) -> Array[Double] {
  [
    for x in xs => second.predict(x) - first.predict(x)
  ]
}

///|
pub fn maximum_prediction_difference(
  first : RegressionModel,
  second : RegressionModel,
  xs : ArrayView[Double],
) -> Double {
  let differences = compare_fit_predictions(first, second, xs)
  let mut maximum = 0.0
  for difference in differences {
    maximum = maximum.max(difference.abs())
  }
  maximum
}

///|
pub fn bootstrap_acceptance_rate(fits : ArrayView[BootstrapFit]) -> Double {
  if fits.length() == 0 {
    return 0.0
  }
  let mut accepted = 0
  for fit in fits {
    if fit.accepted {
      accepted += 1
    }
  }
  accepted.to_double() / fits.length().to_double()
}