///|
/// Summary of a numeric matrix and its conditioning.
pub struct MatrixProfile {
  rows : Int
  columns : Int
  rank_proxy : Int
  frobenius_norm : Double
  maximum_abs : Double
  trace : Double
  symmetry_error : Double
  condition_proxy : Double
}

///|
/// One principal component estimated with robust power iteration.
pub struct PrincipalComponent {
  direction : Array[Double]
  eigenvalue : Double
  explained_ratio : Double
  iterations : Int
  converged : Bool
}

///|
pub fn matrix_row_count(matrix : Array[Array[Double]]) -> Int {
  matrix.length()
}

///|
pub fn matrix_column_count(matrix : Array[Array[Double]]) -> Int {
  if matrix.length() == 0 {
    0
  } else {
    matrix[0].length()
  }
}

///|
pub fn matrix_rectangular(matrix : Array[Array[Double]]) -> Bool {
  let columns = matrix_column_count(matrix)
  for row in matrix {
    if row.length() != columns {
      return false
    }
  }
  true
}

///|
pub fn matrix_non_empty(matrix : Array[Array[Double]]) -> Bool {
  matrix.length() > 0 &&
  matrix_column_count(matrix) > 0 &&
  matrix_rectangular(matrix)
}

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

///|
pub fn matrix_row_or_empty(
  matrix : Array[Array[Double]],
  index : Int,
) -> Array[Double] {
  if index < 0 || index >= matrix.length() {
    []
  } else {
    matrix[index].copy()
  }
}

///|
pub fn matrix_column_or_empty(
  matrix : Array[Array[Double]],
  index : Int,
) -> Array[Double] {
  let result = []
  if index < 0 {
    return result
  }
  for row in matrix {
    if index < row.length() {
      result.push(row[index])
    }
  }
  result
}

