///|
fn safe_ratio(numerator : Int, denominator : Int) -> Double {
  if denominator == 0 {
    0.0
  } else {
    numerator.to_double() / denominator.to_double()
  }
}

///|
fn harmonic_f1(precision : Double, recall : Double) -> Double {
  if precision + recall <= 0.0 {
    0.0
  } else {
    2.0 * precision * recall / (precision + recall)
  }
}

///|
/// Computes one-vs-rest class metrics and macro/weighted aggregates.
pub fn classification_report(
  actual : Array[Int],
  predicted : Array[Int],
  class_count : Int,
) -> Result[ClassificationReport, TreeError] {
  let matrix = match confusion_matrix(actual, predicted, class_count) {
    Ok(value) => value
    Err(error) => return Err(error)
  }
  let total = actual.length()
  let classes : Array[ClassMetrics] = []
  let mut correct = 0
  let mut macro_precision = 0.0
  let mut macro_recall = 0.0
  let mut macro_f1 = 0.0
  let mut weighted_precision = 0.0
  let mut weighted_recall = 0.0
  let mut weighted_f1 = 0.0
  for class_index = 0; class_index < class_count; class_index = class_index + 1 {
    let true_positive = matrix[class_index][class_index]
    correct = correct + true_positive
    let mut support = 0
    let mut predicted_count = 0
    for other = 0; other < class_count; other = other + 1 {
      support = support + matrix[class_index][other]
      predicted_count = predicted_count + matrix[other][class_index]
    }
    let false_positive = predicted_count - true_positive
    let false_negative = support - true_positive
    let true_negative = total - true_positive - false_positive - false_negative
    let precision = safe_ratio(true_positive, true_positive + false_positive)
    let recall = safe_ratio(true_positive, true_positive + false_negative)
    let specificity = safe_ratio(true_negative, true_negative + false_positive)
    let f1 = harmonic_f1(precision, recall)
    classes.push({
      class_index,
      support,
      true_positive,
      false_positive,
      false_negative,
      true_negative,
      precision,
      recall,
      specificity,
      f1,
    })
    macro_precision = macro_precision + precision
    macro_recall = macro_recall + recall
    macro_f1 = macro_f1 + f1
    let weight = support.to_double() / total.to_double()
    weighted_precision = weighted_precision + weight * precision
    weighted_recall = weighted_recall + weight * recall
    weighted_f1 = weighted_f1 + weight * f1
  }
  Ok({
    accuracy: correct.to_double() / total.to_double(),
    macro_precision: macro_precision / class_count.to_double(),
    macro_recall: macro_recall / class_count.to_double(),
    macro_f1: macro_f1 / class_count.to_double(),
    weighted_precision,
    weighted_recall,
    weighted_f1,
    classes,
  })
}

///|
fn validate_probability_rows(
  actual : Array[Int],
  probabilities : Array[Array[Double]],
  class_count : Int,
) -> Result[Unit, TreeError] {
  if class_count <= 0 {
    return Err(InvalidClassCount(class_count))
  }
  if actual.length() != probabilities.length() {
    return Err(
      ProbabilityRowCountMismatch(actual.length(), probabilities.length()),
    )
  }
  if actual.is_empty() {
    return Err(EmptyDataset)
  }
  for row_index = 0
      row_index < probabilities.length()
      row_index = row_index + 1 {
    let label = actual[row_index]
    if label < 0 || label >= class_count {
      return Err(InvalidClassLabel(row_index, label, class_count))
    }
    let row = probabilities[row_index]
    if row.length() != class_count {
      return Err(InvalidProbabilityWidth(row_index, row.length(), class_count))
    }
    let mut sum = 0.0
    for class_index = 0
        class_index < class_count
        class_index = class_index + 1 {
      let probability = row[class_index]
      if !finite_number(probability) || probability < 0.0 || probability > 1.0 {
        return Err(InvalidProbabilityValue(row_index, class_index))
      }
      sum = sum + probability
    }
    let difference = if sum < 1.0 { 1.0 - sum } else { sum - 1.0 }
    if difference > 0.000000001 {
      return Err(InvalidProbabilitySum(row_index))
    }
  }
  Ok(())
}

///|
/// Returns the mean summed squared probability error per row.
pub fn multiclass_brier_score(
  actual : Array[Int],
  probabilities : Array[Array[Double]],
  class_count : Int,
) -> Result[Double, TreeError] {
  match validate_probability_rows(actual, probabilities, class_count) {
    Err(error) => return Err(error)
    Ok(_) => ()
  }
  let mut total = 0.0
  for row_index = 0; row_index < actual.length(); row_index = row_index + 1 {
    for class_index = 0
        class_index < class_count
        class_index = class_index + 1 {
      let expected = if class_index == actual[row_index] { 1.0 } else { 0.0 }
      let difference = probabilities[row_index][class_index] - expected
      total = total + difference * difference
    }
  }
  Ok(total / actual.length().to_double())
}

///|
/// Computes natural-log loss with probabilities clamped at 1e-15.
pub fn multiclass_log_loss(
  actual : Array[Int],
  probabilities : Array[Array[Double]],
  class_count : Int,
) -> Result[Double, TreeError] {
  match validate_probability_rows(actual, probabilities, class_count) {
    Err(error) => return Err(error)
    Ok(_) => ()
  }
  let mut total = 0.0
  for row_index = 0; row_index < actual.length(); row_index = row_index + 1 {
    let raw = probabilities[row_index][actual[row_index]]
    let probability = if raw < 0.000000000000001 {
      0.000000000000001
    } else {
      raw
    }
    total = total - log2_positive(probability) * 0.6931471805599453
  }
  Ok(total / actual.length().to_double())
}

///|
fn maximum_absolute_error(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  let mut maximum = 0.0
  for index = 0; index < actual.length(); index = index + 1 {
    let difference = actual[index] - predicted[index]
    let absolute = if difference < 0.0 { -difference } else { difference }
    if absolute > maximum {
      maximum = absolute
    }
  }
  maximum
}

///|
fn mean_signed_error(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  let mut total = 0.0
  for index = 0; index < actual.length(); index = index + 1 {
    total = total + predicted[index] - actual[index]
  }
  total / actual.length().to_double()
}

///|
fn explained_variance_score(
  actual : Array[Double],
  predicted : Array[Double],
) -> Double {
  let errors = Array::makei(actual.length(), fn(index) {
    actual[index] - predicted[index]
  })
  let actual_variance = variance_value(actual)
  let error_variance = variance_value(errors)
  if actual_variance <= 0.000000000001 {
    if error_variance <= 0.000000000001 {
      1.0
    } else {
      0.0
    }
  } else {
    1.0 - error_variance / actual_variance
  }
}

///|
pub fn regression_report(
  actual : Array[Double],
  predicted : Array[Double],
) -> Result[RegressionReport, TreeError] {
  let mse = match mean_squared_error(actual, predicted) {
    Ok(value) => value
    Err(error) => return Err(error)
  }
  let mae = mean_absolute_error(actual, predicted).unwrap()
  let rmse = root_mean_squared_error(actual, predicted).unwrap()
  let r2 = r_squared(actual, predicted).unwrap()
  Ok({
    mse,
    mae,
    rmse,
    r2,
    max_error: maximum_absolute_error(actual, predicted),
    mean_error: mean_signed_error(actual, predicted),
    explained_variance: explained_variance_score(actual, predicted),
  })
}