///|
pub fn matrix_determinant_2x2(matrix : Array[Array[Double]]) -> Double {
  if matrix.length() != 2 || matrix[0].length() != 2 || matrix[1].length() != 2 {
    abort("matrix must be 2x2")
  }
  matrix[0][0] * matrix[1][1] - matrix[0][1] * matrix[1][0]
}

///|
pub fn matrix_inverse_2x2(
  matrix : Array[Array[Double]],
) -> Array[Array[Double]] {
  let determinant = matrix_determinant_2x2(matrix)
  if determinant == 0.0 {
    abort("matrix is singular")
  }
  [
    [matrix[1][1] / determinant, -matrix[0][1] / determinant],
    [-matrix[1][0] / determinant, matrix[0][0] / determinant],
  ]
}

///|
pub fn quadratic_form(
  value : Array[Double],
  matrix : Array[Array[Double]],
) -> Double {
  let transformed = []
  for row in matrix {
    transformed.push(dot_product(row, value))
  }
  dot_product(value, transformed)
}

///|
pub fn mahalanobis_distance_2d(
  value : Array[Double],
  center : Array[Double],
  covariance_matrix : Array[Array[Double]],
) -> Double {
  if value.length() != 2 || center.length() != 2 {
    abort("2D vectors are required")
  }
  let delta = [value[0] - center[0], value[1] - center[1]]
  let inverse = matrix_inverse_2x2(covariance_matrix)
  quadratic_form(delta, inverse).sqrt()
}

///|
pub fn robust_mahalanobis_2d(data : Array[Array[Double]]) -> Array[Double] {
  let center = robust_center(data)
  let covariance_matrix = robust_covariance_matrix(data, 0.1)
  let result = []
  for row in data {
    result.push(mahalanobis_distance_2d(row, center, covariance_matrix))
  }
  result
}

///|
pub fn matrix_center(matrix : Array[Array[Double]]) -> Array[Array[Double]] {
  let centers = matrix_column_means(matrix)
  let result = []
  for row in matrix {
    let output = []
    for index = 0; index < row.length(); index = index + 1 {
      output.push(row[index] - centers[index])
    }
    result.push(output)
  }
  result
}

///|
pub fn covariance_from_centered(
  centered : Array[Array[Double]],
) -> Array[Array[Double]] {
  if centered.length() == 0 {
    return []
  }
  let transposed = matrix_transpose(centered)
  let product = matrix_multiply(transposed, centered)
  let denominator = if centered.length() <= 1 {
    1.0
  } else {
    (centered.length() - 1).to_double()
  }
  matrix_scale(product, 1.0 / denominator)
}

///|
pub fn robust_covariance_from_centered(
  data : Array[Array[Double]],
) -> Array[Array[Double]] {
  let centered = standardize_columns(data)
  covariance_from_centered(centered)
}

///|
pub fn matrix_max_abs_difference(
  left : Array[Array[Double]],
  right : Array[Array[Double]],
) -> Double {
  if left.length() != right.length() {
    return 0.0
  }
  let mut maximum = 0.0
  for row = 0; row < left.length(); row = row + 1 {
    if left[row].length() != right[row].length() {
      return 0.0
    }
    for column_index = 0
        column_index < left[row].length()
        column_index = column_index + 1 {
      let difference = abs_double(
        left[row][column_index] - right[row][column_index],
      )
      if difference > maximum {
        maximum = difference
      }
    }
  }
  maximum
}

///|
pub fn matrix_clip_diagonal(
  matrix : Array[Array[Double]],
  minimum : Double,
) -> Array[Array[Double]] {
  if minimum <= 0.0 {
    abort("minimum must be positive")
  }
  let result = []
  for row = 0; row < matrix.length(); row = row + 1 {
    let output = []
    for column_index = 0
        column_index < matrix[row].length()
        column_index = column_index + 1 {
      if row == column_index {
        output.push(
          if matrix[row][column_index] < minimum {
            minimum
          } else {
            matrix[row][column_index]
          },
        )
      } else {
        output.push(matrix[row][column_index])
      }
    }
    result.push(output)
  }
  result
}

///|
pub fn matrix_to_rows(matrix : Array[Array[Double]]) -> Array[Array[Double]] {
  copy_matrix(matrix)
}

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

///|
pub fn matrix_identity(size : Int) -> Array[Array[Double]] {
  if size < 0 {
    abort("size must not be negative")
  }
  diagonal_matrix(ones(size))
}

///|
fn ones(size : Int) -> Array[Double] {
  let result = []
  for index = 0; index < size; index = index + 1 {
    result.push(1.0)
  }
  result
}

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

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

///|
pub fn robust_projection(
  data : Array[Array[Double]],
  direction : Array[Double],
) -> Array[Double] {
  let result = []
  for row in data {
    result.push(dot_product(row, direction))
  }
  result
}

///|
pub fn projection_outlier_indices(
  data : Array[Array[Double]],
  direction : Array[Double],
  threshold : Double,
) -> Array[Int] {
  outlier_indices_z(robust_projection(data, direction), threshold~)
}