///|
pub fn matrix_fill(
  rows : Int,
  columns : Int,
  value : Double,
) -> Array[Array[Double]] {
  let result = []
  let safe_rows = if rows < 0 { 0 } else { rows }
  let safe_columns = if columns < 0 { 0 } else { columns }
  for _ in 0.. Array[Array[Double]] {
  let result = matrix_copy(matrix)
  let limit = if result.length() < matrix_column_count(result) {
    result.length()
  } else {
    matrix_column_count(result)
  }
  for index = 0; index < limit; index = index + 1 {
    result[index][index] += value
  }
  result
}

///|
pub fn matrix_diagonal_values(matrix : Array[Array[Double]]) -> Array[Double] {
  let result = []
  let limit = if matrix.length() < matrix_column_count(matrix) {
    matrix.length()
  } else {
    matrix_column_count(matrix)
  }
  for index = 0; index < limit; index = index + 1 {
    result.push(matrix[index][index])
  }
  result
}

///|
pub fn matrix_maximum_abs(matrix : Array[Array[Double]]) -> Double {
  let mut result = 0.0
  for row in matrix {
    for value in row {
      if abs_double(value) > result {
        result = abs_double(value)
      }
    }
  }
  result
}

///|
pub fn matrix_symmetry_error(matrix : Array[Array[Double]]) -> Double {
  if !matrix_rectangular(matrix) ||
    matrix.length() != matrix_column_count(matrix) {
    return 1.0e12
  }
  let mut result = 0.0
  for row = 0; row < matrix.length(); row = row + 1 {
    for column = row + 1; column < matrix.length(); column = column + 1 {
      let error = abs_double(matrix[row][column] - matrix[column][row])
      if error > result {
        result = error
      }
    }
  }
  result
}

///|
pub fn matrix_profile(matrix : Array[Array[Double]]) -> MatrixProfile {
  let rows = matrix_row_count(matrix)
  let columns = matrix_column_count(matrix)
  let diagonal = matrix_diagonal_values(matrix)
  let nonzero = []
  for row in matrix {
    for value in row {
      if abs_double(value) > 1.0e-12 {
        nonzero.push(value)
      }
    }
  }
  let minimum = if nonzero.length() == 0 {
    0.0
  } else {
    abs_double(min_value(nonzero))
  }
  let maximum = matrix_maximum_abs(matrix)
  {
    rows,
    columns,
    rank_proxy: if nonzero.length() > rows {
      rows
    } else {
      nonzero.length()
    },
    frobenius_norm: matrix_frobenius_norm(matrix),
    maximum_abs: maximum,
    trace: sum_values(diagonal),
    symmetry_error: matrix_symmetry_error(matrix),
    condition_proxy: if minimum <= 1.0e-12 {
      1.0e12
    } else {
      maximum / minimum
    },
  }
}

///|
pub fn matrix_profile_vector(profile : MatrixProfile) -> Array[Double] {
  [
    profile.rows.to_double(),
    profile.columns.to_double(),
    profile.rank_proxy.to_double(),
    profile.frobenius_norm,
    profile.maximum_abs,
    profile.trace,
    profile.symmetry_error,
    profile.condition_proxy,
  ]
}

///|
pub fn matrix_center_rows(
  matrix : Array[Array[Double]],
) -> Array[Array[Double]] {
  let result = matrix_copy(matrix)
  for row = 0; row < result.length(); row = row + 1 {
    let center = median(result[row])
    for column = 0; column < result[row].length(); column = column + 1 {
      result[row][column] -= center
    }
  }
  result
}

///|
pub fn matrix_center_columns(
  matrix : Array[Array[Double]],
) -> Array[Array[Double]] {
  if !matrix_rectangular(matrix) {
    return []
  }
  let result = matrix_copy(matrix)
  let columns = matrix_column_count(result)
  for column = 0; column < columns; column = column + 1 {
    let values = matrix_column_or_empty(result, column)
    let center = median(values)
    for row = 0; row < result.length(); row = row + 1 {
      result[row][column] -= center
    }
  }
  result
}

///|
pub fn matrix_scale_rows(matrix : Array[Array[Double]]) -> Array[Array[Double]] {
  let result = matrix_copy(matrix)
  for row = 0; row < result.length(); row = row + 1 {
    let scale = mad(result[row]) * 1.4826
    let denominator = if scale <= 1.0e-12 { 1.0 } else { scale }
    for column = 0; column < result[row].length(); column = column + 1 {
      result[row][column] /= denominator
    }
  }
  result
}

///|
pub fn matrix_scale_columns(
  matrix : Array[Array[Double]],
) -> Array[Array[Double]] {
  if !matrix_rectangular(matrix) {
    return []
  }
  let result = matrix_copy(matrix)
  let columns = matrix_column_count(result)
  for column = 0; column < columns; column = column + 1 {
    let values = matrix_column_or_empty(result, column)
    let scale = mad(values) * 1.4826
    let denominator = if scale <= 1.0e-12 { 1.0 } else { scale }
    for row = 0; row < result.length(); row = row + 1 {
      result[row][column] /= denominator
    }
  }
  result
}

///|
pub fn matrix_robust_standardize(
  matrix : Array[Array[Double]],
) -> Array[Array[Double]] {
  matrix_scale_columns(matrix_center_columns(matrix))
}

///|
pub fn matrix_row_dot(left : Array[Double], right : Array[Double]) -> Double {
  let limit = if left.length() < right.length() {
    left.length()
  } else {
    right.length()
  }
  let mut result = 0.0
  for index = 0; index < limit; index = index + 1 {
    result += left[index] * right[index]
  }
  result
}

///|
pub fn matrix_row_norm(row : Array[Double]) -> Double {
  sum_squared(row).sqrt()
}

///|
pub fn matrix_row_normalize(row : Array[Double]) -> Array[Double] {
  let norm = matrix_row_norm(row)
  if norm <= 1.0e-12 {
    row.copy()
  } else {
    let result = []
    for value in row {
      result.push(value / norm)
    }
    result
  }
}

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

///|
pub fn matrix_transpose_vector_product(
  matrix : Array[Array[Double]],
  vector : Array[Double],
) -> Array[Double] {
  let columns = matrix_column_count(matrix)
  let result = []
  for column = 0; column < columns; column = column + 1 {
    let values = matrix_column_or_empty(matrix, column)
    result.push(matrix_row_dot(values, vector))
  }
  result
}

///|
pub fn matrix_gram(matrix : Array[Array[Double]]) -> Array[Array[Double]] {
  if !matrix_rectangular(matrix) {
    return []
  }
  let columns = matrix_column_count(matrix)
  let result = matrix_fill(columns, columns, 0.0)
  for left = 0; left < columns; left = left + 1 {
    let left_column = matrix_column_or_empty(matrix, left)
    for right = 0; right < columns; right = right + 1 {
      result[left][right] = matrix_row_dot(
        left_column,
        matrix_column_or_empty(matrix, right),
      )
    }
  }
  result
}

///|
pub fn matrix_covariance(matrix : Array[Array[Double]]) -> Array[Array[Double]] {
  covariance_from_centered(matrix_center_columns(matrix))
}

///|
pub fn matrix_robust_covariance(
  matrix : Array[Array[Double]],
  tuning : Double,
) -> Array[Array[Double]] {
  let safe_tuning = if tuning <= 0.0 { 1.0 } else { tuning }
  robust_covariance_from_centered(
    matrix_scale(
      matrix_scale_columns(matrix_center_columns(matrix)),
      safe_tuning,
    ),
  )
}

///|
pub fn matrix_project(
  matrix : Array[Array[Double]],
  direction : Array[Double],
) -> Array[Double] {
  matrix_rows_dot(matrix, direction)
}

///|
pub fn matrix_reconstruct_rank_one(
  scores : Array[Double],
  direction : Array[Double],
) -> Array[Array[Double]] {
  let result = []
  for score in scores {
    let row = []
    for value in direction {
      row.push(score * value)
    }
    result.push(row)
  }
  result
}

///|
pub fn matrix_residual_after_projection(
  matrix : Array[Array[Double]],
  direction : Array[Double],
) -> Array[Array[Double]] {
  let projection = matrix_reconstruct_rank_one(
    matrix_project(matrix, direction),
    direction,
  )
  matrix_add(matrix, matrix_scale(projection, -1.0))
}

///|
pub fn matrix_power_component(
  covariance : Array[Array[Double]],
  max_iter : Int,
  tolerance : Double,
) -> PrincipalComponent {
  let columns = matrix_column_count(covariance)
  if columns == 0 || !matrix_rectangular(covariance) {
    return {
      direction: [],
      eigenvalue: 0.0,
      explained_ratio: 0.0,
      iterations: 0,
      converged: false,
    }
  }
  let direction = []
  for index = 0; index < columns; index = index + 1 {
    direction.push(1.0 / columns.to_double().sqrt())
  }
  let mut vector = direction
  let mut eigenvalue = 0.0
  let mut converged = false
  let iterations = if max_iter < 1 { 1 } else { max_iter }
  let limit = if tolerance <= 0.0 { 1.0e-6 } else { tolerance }
  let mut used = 0
  for iteration = 0; iteration < iterations; iteration = iteration + 1 {
    let next = matrix_rows_dot(covariance, vector)
    let norm = matrix_row_norm(next)
    let normalized = if norm <= 1.0e-12 {
      vector.copy()
    } else {
      matrix_row_normalize(next)
    }
    let value = matrix_row_dot(
      normalized,
      matrix_rows_dot(covariance, normalized),
    )
    if abs_double(value - eigenvalue) <= limit {
      converged = true
    }
    eigenvalue = value
    vector = normalized
    used += 1
  }
  let total = matrix_trace(covariance)
  {
    direction: vector,
    eigenvalue,
    explained_ratio: if abs_double(total) <= 1.0e-12 {
      0.0
    } else {
      eigenvalue / total
    },
    iterations: used,
    converged,
  }
}

///|
pub fn matrix_principal_component(
  matrix : Array[Array[Double]],
) -> PrincipalComponent {
  let covariance = matrix_covariance(matrix)
  matrix_power_component(covariance, 50, 1.0e-6)
}

///|
pub fn matrix_projection_scores(matrix : Array[Array[Double]]) -> Array[Double] {
  let component = matrix_principal_component(matrix)
  matrix_project(matrix, component.direction)
}

///|
pub fn matrix_reconstruction_error(matrix : Array[Array[Double]]) -> Double {
  let component = matrix_principal_component(matrix)
  let residual = matrix_residual_after_projection(matrix, component.direction)
  matrix_frobenius_norm(residual)
}

///|
pub fn matrix_leverage_scores(matrix : Array[Array[Double]]) -> Array[Double] {
  let gram = matrix_gram(matrix)
  let diagonal = matrix_diagonal_values(gram)
  let total = sum_values(diagonal)
  if total <= 1.0e-12 {
    let result = []
    for _ in diagonal {
      result.push(0.0)
    }
    result
  } else {
    let result = []
    for value in diagonal {
      result.push(value / total)
    }
    result
  }
}

///|
pub fn matrix_outlier_rows(
  matrix : Array[Array[Double]],
  threshold : Double,
) -> Array[Int] {
  let scores = mahalanobis_scores(matrix)
  let result = []
  for index = 0; index < scores.length(); index = index + 1 {
    if scores[index] > threshold {
      result.push(index)
    }
  }
  result
}

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

///|
pub fn matrix_quality_score(matrix : Array[Array[Double]]) -> Double {
  let scores = matrix_row_quality(matrix)
  if scores.length() == 0 {
    0.0
  } else {
    mean(scores)
  }
}