///|
pub fn matrix(rows~ : Int, cols~ : Int, values~ : Array[Double]) -> Matrix {
  { rows: Int::max(0, rows), cols: Int::max(0, cols), values }
}

///|
pub fn matrix_zeros(rows : Int, cols : Int) -> Matrix {
  let output : Array[Double] = []
  for _ in 0.. Matrix {
  let safe = Int::max(0, size)
  let output = matrix_zeros(safe, safe)
  for i in 0.. Bool {
  value.rows >= 0 &&
  value.cols >= 0 &&
  value.values.length() == value.rows * value.cols &&
  value.values.all(is_finite)
}

///|
fn matrix_index(value : Matrix, row : Int, col : Int) -> Int? {
  if row < 0 || col < 0 || row >= value.rows || col >= value.cols {
    None
  } else {
    Some(row * value.cols + col)
  }
}

///|
pub fn Matrix::get(self : Matrix, row : Int, col : Int) -> Double? {
  match matrix_index(self, row, col) {
    Some(index) => Some(self.values[index])
    None => None
  }
}

///|
pub fn Matrix::set(self : Matrix, row : Int, col : Int, value : Double) -> Bool {
  match matrix_index(self, row, col) {
    Some(index) => {
      self.values[index] = value
      true
    }
    None => false
  }
}

///|
pub fn Matrix::row(self : Matrix, row : Int) -> Array[Double] {
  let output : Array[Double] = []
  match matrix_index(self, row, 0) {
    Some(_) =>
      for col in 0.. ()
  }
  output
}

///|
pub fn Matrix::column(self : Matrix, col : Int) -> Array[Double] {
  let output : Array[Double] = []
  if col >= 0 && col < self.cols {
    for row in 0.. Matrix {
  let output = matrix_zeros(self.cols, self.rows)
  for row in 0.. Double) -> Matrix {
  { ..self, values: self.values.map(transform) }
}

///|
pub fn matrix_add(left : Matrix, right : Matrix) -> Matrix? {
  if left.rows != right.rows || left.cols != right.cols {
    None
  } else {
    let values : Array[Double] = []
    for i in 0.. Matrix? {
  if left.rows != right.rows || left.cols != right.cols {
    None
  } else {
    let values : Array[Double] = []
    for i in 0.. Matrix {
  self.map(fn(value) { value * factor })
}

///|
pub fn matrix_multiply(left : Matrix, right : Matrix) -> Matrix? {
  if left.cols != right.rows {
    None
  } else {
    let output = matrix_zeros(left.rows, right.cols)
    for row in 0.. Array[Double] {
  let output : Array[Double] = []
  for row in 0.. Double {
  let mut total = 0.0
  for i in 0.. Double {
  let mut total = 0.0
  for item in value.values {
    if is_finite(item) {
      total += item * item
    }
  }
  total.sqrt()
}

///|
pub fn matrix_determinant(value : Matrix) -> Double {
  if value.rows != value.cols || value.rows == 0 {
    0.0
  } else if value.rows == 1 {
    value.values[0]
  } else if value.rows == 2 {
    value.values[0] * value.values[3] - value.values[1] * value.values[2]
  } else {
    let mut total = 0.0
    for col in 0.. Matrix {
  let output : Array[Double] = []
  for row in 0.. Array[Double]? {
  if left.rows != left.cols || right.length() != left.rows {
    None
  } else {
    let n = left.rows
    let a = left.values.copy()
    let b = right.to_owned()
    let mut valid = true
    for pivot in 0.. pivot_size {
          pivot_size = candidate
          pivot_row = row
        }
      }
      if pivot_size < 1.0e-12 || !is_finite(pivot_size) {
        valid = false
      } else {
        if pivot_row != pivot {
          for col in pivot..= 0 {
        let mut value = b[row]
        for col in (row + 1).. Matrix {
  let safe_degree = Int::max(0, degree)
  let values : Array[Double] = []
  for coordinate in coordinates {
    let mut power = 1.0
    for _ in 0..<=safe_degree {
      values.push(power)
      power *= coordinate
    }
  }
  { rows: coordinates.length(), cols: safe_degree + 1, values }
}

///|
pub fn polynomial_evaluate(
  coefficients : ArrayView[Double],
  coordinate : Double,
) -> Double {
  let mut result = 0.0
  for coefficient in coefficients.rev() {
    result = result * coordinate + coefficient
  }
  result
}

///|
pub fn polynomial_fit(
  coordinates : ArrayView[Double],
  values : ArrayView[Double],
  degree : Int,
) -> LeastSquaresResult {
  let count = Int::min(coordinates.length(), values.length())
  let design = polynomial_design(coordinates[:count], degree)
  let transpose = design.transpose()
  let normal = matrix_multiply(transpose, design)
  let rhs_matrix = matrix(
    rows=transpose.rows,
    cols=1,
    values=transpose.multiply_vector(values[:count]),
  )
  let rhs = rhs_matrix.column(0)
  let coefficients = match normal {
    Some(system) => matrix_solve(system, rhs).unwrap_or([])
    None => []
  }
  let mut residual_sum_squares = 0.0
  for i in 0.. Matrix {
  let output = matrix_zeros(values.length(), values.length())
  let mean = if values.length() == 0 {
    0.0
  } else {
    sum_finite(values) / Double::from_int(values.length())
  }
  for row in 0..