///|
/// Sparse linear algebra storage used by the solver kernels.
///
/// The matrix is stored in compressed sparse column (CSC) form because every
/// pricing / basis update step of the revised simplex walks columns, not rows.
///
/// Invariants (enforced at construction, asserted by tests):
/// – `col_ptr` has length `ncols + 1`, starts at 0 and is non-decreasing;
/// – `row_idx` and `values` have the same length (the number of stored entries);
/// – row indices inside one column are strictly increasing;
/// – no structural zeros are stored.

///|
/// A sparse matrix in compressed sparse column form.
pub struct SparseMatrix {
  nrows : Int32
  ncols : Int32
  col_ptr : Array[Int32]
  row_idx : Array[Int32]
  values : Array[Double]
} derive(Debug)

///|
/// Builds a matrix from per-column `(row, value)` pairs.
///
/// Entries are sorted by row index, duplicate rows inside a column are summed,
/// and exact zeros are dropped. Passing a row index outside `[0, nrows)` is a
/// programming error and aborts.
pub fn SparseMatrix::from_columns(
  nrows : Int32,
  columns : Array[Array[(Int32, Double)]],
) -> SparseMatrix {
  if nrows < 0 {
    abort("SparseMatrix::from_columns: nrows must be non-negative")
  }
  let ncols = columns.length()
  let col_ptr : Array[Int32] = [0]
  let row_idx : Array[Int32] = []
  let values : Array[Double] = []
  for j = 0; j < ncols; j = j + 1 {
    let sorted = sort_entries(columns[j])
    let col_start = row_idx.length()
    for k = 0; k < sorted.length(); k = k + 1 {
      let (r, v) = sorted[k]
      if r < 0 || r >= nrows {
        abort("SparseMatrix::from_columns: row index out of range")
      }
      // Structural zeros are dropped: the stored pattern must stay exact.
      if v != 0.0 {
        if row_idx.length() > col_start && row_idx[row_idx.length() - 1] == r {
          // Same row twice inside *this* column: merge.
          let last = values.length() - 1
          values[last] = values[last] + v
        } else {
          row_idx.push(r)
          values.push(v)
        }
      }
    }
    col_ptr.push(row_idx.length())
  }
  { nrows, ncols, col_ptr, row_idx, values, }
}

///|
/// Builds an all-zero matrix with the given shape.
pub fn SparseMatrix::zeros(nrows : Int32, ncols : Int32) -> SparseMatrix {
  let col_ptr : Array[Int32] = []
  for j = 0; j <= ncols; j = j + 1 {
    col_ptr.push(0)
  }
  { nrows, ncols, col_ptr, row_idx: [], values: [], }
}

///|
/// Number of rows.
pub fn SparseMatrix::nrows(self : SparseMatrix) -> Int32 {
  self.nrows
}

///|
/// Number of columns.
pub fn SparseMatrix::ncols(self : SparseMatrix) -> Int32 {
  self.ncols
}

///|
/// Number of stored (non-zero) entries.
pub fn SparseMatrix::nnz(self : SparseMatrix) -> Int32 {
  self.values.length()
}

///|
/// Number of stored entries in column `j`.
pub fn SparseMatrix::col_nnz(self : SparseMatrix, j : Int32) -> Int32 {
  self.col_ptr[j + 1] - self.col_ptr[j]
}

///|
/// `(row, value)` pairs of column `j`, row index ascending.
pub fn SparseMatrix::column(
  self : SparseMatrix,
  j : Int32,
) -> Array[(Int32, Double)] {
  let out : Array[(Int32, Double)] = []
  for k = self.col_ptr[j]; k < self.col_ptr[j + 1]; k = k + 1 {
    out.push((self.row_idx[k], self.values[k]))
  }
  out
}

