///|
/// A calibration bin for probabilistic binary predictions.
pub struct CalibrationBin {
  lower : Double
  upper : Double
  count : Int
  positives : Int
  mean_probability : Double
  observed_rate : Double
  gap : Double
}

///|
/// Classification metrics for one threshold.
pub struct ClassificationReport {
  threshold : Double
  matrix : ConfusionMatrix
  accuracy : Double
  balanced_accuracy : Double
  precision : Double
  recall : Double
  specificity : Double
  f1 : Double
  matthews : Double
  brier : Double
  calibration_error : Double
}

///|
/// Cost-sensitive threshold rule.
pub struct ThresholdRule {
  false_positive_cost : Double
  false_negative_cost : Double
  minimum_recall : Double
  minimum_precision : Double
}

///|
pub fn classification_threshold_rule(
  false_positive_cost : Double,
  false_negative_cost : Double,
  minimum_recall : Double,
  minimum_precision : Double,
) -> ThresholdRule {
  {
    false_positive_cost: if false_positive_cost < 0.0 {
      0.0
    } else {
      false_positive_cost
    },
    false_negative_cost: if false_negative_cost < 0.0 {
      0.0
    } else {
      false_negative_cost
    },
    minimum_recall: if minimum_recall < 0.0 {
      0.0
    } else if minimum_recall > 1.0 {
      1.0
    } else {
      minimum_recall
    },
    minimum_precision: if minimum_precision < 0.0 {
      0.0
    } else if minimum_precision > 1.0 {
      1.0
    } else {
      minimum_precision
    },
  }
}

///|
pub fn classification_probability(value : Double) -> Double {
  if value < 0.0 {
    0.0
  } else if value > 1.0 {
    1.0
  } else {
    value
  }
}

///|
pub fn classification_labels(
  probabilities : Array[Double],
  threshold : Double,
) -> Array[Bool] {
  let result = []
  let cutoff = classification_probability(threshold)
  for probability in probabilities {
    result.push(classification_probability(probability) >= cutoff)
  }
  result
}

///|
pub fn classification_positive_count(labels : Array[Bool]) -> Int {
  let mut count = 0
  for label in labels {
    if label {
      count += 1
    }
  }
  count
}

///|
pub fn classification_prevalence(labels : Array[Bool]) -> Double {
  if labels.length() == 0 {
    0.0
  } else {
    classification_positive_count(labels).to_double() /
    labels.length().to_double()
  }
}

///|
pub fn classification_threshold_report(
  actual : Array[Bool],
  probabilities : Array[Double],
  threshold : Double,
  bins : Int,
) -> ClassificationReport {
  let predicted = classification_labels(probabilities, threshold)
  let matrix = confusion_matrix(actual, predicted)
  let calibration = classification_calibration_bins(actual, probabilities, bins)
  {
    threshold: classification_probability(threshold),
    matrix,
    accuracy: accuracy(matrix),
    balanced_accuracy: balanced_accuracy(matrix),
    precision: precision(matrix),
    recall: recall(matrix),
    specificity: specificity(matrix),
    f1: f1_score(matrix),
    matthews: matthews_correlation(matrix),
    brier: brier_score(actual, probabilities),
    calibration_error: classification_expected_calibration_error(calibration),
  }
}

///|
pub fn classification_calibration_bins(
  actual : Array[Bool],
  probabilities : Array[Double],
  bins : Int,
) -> Array[CalibrationBin] {
  let count = if bins < 1 { 1 } else { bins }
  let result = []
  for bin = 0; bin < count; bin = bin + 1 {
    let lower = bin.to_double() / count.to_double()
    let upper = (bin + 1).to_double() / count.to_double()
    let mut observations = 0
    let mut positives = 0
    let mut probability_total = 0.0
    let limit = if actual.length() < probabilities.length() {
      actual.length()
    } else {
      probabilities.length()
    }
    for index = 0; index < limit; index = index + 1 {
      let probability = classification_probability(probabilities[index])
      let in_bin = if bin == count - 1 {
        probability >= lower && probability <= upper
      } else {
        probability >= lower && probability < upper
      }
      if in_bin {
        observations += 1
        if actual[index] {
          positives += 1
        }
        probability_total += probability
      }
    }
    let mean_probability = if observations == 0 {
      0.0
    } else {
      probability_total / observations.to_double()
    }
    let observed_rate = if observations == 0 {
      0.0
    } else {
      positives.to_double() / observations.to_double()
    }
    result.push({
      lower,
      upper,
      count: observations,
      positives,
      mean_probability,
      observed_rate,
      gap: abs_double(mean_probability - observed_rate),
    })
  }
  result
}

///|
pub fn classification_expected_calibration_error(
  bins : Array[CalibrationBin],
) -> Double {
  let mut total = 0
  for bin in bins {
    total += bin.count
  }
  if total == 0 {
    return 0.0
  }
  let mut error = 0.0
  for bin in bins {
    error += bin.count.to_double() * bin.gap
  }
  error / total.to_double()
}

