///|
pub struct AnomalyDetector {
  mut center : Double
  mut scale : Double
  threshold : Double
  mut fitted : Bool
}

///|
pub fn AnomalyDetector::new(threshold? : Double = 3.5) -> AnomalyDetector {
  if threshold <= 0.0 {
    abort("threshold must be positive")
  }
  { center: 0.0, scale: 0.0, threshold, fitted: false }
}

///|
pub fn AnomalyDetector::fit(
  self : AnomalyDetector,
  data : Array[Double],
) -> Unit {
  self.center = median(data)
  self.scale = mad(data)
  self.fitted = data.length() > 0
}

///|
pub fn AnomalyDetector::score(self : AnomalyDetector, value : Double) -> Double {
  robust_z_score(value, self.center, self.scale)
}

///|
pub fn AnomalyDetector::is_anomaly(
  self : AnomalyDetector,
  value : Double,
) -> Bool {
  abs_double(self.score(value)) > self.threshold
}

///|
pub fn AnomalyDetector::score_many(
  self : AnomalyDetector,
  data : Array[Double],
) -> Array[Double] {
  let result = []
  for value in data {
    result.push(self.score(value))
  }
  result
}

///|
pub fn AnomalyDetector::flags(
  self : AnomalyDetector,
  data : Array[Double],
) -> Array[Bool] {
  let result = []
  for value in data {
    result.push(self.is_anomaly(value))
  }
  result
}

///|
pub fn AnomalyDetector::refit_without_anomalies(
  self : AnomalyDetector,
  data : Array[Double],
) -> Unit {
  let clean = []
  for value in data {
    if !self.is_anomaly(value) {
      clean.push(value)
    }
  }
  if clean.length() > 0 {
    self.fit(clean)
  }
}

///|
pub struct ConfusionMatrix {
  true_positive : Int
  false_positive : Int
  true_negative : Int
  false_negative : Int
}

///|
pub fn confusion_matrix(
  actual : Array[Bool],
  predicted : Array[Bool],
) -> ConfusionMatrix {
  if actual.length() != predicted.length() {
    return {
      true_positive: 0,
      false_positive: 0,
      true_negative: 0,
      false_negative: 0,
    }
  }
  let mut true_positive = 0
  let mut false_positive = 0
  let mut true_negative = 0
  let mut false_negative = 0
  for index = 0; index < actual.length(); index = index + 1 {
    if actual[index] && predicted[index] {
      true_positive += 1
    } else if !actual[index] && predicted[index] {
      false_positive += 1
    } else if !actual[index] && !predicted[index] {
      true_negative += 1
    } else {
      false_negative += 1
    }
  }
  { true_positive, false_positive, true_negative, false_negative }
}

///|
pub fn precision(matrix : ConfusionMatrix) -> Double {
  let denominator = matrix.true_positive + matrix.false_positive
  if denominator == 0 {
    0.0
  } else {
    matrix.true_positive.to_double() / denominator.to_double()
  }
}

///|
pub fn recall(matrix : ConfusionMatrix) -> Double {
  let denominator = matrix.true_positive + matrix.false_negative
  if denominator == 0 {
    0.0
  } else {
    matrix.true_positive.to_double() / denominator.to_double()
  }
}

///|
pub fn specificity(matrix : ConfusionMatrix) -> Double {
  let denominator = matrix.true_negative + matrix.false_positive
  if denominator == 0 {
    0.0
  } else {
    matrix.true_negative.to_double() / denominator.to_double()
  }
}

///|
pub fn f1_score(matrix : ConfusionMatrix) -> Double {
  let p = precision(matrix)
  let r = recall(matrix)
  if p + r == 0.0 {
    0.0
  } else {
    2.0 * p * r / (p + r)
  }
}

///|
pub fn balanced_accuracy(matrix : ConfusionMatrix) -> Double {
  (recall(matrix) + specificity(matrix)) / 2.0
}

///|
pub fn matthews_correlation(matrix : ConfusionMatrix) -> Double {
  let numerator = (matrix.true_positive * matrix.true_negative -
  matrix.false_positive * matrix.false_negative).to_double()
  let left = (matrix.true_positive + matrix.false_positive).to_double() *
    (matrix.true_positive + matrix.false_negative).to_double()
  let right = (matrix.true_negative + matrix.false_positive).to_double() *
    (matrix.true_negative + matrix.false_negative).to_double()
  let denominator = (left * right).sqrt()
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
pub fn threshold_sweep(
  scores : Array[Double],
  actual : Array[Bool],
  thresholds : Array[Double],
) -> Array[Array[Double]] {
  if scores.length() != actual.length() {
    return []
  }
  let result = []
  for threshold in thresholds {
    let predicted = []
    for score in scores {
      predicted.push(abs_double(score) > threshold)
    }
    let matrix = confusion_matrix(actual, predicted)
    result.push([threshold, precision(matrix), recall(matrix), f1_score(matrix)])
  }
  result
}

///|
pub fn best_f1_threshold(
  scores : Array[Double],
  actual : Array[Bool],
  thresholds : Array[Double],
) -> Double {
  let sweep = threshold_sweep(scores, actual, thresholds)
  if sweep.length() == 0 {
    return 0.0
  }
  let mut best = sweep[0][0]
  let mut best_score = sweep[0][3]
  for row in sweep {
    if row[3] > best_score {
      best = row[0]
      best_score = row[3]
    }
  }
  best
}

///|
pub fn precision_recall_curve(
  scores : Array[Double],
  actual : Array[Bool],
  thresholds : Array[Double],
) -> Array[Array[Double]] {
  let sweep = threshold_sweep(scores, actual, thresholds)
  let result = []
  for row in sweep {
    result.push([row[2], row[1], row[0]])
  }
  result
}

///|
pub fn anomaly_rate(flags : Array[Bool]) -> Double {
  if flags.length() == 0 {
    return 0.0
  }
  let mut count = 0
  for flag in flags {
    if flag {
      count += 1
    }
  }
  count.to_double() / flags.length().to_double()
}

///|
pub fn top_anomaly_indices(scores : Array[Double], count : Int) -> Array[Int] {
  if count < 0 {
    abort("count must not be negative")
  }
  let indices = []
  for index = 0; index < scores.length(); index = index + 1 {
    indices.push((index, abs_double(scores[index])))
  }
  indices.sort_by((left, right) => {
    if left.1 > right.1 {
      -1
    } else if left.1 < right.1 {
      1
    } else {
      0
    }
  })
  let result = []
  let limit = if count > indices.length() { indices.length() } else { count }
  for index = 0; index < limit; index = index + 1 {
    result.push(indices[index].0)
  }
  result
}