///|
/// Entry at `(i, j)`; missing entries read as `0.0`.
pub fn SparseMatrix::get(self : SparseMatrix, i : Int32, j : Int32) -> Double {
  for k = self.col_ptr[j]; k < self.col_ptr[j + 1]; k = k + 1 {
    if self.row_idx[k] == i {
      return self.values[k]
    }
  }
  0.0
}

///|
/// Raw column pointers, for kernels whose inner loops must not allocate.
///
/// With `row_indices` and `values_view` this exposes the storage so a solver can
/// walk the entries of a column in place: `column(j)` builds a fresh
/// `(row, value)` list per call, which is convenient but costs an allocation per
/// entry per pivot in a hot loop. The invariants documented at the top of this
/// file still hold for everything reached this way; the views must not outlive a
/// call that rebuilds the matrix.
pub fn SparseMatrix::col_pointers(self : SparseMatrix) -> ArrayView[Int32] {
  self.col_ptr[:]
}

///|
/// Raw row indices of the stored entries.
pub fn SparseMatrix::row_indices(self : SparseMatrix) -> ArrayView[Int32] {
  self.row_idx[:]
}

///|
/// Raw values of the stored entries, aligned with `row_indices`.
pub fn SparseMatrix::values_view(self : SparseMatrix) -> ArrayView[Double] {
  self.values[:]
}

///|
/// Transposes the matrix. Row `i` of `self` becomes column `i` of the result, so
/// the result has `ncols` rows and `nrows` columns.
pub fn SparseMatrix::transpose(self : SparseMatrix) -> SparseMatrix {
  let cols : Array[Array[(Int32, Double)]] = []
  for i = 0; i < self.nrows; i = i + 1 {
    cols.push([])
  }
  for j = 0; j < self.ncols; j = j + 1 {
    for k = self.col_ptr[j]; k < self.col_ptr[j + 1]; k = k + 1 {
      cols[self.row_idx[k]].push((j, self.values[k]))
    }
  }
  SparseMatrix::from_columns(self.ncols, cols)
}

///|
/// Row-major dense copy. Intended for tests and for the reference solver, not
/// for production paths.
pub fn SparseMatrix::to_dense(self : SparseMatrix) -> Array[Array[Double]] {
  let rows : Array[Array[Double]] = []
  for i = 0; i < self.nrows; i = i + 1 {
    let row : Array[Double] = []
    for j = 0; j < self.ncols; j = j + 1 {
      row.push(0.0)
    }
    rows.push(row)
  }
  for j = 0; j < self.ncols; j = j + 1 {
    for k = self.col_ptr[j]; k < self.col_ptr[j + 1]; k = k + 1 {
      rows[self.row_idx[k]][j] = self.values[k]
    }
  }
  rows
}

///|
/// Matrix–vector product `self * x`.
pub fn SparseMatrix::mul_vec(
  self : SparseMatrix,
  x : Array[Double],
) -> Array[Double] {
  let y : Array[Double] = []
  for i = 0; i < self.nrows; i = i + 1 {
    y.push(0.0)
  }
  for j = 0; j < self.ncols; j = j + 1 {
    let xj = x[j]
    if xj != 0.0 {
      for k = self.col_ptr[j]; k < self.col_ptr[j + 1]; k = k + 1 {
        let i = self.row_idx[k]
        y[i] = y[i] + self.values[k] * xj
      }
    }
  }
  y
}

///|
/// Insertion sort of `(row, value)` pairs by row index.
///
/// Insertion sort is intentional: columns of an LP are short, and this keeps the
/// reference path allocation free and obviously correct.
fn sort_entries(entries : Array[(Int32, Double)]) -> Array[(Int32, Double)] {
  let out = entries.copy()
  let n = out.length()
  for i = 1; i < n; i = i + 1 {
    let cur = out[i]
    let mut j = i - 1
    while j >= 0 && out[j].0 > cur.0 {
      out[j + 1] = out[j]
      j = j - 1
    }
    out[j + 1] = cur
  }
  out
}