///|
pub fn classification_max_calibration_error(
  bins : Array[CalibrationBin],
) -> Double {
  let mut result = 0.0
  for bin in bins {
    if bin.gap > result {
      result = bin.gap
    }
  }
  result
}

///|
pub fn classification_log_loss(
  actual : Array[Bool],
  probabilities : Array[Double],
) -> Double {
  let limit = if actual.length() < probabilities.length() {
    actual.length()
  } else {
    probabilities.length()
  }
  if limit == 0 {
    return 0.0
  }
  let mut total = 0.0
  for index = 0; index < limit; index = index + 1 {
    let probability = classification_probability(probabilities[index])
    let clipped = if probability < 1.0e-12 {
      1.0e-12
    } else if probability > 1.0 - 1.0e-12 {
      1.0 - 1.0e-12
    } else {
      probability
    }
    if actual[index] {
      total -= drift_log(clipped)
    } else {
      total -= drift_log(1.0 - clipped)
    }
  }
  total / limit.to_double()
}

///|
pub fn classification_fbeta(matrix : ConfusionMatrix, beta : Double) -> Double {
  let weight = if beta < 0.0 { 0.0 } else { beta * beta }
  let p = precision(matrix)
  let r = recall(matrix)
  let denominator = weight * p + r
  if denominator == 0.0 {
    0.0
  } else {
    (1.0 + weight) * p * r / denominator
  }
}

///|
pub fn classification_cost(
  matrix : ConfusionMatrix,
  rule : ThresholdRule,
) -> Double {
  matrix.false_positive.to_double() * rule.false_positive_cost +
  matrix.false_negative.to_double() * rule.false_negative_cost
}

///|
pub fn classification_reports(
  actual : Array[Bool],
  probabilities : Array[Double],
  thresholds : Array[Double],
) -> Array[ClassificationReport] {
  let result = []
  for threshold in thresholds {
    result.push(
      classification_threshold_report(actual, probabilities, threshold, 10),
    )
  }
  result
}

///|
pub fn classification_best_f1_threshold(
  actual : Array[Bool],
  probabilities : Array[Double],
  thresholds : Array[Double],
) -> Double {
  if thresholds.length() == 0 {
    return 0.5
  }
  let mut best = classification_probability(thresholds[0])
  let mut score = -1.0
  for threshold in thresholds {
    let matrix = confusion_matrix(
      actual,
      classification_labels(probabilities, threshold),
    )
    let current = f1_score(matrix)
    if current > score {
      score = current
      best = classification_probability(threshold)
    }
  }
  best
}

///|
pub fn classification_best_cost_threshold(
  actual : Array[Bool],
  probabilities : Array[Double],
  thresholds : Array[Double],
  rule : ThresholdRule,
) -> Double {
  if thresholds.length() == 0 {
    return 0.5
  }
  let mut best = classification_probability(thresholds[0])
  let mut score = 1.0e300
  for threshold in thresholds {
    let matrix = confusion_matrix(
      actual,
      classification_labels(probabilities, threshold),
    )
    let current_recall = recall(matrix)
    let current_precision = precision(matrix)
    let current_cost = classification_cost(matrix, rule)
    if current_recall >= rule.minimum_recall &&
      current_precision >= rule.minimum_precision &&
      current_cost < score {
      score = current_cost
      best = classification_probability(threshold)
    }
  }
  best
}

///|
pub fn classification_threshold_grid(steps : Int) -> Array[Double] {
  let count = if steps < 1 { 1 } else { steps }
  let result = []
  for index = 0; index <= count; index = index + 1 {
    result.push(index.to_double() / count.to_double())
  }
  result
}

///|
pub fn classification_top_k_labels(
  probabilities : Array[Double],
  k : Int,
) -> Array[Bool] {
  let result = []
  for _ in probabilities {
    result.push(false)
  }
  let ranks = rank_of_values(probabilities)
  let limit = if k < 0 {
    0
  } else if k > probabilities.length() {
    probabilities.length()
  } else {
    k
  }
  for index = 0; index < limit && index < ranks.length(); index = index + 1 {
    let position = ranks[index]
    if position >= 0 && position < result.length() {
      result[position] = true
    }
  }
  result
}

///|
pub fn rank_of_values(data : Array[Double]) -> Array[Int] {
  let order = []
  for index = 0; index < data.length(); index = index + 1 {
    order.push(index)
  }
  for left = 0; left < order.length(); left = left + 1 {
    for right = left + 1; right < order.length(); right = right + 1 {
      if data[order[right]] > data[order[left]] {
        let value = order[left]
        order[left] = order[right]
        order[right] = value
      }
    }
  }
  order
}

