///|
pub fn column(data : Array[Array[Double]], index : Int) -> Array[Double] {
  let result = []
  for row in data {
    if index < 0 || index >= row.length() {
      abort("column index out of bounds")
    }
    result.push(row[index])
  }
  result
}

///|
pub fn robust_covariance_matrix(
  data : Array[Array[Double]],
  trim_percent : Double,
) -> Array[Array[Double]] {
  if data.length() == 0 {
    return []
  }
  let dimensions = data[0].length()
  let result = []
  for left = 0; left < dimensions; left = left + 1 {
    let row = []
    let x = column(data, left)
    for right = 0; right < dimensions; right = right + 1 {
      let y = column(data, right)
      row.push(robust_covariance(x, y, trim_percent))
    }
    result.push(row)
  }
  result
}

///|
pub fn robust_correlation_matrix(
  data : Array[Array[Double]],
  trim_percent : Double,
) -> Array[Array[Double]] {
  if data.length() == 0 {
    return []
  }
  let dimensions = data[0].length()
  let result = []
  for left = 0; left < dimensions; left = left + 1 {
    let row = []
    let x = column(data, left)
    for right = 0; right < dimensions; right = right + 1 {
      let y = column(data, right)
      row.push(
        pearson_correlation(
          winsorize(x, trim_percent),
          winsorize(y, trim_percent),
        ),
      )
    }
    result.push(row)
  }
  result
}

///|
pub fn diagonal_matrix(values : Array[Double]) -> Array[Array[Double]] {
  let result = []
  for row = 0; row < values.length(); row = row + 1 {
    let output = []
    for column_index = 0
        column_index < values.length()
        column_index = column_index + 1 {
      output.push(if row == column_index { values[row] } else { 0.0 })
    }
    result.push(output)
  }
  result
}

///|
pub fn matrix_transpose(matrix : Array[Array[Double]]) -> Array[Array[Double]] {
  if matrix.length() == 0 {
    return []
  }
  let columns = matrix[0].length()
  let result = []
  for column_index = 0; column_index < columns; column_index = column_index + 1 {
    result.push(column(matrix, column_index))
  }
  result
}

///|
pub fn matrix_multiply(
  left : Array[Array[Double]],
  right : Array[Array[Double]],
) -> Array[Array[Double]] {
  if left.length() == 0 || right.length() == 0 {
    return []
  }
  if left[0].length() != right.length() {
    abort("incompatible matrix dimensions")
  }
  let right_t = matrix_transpose(right)
  let result = []
  for row in left {
    let output = []
    for column_values in right_t {
      output.push(dot_product(row, column_values))
    }
    result.push(output)
  }
  result
}

///|
pub fn matrix_add(
  left : Array[Array[Double]],
  right : Array[Array[Double]],
) -> Array[Array[Double]] {
  if left.length() != right.length() {
    abort("incompatible matrix dimensions")
  }
  let result = []
  for row = 0; row < left.length(); row = row + 1 {
    if left[row].length() != right[row].length() {
      abort("incompatible matrix dimensions")
    }
    let output = []
    for column_index = 0
        column_index < left[row].length()
        column_index = column_index + 1 {
      output.push(left[row][column_index] + right[row][column_index])
    }
    result.push(output)
  }
  result
}

///|
pub fn matrix_scale(
  matrix : Array[Array[Double]],
  factor : Double,
) -> Array[Array[Double]] {
  let result = []
  for row in matrix {
    let output = []
    for value in row {
      output.push(value * factor)
    }
    result.push(output)
  }
  result
}

///|
pub fn matrix_trace(matrix : Array[Array[Double]]) -> Double {
  let mut total = 0.0
  for index = 0; index < matrix.length(); index = index + 1 {
    if index < matrix[index].length() {
      total += matrix[index][index]
    }
  }
  total
}

///|
pub fn mahalanobis_distance_diagonal(
  value : Array[Double],
  center : Array[Double],
  scale : Array[Double],
) -> Double {
  if value.length() != center.length() || value.length() != scale.length() {
    return 0.0
  }
  let mut total = 0.0
  for index = 0; index < value.length(); index = index + 1 {
    if scale[index] <= 0.0 {
      continue
    }
    let normalized = (value[index] - center[index]) / scale[index]
    total += normalized * normalized
  }
  total.sqrt()
}

///|
pub fn robust_center(data : Array[Array[Double]]) -> Array[Double] {
  if data.length() == 0 {
    return []
  }
  let dimensions = data[0].length()
  let result = []
  for index = 0; index < dimensions; index = index + 1 {
    result.push(median(column(data, index)))
  }
  result
}

///|
pub fn robust_scales(data : Array[Array[Double]]) -> Array[Double] {
  if data.length() == 0 {
    return []
  }
  let dimensions = data[0].length()
  let result = []
  for index = 0; index < dimensions; index = index + 1 {
    result.push(mad(column(data, index)))
  }
  result
}

///|
pub fn mahalanobis_scores(data : Array[Array[Double]]) -> Array[Double] {
  let center = robust_center(data)
  let scales = robust_scales(data)
  let result = []
  for row in data {
    result.push(mahalanobis_distance_diagonal(row, center, scales))
  }
  result
}

///|
pub fn multivariate_outlier_indices(
  data : Array[Array[Double]],
  threshold : Double,
) -> Array[Int] {
  if threshold <= 0.0 {
    abort("threshold must be positive")
  }
  let scores = mahalanobis_scores(data)
  let result = []
  for index = 0; index < scores.length(); index = index + 1 {
    if scores[index] > threshold {
      result.push(index)
    }
  }
  result
}

///|
pub fn standardize_columns(data : Array[Array[Double]]) -> Array[Array[Double]] {
  if data.length() == 0 {
    return []
  }
  let center = robust_center(data)
  let scales = robust_scales(data)
  let result = []
  for row in data {
    let output = []
    for index = 0; index < row.length(); index = index + 1 {
      let scale = if scales[index] == 0.0 { 1.0 } else { scales[index] }
      output.push((row[index] - center[index]) / scale)
    }
    result.push(output)
  }
  result
}

///|
pub fn matrix_row_means(matrix : Array[Array[Double]]) -> Array[Double] {
  let result = []
  for row in matrix {
    result.push(mean(row))
  }
  result
}

///|
pub fn matrix_column_means(matrix : Array[Array[Double]]) -> Array[Double] {
  if matrix.length() == 0 {
    return []
  }
  let result = []
  for index = 0; index < matrix[0].length(); index = index + 1 {
    result.push(mean(column(matrix, index)))
  }
  result
}

///|
pub fn matrix_frobenius_norm(matrix : Array[Array[Double]]) -> Double {
  let mut total = 0.0
  for row in matrix {
    total += sum_squared(row)
  }
  total.sqrt()
}

///|
pub fn matrix_is_symmetric(
  matrix : Array[Array[Double]],
  tolerance? : Double = 0.000001,
) -> Bool {
  for row = 0; row < matrix.length(); row = row + 1 {
    for column_index = 0
        column_index < matrix[row].length()
        column_index = column_index + 1 {
      if column_index >= matrix.length() || row >= matrix[column_index].length() {
        return false
      }
      if !nearly_equal(
          matrix[row][column_index],
          matrix[column_index][row],
          tolerance,
        ) {
        return false
      }
    }
  }
  true
}