///|
/// Dense dynamic M x N matrix with double precision.
pub(all) struct DenseMatrix {
  rows : Int
  cols : Int
  data : Array[Double]
} derive(Debug, Eq)

///|
/// Creates an M x N matrix initialized to zero.
pub fn DenseMatrix::zeros(rows : Int, cols : Int) -> DenseMatrix {
  let size = rows * cols
  let data = Array::make(size, 0.0)
  { rows, cols, data }
}

///|
/// Creates an N x N identity matrix.
pub fn DenseMatrix::identity(n : Int) -> DenseMatrix {
  let mat = DenseMatrix::zeros(n, n)
  for i in 0.. DenseMatrix {
  let r = rows_arr.length()
  if r == 0 {
    return DenseMatrix::zeros(0, 0)
  }
  let c = rows_arr[0].length()
  let mat = DenseMatrix::zeros(r, c)
  for i in 0.. DenseMatrix {
  let expected = rows * cols
  let actual_data = if data.length() == expected {
    data.copy()
  } else {
    let d = Array::make(expected, 0.0)
    let min_len = if data.length() < expected {
      data.length()
    } else {
      expected
    }
    for i in 0.. Double {
  if row < 0 || row >= self.rows || col < 0 || col >= self.cols {
    return 0.0
  }
  self.data[row * self.cols + col]
}

///|
/// Sets element at (row, col).
pub fn DenseMatrix::set(
  self : DenseMatrix,
  row : Int,
  col : Int,
  val : Double,
) -> Unit {
  if row >= 0 && row < self.rows && col >= 0 && col < self.cols {
    self.data[row * self.cols + col] = val
  }
}

///|
/// Clones the matrix.
pub fn DenseMatrix::clone(self : DenseMatrix) -> DenseMatrix {
  { rows: self.rows, cols: self.cols, data: self.data.copy() }
}

///|
/// Transposes the matrix into an N x M matrix.
pub fn DenseMatrix::transpose(self : DenseMatrix) -> DenseMatrix {
  let res = DenseMatrix::zeros(self.cols, self.rows)
  for i in 0.. DenseMatrix? {
  if self.rows != other.rows || self.cols != other.cols {
    return None
  }
  let res = DenseMatrix::zeros(self.rows, self.cols)
  let len = self.data.length()
  for i in 0.. DenseMatrix? {
  if self.rows != other.rows || self.cols != other.cols {
    return None
  }
  let res = DenseMatrix::zeros(self.rows, self.cols)
  let len = self.data.length()
  for i in 0.. DenseMatrix {
  let res = DenseMatrix::zeros(self.rows, self.cols)
  let len = self.data.length()
  for i in 0.. DenseMatrix? {
  if self.cols != other.rows {
    return None
  }
  let res = DenseMatrix::zeros(self.rows, other.cols)
  for i in 0.. 1.0e-15 {
        for j in 0.. Array[Double]? {
  if self.cols != vec.length() {
    return None
  }
  let res = Array::make(self.rows, 0.0)
  for i in 0.. Double {
  let mut sum = 0.0
  for v in self.data {
    sum += v * v
  }
  sum.sqrt()
}

///|
/// Extracts a submatrix from (r0, c0) to (r1-1, c1-1).
pub fn DenseMatrix::submatrix(
  self : DenseMatrix,
  r0 : Int,
  c0 : Int,
  r1 : Int,
  c1 : Int,
) -> DenseMatrix {
  let r = if r1 > r0 { r1 - r0 } else { 0 }
  let c = if c1 > c0 { c1 - c0 } else { 0 }
  let res = DenseMatrix::zeros(r, c)
  for i in 0.. Unit {
  for i in 0.. Unit {
  if r1 == r2 {
    return
  }
  for j in 0.. Array[Double]? {
  if self.rows != self.cols || self.rows != b.length() {
    return None
  }
  let n = self.rows
  let a = self.clone()
  let x = b.copy()

  // Forward elimination
  for k in 0.. max_val {
        max_val = val
        max_row = i
      }
    }

    if max_val < 1.0e-12 {
      return None // Singular matrix
    }

    // Swap row k and max_row
    if max_row != k {
      a.swap_rows(k, max_row)
      let tmp = x[k]
      x[k] = x[max_row]
      x[max_row] = tmp
    }

    // Eliminate below
    let pivot = a.get(k, k)
    for i in (k + 1)..= 0 {
    let mut sum = x[i]
    for j in (i + 1).. DenseMatrix? {
  if self.rows != self.cols {
    return None
  }
  let n = self.rows
  let a = self.clone()
  let inv = DenseMatrix::identity(n)

  for k in 0.. max_val {
        max_val = val
        max_row = i
      }
    }

    if max_val < 1.0e-12 {
      return None
    }

    if max_row != k {
      a.swap_rows(k, max_row)
      inv.swap_rows(k, max_row)
    }

    let pivot = a.get(k, k)
    for j in 0.. Double {
  if self.rows != self.cols {
    return 0.0
  }
  let n = self.rows
  let a = self.clone()
  let mut det = 1.0

  for k in 0.. max_val {
        max_val = val
        max_row = i
      }
    }

    if max_val < 1.0e-14 {
      return 0.0
    }

    if max_row != k {
      a.swap_rows(k, max_row)
      det = -det
    }

    let pivot = a.get(k, k)
    det *= pivot

    for i in (k + 1).. DenseMatrix? {
  if self.rows != self.cols {
    return None
  }
  let n = self.rows
  let l = DenseMatrix::zeros(n, n)

  for i in 0.. (DenseMatrix, DenseMatrix)? {
  let m = self.rows
  let n = self.cols
  if m < n {
    return None
  }
  let q = DenseMatrix::identity(m)
  let r = self.clone()

  for k in 0.. 1.0e-14 {
      let x0 = r.get(k, k)
      let alpha = if x0 >= 0.0 { -norm } else { norm }
      let u0 = x0 - alpha
      let mut u_norm_sq = u0 * u0
      for i in (k + 1).. 1.0e-14 {
        // v is normalized Householder vector
        let v = Array::make(len, 0.0)
        v[0] = u0 / u_norm
        for i in 1..= N) using One-Sided Jacobi rotations.
/// Returns (U, S, V) such that A = U * diag(S) * V^T.
pub fn DenseMatrix::svd(
  self : DenseMatrix,
  max_iterations? : Int = 50,
  epsilon? : Double = 1.0e-12,
) -> (DenseMatrix, Array[Double], DenseMatrix)? {
  let m = self.rows
  let n = self.cols
  if m < n {
    // If m < n, compute SVD of A^T
    let at = self.transpose()
    match at.svd(max_iterations~, epsilon~) {
      Some((u, s, v)) => return Some((v, s, u))
      None => return None
    }
  }

  let a = self.clone()
  let v = DenseMatrix::identity(n)

  // One-sided Jacobi algorithm
  let mut iter = 0
  let mut converged = false
  while iter < max_iterations && !converged {
    converged = true
    for p in 0..<(n - 1) {
      for q in (p + 1).. epsilon * (app * aqq).sqrt() && apq.abs() > epsilon {
          converged = false
          let tau = (aqq - app) / (2.0 * apq)
          let t = if tau >= 0.0 {
            1.0 / (tau + (1.0 + tau * tau).sqrt())
          } else {
            -1.0 / (-tau + (1.0 + tau * tau).sqrt())
          }
          let c = 1.0 / (1.0 + t * t).sqrt()
          let s = t * c

          // Rotate columns p and q of A
          for i in 0.. 1.0e-14 {
      for i in 0.. s_vals[max_idx] {
        max_idx = k
      }
    }
    if max_idx != i {
      // Swap singular values
      let tmp_s = s_vals[i]
      s_vals[i] = s_vals[max_idx]
      s_vals[max_idx] = tmp_s

      // Swap columns of U
      for r in 0.. DenseMatrix? {
  match self.svd() {
    Some((u, s, v)) => {
      let m = self.rows
      let n = self.cols
      let pinv = DenseMatrix::zeros(n, m)
      let k = s.length()

      // V * S^+ * U^T
      for i in 0.. tolerance {
              sum += v.get(i, l) * (1.0 / sigma) * u.get(j, l)
            }
          }
          pinv.set(i, j, sum)
        }
      }
      Some(pinv)
    }
    None => None
  }
}

///|
/// Solves linear least squares problem min ||A * x - b||_2 using SVD / Pseudo-Inverse.
pub fn DenseMatrix::solve_least_squares(
  self : DenseMatrix,
  b : Array[Double],
) -> Array[Double]? {
  match self.pseudo_inverse() {
    Some(pinv) => pinv.mul_vec(b)
    None => None
  }
}