///|
pub fn classification_precision_at_k(
  actual : Array[Bool],
  probabilities : Array[Double],
  k : Int,
) -> Double {
  let predicted = classification_top_k_labels(probabilities, k)
  precision(confusion_matrix(actual, predicted))
}

///|
pub fn classification_recall_at_k(
  actual : Array[Bool],
  probabilities : Array[Double],
  k : Int,
) -> Double {
  let predicted = classification_top_k_labels(probabilities, k)
  recall(confusion_matrix(actual, predicted))
}

///|
pub fn classification_lift_at_k(
  actual : Array[Bool],
  probabilities : Array[Double],
  k : Int,
) -> Double {
  let baseline = classification_prevalence(actual)
  let observed = classification_precision_at_k(actual, probabilities, k)
  if baseline == 0.0 {
    0.0
  } else {
    observed / baseline
  }
}

///|
pub fn classification_gain_curve(
  actual : Array[Bool],
  probabilities : Array[Double],
  steps : Int,
) -> Array[Array[Double]] {
  let result = []
  let count = if steps < 1 { 1 } else { steps }
  for index = 1; index <= count; index = index + 1 {
    let k = probabilities.length() * index / count
    let predicted = classification_top_k_labels(probabilities, k)
    let matrix = confusion_matrix(actual, predicted)
    result.push([
      k.to_double() / probabilities.length().to_double(),
      recall(matrix),
      precision(matrix),
      classification_lift_at_k(actual, probabilities, k),
    ])
  }
  result
}

///|
pub fn classification_positive_rate_by_threshold(
  probabilities : Array[Double],
  thresholds : Array[Double],
) -> Array[Double] {
  let result = []
  for threshold in thresholds {
    result.push(
      classification_prevalence(classification_labels(probabilities, threshold)),
    )
  }
  result
}

///|
pub fn classification_monotone_probabilities(
  probabilities : Array[Double],
) -> Bool {
  for index = 1; index < probabilities.length(); index = index + 1 {
    if classification_probability(probabilities[index]) <
      classification_probability(probabilities[index - 1]) {
      return false
    }
  }
  true
}

///|
pub fn classification_reliability_curve(
  actual : Array[Bool],
  probabilities : Array[Double],
  bins : Int,
) -> Array[Array[Double]] {
  let result = []
  for bin in classification_calibration_bins(actual, probabilities, bins) {
    result.push([
      bin.mean_probability,
      bin.observed_rate,
      bin.count.to_double(),
      bin.gap,
    ])
  }
  result
}

///|
pub fn classification_report_vector(
  report : ClassificationReport,
) -> Array[Double] {
  [
    report.threshold,
    report.accuracy,
    report.balanced_accuracy,
    report.precision,
    report.recall,
    report.specificity,
    report.f1,
    report.matthews,
    report.brier,
    report.calibration_error,
  ]
}

///|
pub fn classification_report_lines(
  report : ClassificationReport,
) -> Array[String] {
  [
    "threshold=" + report.threshold.to_string(),
    "accuracy=" + report.accuracy.to_string(),
    "balanced_accuracy=" + report.balanced_accuracy.to_string(),
    "precision=" + report.precision.to_string(),
    "recall=" + report.recall.to_string(),
    "specificity=" + report.specificity.to_string(),
    "f1=" + report.f1.to_string(),
    "matthews=" + report.matthews.to_string(),
    "brier=" + report.brier.to_string(),
    "calibration_error=" + report.calibration_error.to_string(),
  ]
}

///|
pub fn classification_report_string(report : ClassificationReport) -> String {
  classification_report_lines(report).join("\n")
}

///|
pub fn classification_confidence_margin(
  probabilities : Array[Double],
  threshold : Double,
) -> Array[Double] {
  let result = []
  for probability in probabilities {
    result.push(abs_double(classification_probability(probability) - threshold))
  }
  result
}

///|
pub fn classification_uncertain_indices(
  probabilities : Array[Double],
  threshold : Double,
  margin : Double,
) -> Array[Int] {
  let result = []
  let limit = if margin < 0.0 { -margin } else { margin }
  for index = 0; index < probabilities.length(); index = index + 1 {
    if abs_double(classification_probability(probabilities[index]) - threshold) <=
      limit {
      result.push(index)
    }
  }
  result
}

///|
pub fn classification_disagreement(
  first : Array[Double],
  second : Array[Double],
  threshold : Double,
) -> Array[Int] {
  let result = []
  let limit = if first.length() < second.length() {
    first.length()
  } else {
    second.length()
  }
  for index = 0; index < limit; index = index + 1 {
    let left = classification_probability(first[index]) >= threshold
    let right = classification_probability(second[index]) >= threshold
    if left != right {
      result.push(index)
    }
  }
  result
}

///|
pub fn classification_report_quality(report : ClassificationReport) -> Double {
  (
    report.f1 +
    report.balanced_accuracy +
    (1.0 - report.brier) +
    (1.0 - report.calibration_error)
  ) /
  4.0
}