///|
/// Symbolic matrix expression nodes and front doors.
///
/// Current Limits:
/// - This layer models symbolic matrix syntax but only evaluates operations implemented in `symmatrices`.
/// - Unsupported symbolic rewrites stay unevaluated instead of introducing a separate assumption system.
///
/// - Does: Represents symbolic matrix syntax such as symbols, products, sums, block layouts, transpose, inverse, and integer powers.
/// - Input: Each variant stores matrix metadata or child `MatrixExpr` nodes.
/// - Returns: One `MatrixExpr` tree.
/// - Limits: This is a syntax layer, not a full theorem prover; semantically equivalent matrix expressions can remain structurally different until a caller evaluates or collapses them.
pub enum MatrixExpr {
  Concrete(Matrix)
  Symbol(String, Int, Int)
  Identity(Int)
  Zero(Int, Int)
  Add(Array[MatrixExpr])
  Mul(Array[MatrixExpr])
  BlockDiag(Array[MatrixExpr])
  Block(Array[Array[MatrixExpr]])
  Transpose(MatrixExpr)
  Inverse(MatrixExpr)
  Pow(MatrixExpr, Int)
}

///|
fn reverse_matrix_exprs(items : Array[MatrixExpr]) -> Array[MatrixExpr] {
  let out : Array[MatrixExpr] = []
  for i in 0.. Array[Array[MatrixExpr]] {
  if rows.is_empty() {
    return []
  }
  let cols = rows[0].length()
  let out : Array[Array[MatrixExpr]] = []
  for j in 0.. Bool {
  if lhs.length() != rhs.length() {
    return false
  }
  for i in 0.. return false }
      let rhs_shape = rhs[i][j].shape() catch { _ => return false }
      if lhs_shape != rhs_shape {
        return false
      }
    }
  }
  true
}

///|
fn block_mul_compatible(
  lhs : Array[Array[MatrixExpr]],
  rhs : Array[Array[MatrixExpr]],
) -> Bool {
  if lhs.is_empty() || rhs.is_empty() || lhs[0].length() != rhs.length() {
    return false
  }
  for k in 0.. return false }
    let right_shape = rhs[k][0].shape() catch { _ => return false }
    if left_shape.1 != right_shape.0 {
      return false
    }
  }
  true
}

///|
fn block_add_rows(
  lhs : Array[Array[MatrixExpr]],
  rhs : Array[Array[MatrixExpr]],
) -> Array[Array[MatrixExpr]] {
  let out : Array[Array[MatrixExpr]] = []
  for i in 0.. Array[Array[MatrixExpr]] {
  let out : Array[Array[MatrixExpr]] = []
  for i in 0.. Array[MatrixExpr] {
  let out : Array[MatrixExpr] = []
  for item in items {
    match item {
      MatrixExpr::Add(inner) =>
        for nested in normalize_add_items(inner) {
          out.push(nested)
        }
      MatrixExpr::Zero(_, _) => ()
      _ => out.push(item)
    }
  }
  out
}

///|
fn normalize_mul_items(items : Array[MatrixExpr]) -> Array[MatrixExpr] {
  let out : Array[MatrixExpr] = []
  for item in items {
    match item {
      MatrixExpr::Mul(inner) =>
        for nested in normalize_mul_items(inner) {
          out.push(nested)
        }
      MatrixExpr::Identity(_) => ()
      _ => out.push(item)
    }
  }
  out
}

///|
fn normalize_block_diag_items(items : Array[MatrixExpr]) -> Array[MatrixExpr] {
  let out : Array[MatrixExpr] = []
  for item in items {
    match item {
      MatrixExpr::BlockDiag(inner) =>
        for nested in normalize_block_diag_items(inner) {
          out.push(nested)
        }
      _ => out.push(item)
    }
  }
  out
}

///|
fn all_concrete(items : Array[MatrixExpr]) -> Bool {
  items.all(fn(item) {
    match item {
      MatrixExpr::Concrete(_) => true
      _ => false
    }
  })
}

///|
fn concrete_add(items : Array[MatrixExpr]) -> MatrixExpr {
  let mut acc = match items[0] {
    MatrixExpr::Concrete(m) => m
    _ => abort("unreachable")
  }
  for item in items[1:] {
    match item {
      MatrixExpr::Concrete(m) => acc = acc + m
      _ => abort("unreachable")
    }
  }
  MatrixExpr::Concrete(acc)
}

///|
fn concrete_mul(items : Array[MatrixExpr]) -> MatrixExpr {
  let mut acc = match items[0] {
    MatrixExpr::Concrete(m) => m
    _ => abort("unreachable")
  }
  for item in items[1:] {
    match item {
      MatrixExpr::Concrete(m) => acc = acc * m
      _ => abort("unreachable")
    }
  }
  MatrixExpr::Concrete(acc)
}

///|
fn matrix_expr_same(lhs : MatrixExpr, rhs : MatrixExpr) -> Bool {
  (try? lhs.shape()) == (try? rhs.shape()) && lhs.to_string() == rhs.to_string()
}

///|
fn matrix_expr_power_view(expr : MatrixExpr) -> (MatrixExpr, Int) {
  match expr {
    MatrixExpr::Pow(base, exp) => (base, exp)
    _ => (expr, 1)
  }
}

///|
fn matrix_expr_cancel_or_merge_pair(
  lhs : MatrixExpr,
  rhs : MatrixExpr,
) -> MatrixExpr? {
  match (lhs, rhs) {
    (MatrixExpr::Inverse(inner), other) =>
      if matrix_expr_same(inner, other) {
        let shape = try? inner.shape()
        match shape {
          Ok((rows, cols)) =>
            if rows == cols {
              Some(identity_expr(rows))
            } else {
              None
            }
          Err(_) => None
        }
      } else {
        None
      }
    (other, MatrixExpr::Inverse(inner)) =>
      if matrix_expr_same(other, inner) {
        let shape = try? inner.shape()
        match shape {
          Ok((rows, cols)) =>
            if rows == cols {
              Some(identity_expr(rows))
            } else {
              None
            }
          Err(_) => None
        }
      } else {
        None
      }
    _ => {
      let left = matrix_expr_power_view(lhs)
      let right = matrix_expr_power_view(rhs)
      if matrix_expr_same(left.0, right.0) {
        Some(matrix_expr_pow(left.0, left.1 + right.1))
      } else {
        None
      }
    }
  }
}

///|
fn blockdiag_identity_like(items : Array[MatrixExpr]) -> MatrixExpr? {
  let identities : Array[MatrixExpr] = []
  for item in items {
    let shape = try? item.shape()
    match shape {
      Ok((rows, cols)) => {
        if rows != cols {
          return None
        }
        identities.push(identity_expr(rows))
      }
      Err(_) => return None
    }
  }
  Some(block_diag_matrix_expr(identities))
}

///|
fn block_identity_like(rows : Array[Array[MatrixExpr]]) -> MatrixExpr? {
  if rows.is_empty() {
    return Some(zero_matrix_expr(0, 0))
  }
  if rows.length() != rows[0].length() {
    return None
  }
  let heights : Array[Int] = []
  let widths : Array[Int] = []
  for i in 0..
          if height == 0 {
            height = rows0
          } else if height != rows0 {
            return None
          }
        Err(_) => return None
      }
    }
    heights.push(height)
  }
  for j in 0..
          if width == 0 {
            width = cols0
          } else if width != cols0 {
            return None
          }
        Err(_) => return None
      }
    }
    widths.push(width)
  }
  for i in 0.. MatrixExpr {
  MatrixExpr::Symbol(name, rows, cols)
}

///|
/// Lift a concrete dense matrix into the symbolic matrix expression layer.
///
/// - Does: Wraps a `Matrix` as `MatrixExpr::Concrete`.
/// - Input: A dense `Matrix`.
/// - Returns: A `MatrixExpr`.
/// - Limits: This is a structural wrapper only; it does not copy or simplify entries beyond the underlying value semantics.
pub fn matrix_expr(m : Matrix) -> MatrixExpr {
  MatrixExpr::Concrete(m)
}

///|
/// Create a symbolic identity matrix expression.
///
/// - Does: Builds an identity-node placeholder.
/// - Input: A single size `n`.
/// - Returns: `MatrixExpr::Identity(n)`.
/// - Limits: This front door does not validate that `n` is non-negative.
pub fn identity_expr(n : Int) -> MatrixExpr {
  MatrixExpr::Identity(n)
}

///|
/// Create a symbolic zero matrix expression.
///
/// - Does: Builds a zero-matrix placeholder with explicit shape.
/// - Input: Row and column counts.
/// - Returns: `MatrixExpr::Zero(rows, cols)`.
/// - Limits: This front door does not validate that the provided dimensions are non-negative.
pub fn zero_matrix_expr(rows : Int, cols : Int) -> MatrixExpr {
  MatrixExpr::Zero(rows, cols)
}

///|
/// Build a symbolic matrix sum.
///
/// - Does: Flattens nested adds, drops symbolic zero terms, and eagerly combines concrete and block-diagonal operands when possible.
/// - Input: An `Array[MatrixExpr]`.
/// - Returns: A simplified `MatrixExpr`.
/// - Limits: If shapes are inconsistent, the expression may stay unevaluated until callers ask for `shape()` or `eval()`.
pub fn matrix_expr_add(items : Array[MatrixExpr]) -> MatrixExpr {
  if items.is_empty() {
    return MatrixExpr::Zero(0, 0)
  }
  let normalized = normalize_add_items(items)
  if normalized.length() > 1 && all_concrete(normalized) {
    return concrete_add(normalized)
  }
  if normalized.length() > 1 &&
    normalized.all(fn(item) {
      match item {
        MatrixExpr::BlockDiag(_) => true
        _ => false
      }
    }) {
    let first = match normalized[0] {
      MatrixExpr::BlockDiag(items) => items
      _ => []
    }
    if normalized.all(fn(item) {
        match item {
          MatrixExpr::BlockDiag(parts) => parts.length() == first.length()
          _ => false
        }
      }) {
      let blocks : Array[MatrixExpr] = []
      for i in 0.. pieces.push(parts[i])
            _ => ()
          }
        }
        blocks.push(matrix_expr_add(pieces))
      }
      return matrix_expr_block_diag(blocks)
    }
  }
  match normalized.length() {
    0 =>
      match items[0] {
        MatrixExpr::Zero(rows, cols) => MatrixExpr::Zero(rows, cols)
        _ => items[0]
      }
    1 => normalized[0]
    _ => MatrixExpr::Add(normalized)
  }
}

///|
/// Build a symbolic matrix product.
///
/// - Does: Flattens nested products, removes identities, folds adjacent concrete factors, and performs a small amount of inverse/power cancellation.
/// - Input: An `Array[MatrixExpr]`.
/// - Returns: A simplified `MatrixExpr`.
/// - Limits: Incompatible dimensions are not rejected here unless they are needed for a zero or identity shortcut; callers should use `shape()` or `eval()` for strict checking.
pub fn matrix_expr_mul(items : Array[MatrixExpr]) -> MatrixExpr {
  if items.is_empty() {
    return MatrixExpr::Identity(0)
  }
  let normalized = normalize_mul_items(items)
  if normalized.length() > 1 && all_concrete(normalized) {
    return concrete_mul(normalized)
  }
  if normalized.any(fn(item) {
      match item {
        MatrixExpr::Zero(_, _) => true
        _ => false
      }
    }) {
    let shape = try? MatrixExpr::Mul(normalized).shape()
    match shape {
      Ok((rows, cols)) => return MatrixExpr::Zero(rows, cols)
      Err(_) => ()
    }
  }
  let reduced : Array[MatrixExpr] = []
  for item in normalized {
    if reduced.is_empty() {
      reduced.push(item)
      continue
    }
    let last = reduced[reduced.length() - 1]
    match matrix_expr_cancel_or_merge_pair(last, item) {
      Some(merged) => {
        ignore(reduced.pop())
        match merged {
          MatrixExpr::Identity(_) => ()
          _ => reduced.push(merged)
        }
      }
      None =>
        match (last, item) {
          (MatrixExpr::Concrete(lhs), MatrixExpr::Concrete(rhs)) => {
            ignore(reduced.pop())
            reduced.push(MatrixExpr::Concrete(lhs * rhs))
          }
          _ => reduced.push(item)
        }
    }
  }
  if reduced.length() > 1 &&
    reduced.all(fn(item) {
      match item {
        MatrixExpr::BlockDiag(_) => true
        _ => false
      }
    }) {
    let first = match reduced[0] {
      MatrixExpr::BlockDiag(items) => items
      _ => []
    }
    if reduced.all(fn(item) {
        match item {
          MatrixExpr::BlockDiag(parts) => parts.length() == first.length()
          _ => false
        }
      }) {
      let blocks : Array[MatrixExpr] = []
      for i in 0.. pieces.push(parts[i])
            _ => ()
          }
        }
        blocks.push(matrix_expr_mul(pieces))
      }
      return matrix_expr_block_diag(blocks)
    }
  }
  match reduced.length() {
    0 =>
      try MatrixExpr::Mul(items).shape() catch {
        _ =>
          match items[0] {
            MatrixExpr::Identity(n) => MatrixExpr::Identity(n)
            _ => items[0]
          }
      } noraise {
        (rows, cols) =>
          if rows == cols {
            MatrixExpr::Identity(rows)
          } else {
            items[0]
          }
      }
    1 => reduced[0]
    _ => MatrixExpr::Mul(reduced)
  }
}

///|
/// Build a symbolic block-diagonal expression.
///
/// - Does: Flattens nested block-diagonal nodes and preserves the remaining pieces structurally.
/// - Input: An `Array[MatrixExpr]`.
/// - Returns: A simplified `MatrixExpr`.
/// - Limits: Empty input returns a `0 x 0` zero expression.
pub fn matrix_expr_block_diag(items : Array[MatrixExpr]) -> MatrixExpr {
  if items.is_empty() {
    return MatrixExpr::Zero(0, 0)
  }
  let normalized = normalize_block_diag_items(items)
  match normalized.length() {
    0 =>
      match items[0] {
        MatrixExpr::Zero(rows, cols) => MatrixExpr::Zero(rows, cols)
        _ => items[0]
      }
    1 => normalized[0]
    _ => MatrixExpr::BlockDiag(normalized)
  }
}

///|
/// Build a symbolic block matrix expression.
///
/// - Does: Stores the block layout directly, or unwraps a `1 x 1` block grid back into its only entry.
/// - Input: Rows of block expressions.
/// - Returns: A `MatrixExpr`.
/// - Limits: Rectangularity and block-shape compatibility are validated later by `shape()` or `eval()`.
pub fn matrix_expr_block(rows : Array[Array[MatrixExpr]]) -> MatrixExpr {
  if rows.length() == 1 && rows[0].length() == 1 {
    return rows[0][0]
  }
  MatrixExpr::Block(rows)
}

///|
/// Convenience alias for `matrix_expr_add`.
///
/// - Does: Delegates to `matrix_expr_add`.
/// - Input: An `Array[MatrixExpr]`.
/// - Returns: A simplified `MatrixExpr`.
/// - Limits: Incompatible shapes are still detected later by `shape()` or `eval()`.
pub fn mat_add(items : Array[MatrixExpr]) -> MatrixExpr {
  matrix_expr_add(items)
}

///|
/// Convenience alias for `matrix_expr_mul`.
///
/// - Does: Delegates to `matrix_expr_mul`.
/// - Input: An `Array[MatrixExpr]`.
/// - Returns: A simplified `MatrixExpr`.
/// - Limits: Incompatible dimensions are still detected later by `shape()` or `eval()`.
pub fn mat_mul(items : Array[MatrixExpr]) -> MatrixExpr {
  matrix_expr_mul(items)
}

///|
/// Convenience alias for `matrix_expr_block`.
///
/// - Does: Delegates to `matrix_expr_block`.
/// - Input: Rows of block expressions.
/// - Returns: A `MatrixExpr`.
/// - Limits: Block-layout validation is still deferred to `shape()` or `eval()`.
pub fn block_matrix_expr(rows : Array[Array[MatrixExpr]]) -> MatrixExpr {
  matrix_expr_block(rows)
}

///|
/// Convenience alias for `matrix_expr_block_diag`.
///
/// - Does: Delegates to `matrix_expr_block_diag`.
/// - Input: An `Array[MatrixExpr]`.
/// - Returns: A `MatrixExpr`.
/// - Limits: Empty input still becomes a `0 x 0` zero expression.
pub fn block_diag_matrix_expr(items : Array[MatrixExpr]) -> MatrixExpr {
  matrix_expr_block_diag(items)
}

///|
/// Transpose a symbolic matrix expression.
///
/// - Does: Pushes transpose through supported nodes and simplifies double-transpose cases.
/// - Input: Any `MatrixExpr`.
/// - Returns: A `MatrixExpr`.
/// - Limits: Unsupported rewrites remain as explicit `Transpose(...)` nodes.
pub fn matrix_expr_transpose(expr : MatrixExpr) -> MatrixExpr {
  match expr {
    MatrixExpr::Transpose(inner) => inner
    MatrixExpr::Identity(n) => MatrixExpr::Identity(n)
    MatrixExpr::Zero(rows, cols) => MatrixExpr::Zero(cols, rows)
    MatrixExpr::Add(items) =>
      mat_add(items.map(fn(item) { matrix_expr_transpose(item) }))
    MatrixExpr::Mul(items) =>
      mat_mul(
        reverse_matrix_exprs(items).map(fn(item) { matrix_expr_transpose(item) }),
      )
    MatrixExpr::BlockDiag(items) =>
      block_diag_matrix_expr(
        items.map(fn(item) { matrix_expr_transpose(item) }),
      )
    MatrixExpr::Block(rows) => block_matrix_expr(transpose_block_rows(rows))
    _ => MatrixExpr::Transpose(expr)
  }
}

///|
/// Invert a symbolic matrix expression.
///
/// - Does: Simplifies inverse-on-inverse, inverse-of-transpose, inverse-of-products, and a small set of block patterns.
/// - Input: Any `MatrixExpr`.
/// - Returns: A `MatrixExpr`.
/// - Limits: General block-matrix inversion is not expanded here; unsupported cases remain as `Inverse(...)`.
pub fn matrix_expr_inverse(expr : MatrixExpr) -> MatrixExpr {
  match expr {
    MatrixExpr::Identity(n) => MatrixExpr::Identity(n)
    MatrixExpr::Inverse(inner) => inner
    MatrixExpr::Transpose(inner) =>
      matrix_expr_transpose(matrix_expr_inverse(inner))
    MatrixExpr::BlockDiag(items) =>
      block_diag_matrix_expr(items.map(fn(item) { matrix_expr_inverse(item) }))
    MatrixExpr::Block(rows) =>
      if rows.length() == 1 && rows[0].length() == 1 {
        matrix_expr_inverse(rows[0][0])
      } else {
        MatrixExpr::Inverse(expr)
      }
    MatrixExpr::Mul(items) =>
      mat_mul(
        reverse_matrix_exprs(items).map(fn(item) { matrix_expr_inverse(item) }),
      )
    MatrixExpr::Pow(inner, exp) => matrix_expr_pow(inner, -exp)
    _ => MatrixExpr::Inverse(expr)
  }
}

///|
/// Raise a symbolic matrix expression to an integer power.
///
/// - Does: Simplifies identity, zero, transpose, inverse, and nested-power cases for integer exponents.
/// - Input: A `MatrixExpr` and an integer exponent.
/// - Returns: A `MatrixExpr`.
/// - Limits: Non-square inputs are only rejected when the zero-power shortcut needs the shape; other invalid powers can stay unevaluated until `eval()`.
pub fn matrix_expr_pow(expr : MatrixExpr, exp : Int) -> MatrixExpr {
  if exp == 0 {
    let shape = try? expr.shape()
    match shape {
      Ok((rows, cols)) => if rows == cols { return identity_expr(rows) }
      Err(_) => ()
    }
  }
  if exp == 1 {
    return expr
  }
  if exp == -1 {
    return matrix_expr_inverse(expr)
  }
  match expr {
    MatrixExpr::Identity(n) => MatrixExpr::Identity(n)
    MatrixExpr::Zero(rows, cols) if rows == cols && exp > 0 =>
      MatrixExpr::Zero(rows, cols)
    MatrixExpr::Transpose(inner) =>
      matrix_expr_transpose(matrix_expr_pow(inner, exp))
    MatrixExpr::Pow(inner, inner_exp) => MatrixExpr::Pow(inner, inner_exp * exp)
    MatrixExpr::BlockDiag(items) if exp > 0 =>
      block_diag_matrix_expr(items.map(fn(item) { matrix_expr_pow(item, exp) }))
    _ => MatrixExpr::Pow(expr, exp)
  }
}

///|
/// Compute the shape of a symbolic matrix expression.
///
/// - Does: Checks the structural dimensions implied by the expression tree.
/// - Input: Any `MatrixExpr`.
/// - Returns: A pair `(rows, cols)`.
/// - Limits: Raises `MatrixError::ShapeError` when add/mul/block layouts are inconsistent.
pub fn MatrixExpr::shape(self : MatrixExpr) -> (Int, Int) raise MatrixError {
  match self {
    MatrixExpr::Concrete(m) => m.shape()
    MatrixExpr::Symbol(_, rows, cols) => (rows, cols)
    MatrixExpr::Identity(n) => (n, n)
    MatrixExpr::Zero(rows, cols) => (rows, cols)
    MatrixExpr::Add(items) => {
      guard !items.is_empty() else {
        raise MatrixError::ShapeError("empty MatrixExpr Add")
      }
      let shape = items[0].shape()
      for item in items[1:] {
        if item.shape() != shape {
          raise MatrixError::ShapeError(
            "matrix expression add expects equal shapes",
          )
        }
      }
      shape
    }
    MatrixExpr::Mul(items) => {
      guard !items.is_empty() else {
        raise MatrixError::ShapeError("empty MatrixExpr Mul")
      }
      let first = items[0].shape()
      let mut last = first
      for item in items[1:] {
        let current = item.shape()
        if last.1 != current.0 {
          raise MatrixError::ShapeError(
            "matrix expression mul expects aligned dimensions",
          )
        }
        last = current
      }
      (first.0, last.1)
    }
    MatrixExpr::BlockDiag(items) => {
      let mut rows = 0
      let mut cols = 0
      for item in items {
        let shape = item.shape()
        rows += shape.0
        cols += shape.1
      }
      (rows, cols)
    }
    MatrixExpr::Block(rows0) => {
      guard !rows0.is_empty() else { return (0, 0) }
      let mut rows = 0
      let col_widths : Array[Int] = []
      let first_width = rows0[0].length()
      for row in rows0 {
        if row.length() != first_width {
          raise MatrixError::ShapeError(
            "block matrix expression expects rectangular block layout",
          )
        }
      }
      for j in 0.. {
      let shape = inner.shape()
      (shape.1, shape.0)
    }
    MatrixExpr::Inverse(inner) | MatrixExpr::Pow(inner, _) => inner.shape()
  }
}

///|
/// Render a symbolic matrix expression as text.
///
/// - Does: Produces a stable textual representation of the expression tree.
/// - Input: One `MatrixExpr`.
/// - Returns: A `String`.
/// - Limits: The text reflects structural form rather than full algebraic equivalence.
pub impl Show for MatrixExpr with to_string(self) {
  match self {
    MatrixExpr::Concrete(m) => m.to_string()
    MatrixExpr::Symbol(name, _, _) => name
    MatrixExpr::Identity(n) => "Identity(" + n.to_string() + ")"
    MatrixExpr::Zero(rows, cols) =>
      "ZeroMatrix(" + rows.to_string() + ", " + cols.to_string() + ")"
    MatrixExpr::Add(items) =>
      "(" + items.map(fn(item) { item.to_string() }).join(" + ") + ")"
    MatrixExpr::Mul(items) =>
      "(" + items.map(fn(item) { item.to_string() }).join(" * ") + ")"
    MatrixExpr::BlockDiag(items) =>
      "BlockDiag(" + items.map(fn(item) { item.to_string() }).join(", ") + ")"
    MatrixExpr::Block(rows) =>
      "Block(" +
      rows
      .map(fn(row) {
        "[" + row.map(fn(item) { item.to_string() }).join(", ") + "]"
      })
      .join(", ") +
      ")"
    MatrixExpr::Transpose(inner) => inner.to_string() + ".T"
    MatrixExpr::Inverse(inner) => inner.to_string() + "^-1"
    MatrixExpr::Pow(inner, exp) => inner.to_string() + "^" + exp.to_string()
  }
}

///|
/// Write the textual form of a symbolic matrix expression into a logger.
///
/// - Does: Delegates to `to_string()` and writes the result.
/// - Input: One `MatrixExpr` and a `Logger`.
/// - Returns: `Unit`.
/// - Limits: Output stability matters for tests and diagnostics.
pub impl Show for MatrixExpr with output(self, logger : &Logger) -> Unit {
  logger.write_string(self.to_string())
}

///|
/// Evaluate a symbolic matrix expression against a binding environment.
///
/// - Does: Replaces symbolic matrices from `env` and executes supported matrix operations to produce a concrete dense matrix.
/// - Input: A `MatrixExpr` and `Map[String, Matrix]` environment for symbol bindings.
/// - Returns: A dense `Matrix`.
/// - Limits: Raises `MatrixError::ValueError` for missing symbol bindings and shape-related matrix errors when an operation cannot be evaluated.
pub fn MatrixExpr::eval(
  self : MatrixExpr,
  env : Map[String, Matrix],
) -> Matrix raise MatrixError {
  match self {
    MatrixExpr::Concrete(m) => m
    MatrixExpr::Symbol(name, rows, cols) =>
      match env.get(name) {
        Some(value) => {
          if value.shape() != (rows, cols) {
            raise MatrixError::ShapeError(
              "matrix symbol binding has wrong shape",
            )
          }
          value
        }
        None => raise MatrixError::ValueError("missing matrix symbol binding")
      }
    MatrixExpr::Identity(n) => eye(n)
    MatrixExpr::Zero(rows, cols) => zeros(rows, cols=Some(cols))
    MatrixExpr::Add(items) => {
      guard !items.is_empty() else { return zeros(0, cols=Some(0)) }
      let mut acc = items[0].eval(env)
      for item in items[1:] {
        let next = item.eval(env)
        if acc.shape() != next.shape() {
          raise MatrixError::ShapeError(
            "matrix expression add expects equal shapes",
          )
        }
        acc = acc + next
      }
      acc
    }
    MatrixExpr::Mul(items) => {
      guard !items.is_empty() else { return zeros(0, cols=Some(0)) }
      let mut acc = items[0].eval(env)
      for item in items[1:] {
        let next = item.eval(env)
        if acc.cols != next.rows {
          raise MatrixError::ShapeError(
            "matrix expression mul expects aligned dimensions",
          )
        }
        acc = acc * next
      }
      acc
    }
    MatrixExpr::BlockDiag(items) => {
      let mats : Array[Matrix] = []
      for item in items {
        mats.push(item.eval(env))
      }
      block_diag(mats)
    }
    MatrixExpr::Block(rows) => {
      let mats : Array[Array[Matrix]] = []
      for row in rows {
        let mat_row : Array[Matrix] = []
        for item in row {
          mat_row.push(item.eval(env))
        }
        mats.push(mat_row)
      }
      block_matrix(mats)
    }
    MatrixExpr::Transpose(inner) => inner.eval(env).transpose()
    MatrixExpr::Inverse(inner) => inner.eval(env).inv()
    MatrixExpr::Pow(inner, exp) => {
      let base = inner.eval(env)
      if !base.is_square() {
        raise MatrixError::NonSquareMatrixError(
          "matrix power requires a square matrix",
        )
      }
      if exp == 0 {
        return eye(base.rows)
      }
      if exp < 0 {
        let inv = base.inv()
        let mut out = inv
        for _ in 1..<-exp {
          out = out * inv
        }
        return out
      }
      let mut out = base
      for _ in 1.. Matrix raise MatrixError {
  self.eval(env)
}

///|
/// Convert a symbolic matrix expression into an explicit dense matrix.
///
/// - Does: Alias of `eval`, intended for callers that want a fully explicit matrix result.
/// - Input: A `MatrixExpr` and symbol environment.
/// - Returns: A dense `Matrix`.
/// - Limits: Propagates the same evaluation errors as `eval`.
pub fn MatrixExpr::as_explicit(
  self : MatrixExpr,
  env : Map[String, Matrix],
) -> Matrix raise MatrixError {
  self.eval(env)
}

///|
/// Collapse block-oriented symbolic wrappers where local rewrites are available.
///
/// - Does: Recursively simplifies nested block, block-diagonal, identity, transpose, inverse, and product forms.
/// - Input: Any `MatrixExpr`.
/// - Returns: A simplified `MatrixExpr`.
/// - Limits: This is not a full symbolic canonicalizer; many semantically equivalent expressions remain structurally different.
pub fn MatrixExpr::block_collapse(self : MatrixExpr) -> MatrixExpr {
  match self {
    MatrixExpr::Add(items) => {
      let collapsed = items.map(fn(item) { item.block_collapse() })
      if collapsed.any(fn(item) {
          match item {
            MatrixExpr::BlockDiag(_) => true
            _ => false
          }
        }) {
        let normalized : Array[MatrixExpr] = []
        let mut changed = false
        let mut template : MatrixExpr? = None
        for item in collapsed {
          match item {
            MatrixExpr::BlockDiag(_) => {
              template = Some(item)
              break
            }
            _ => ()
          }
        }
        match template {
          Some(MatrixExpr::BlockDiag(parts)) =>
            for item in collapsed {
              match item {
                MatrixExpr::Identity(_) =>
                  match blockdiag_identity_like(parts) {
                    Some(identity) => {
                      normalized.push(identity)
                      changed = true
                    }
                    None => normalized.push(item)
                  }
                _ => normalized.push(item)
              }
            }
          _ => ()
        }
        if changed {
          return mat_add(normalized).block_collapse()
        }
      }
      if collapsed.all(fn(item) {
          match item {
            MatrixExpr::BlockDiag(_) => true
            _ => false
          }
        }) {
        let first = match collapsed[0] {
          MatrixExpr::BlockDiag(items) => items
          _ => []
        }
        let aligned = collapsed.all(fn(item) {
          match item {
            MatrixExpr::BlockDiag(parts) => parts.length() == first.length()
            _ => false
          }
        })
        if aligned {
          let blocks : Array[MatrixExpr] = []
          for i in 0.. pieces.push(parts[i])
                _ => ()
              }
            }
            blocks.push(mat_add(pieces))
          }
          return block_diag_matrix_expr(blocks)
        }
      }
      if collapsed.all(fn(item) {
          match item {
            MatrixExpr::Block(_) => true
            _ => false
          }
        }) {
        let rows = match collapsed[0] {
          MatrixExpr::Block(rows) => rows
          _ => []
        }
        let aligned = collapsed.all(fn(item) {
          match item {
            MatrixExpr::Block(other) => block_same_layout(rows, other)
            _ => false
          }
        })
        if aligned {
          let mut acc = rows
          for item in collapsed[1:] {
            match item {
              MatrixExpr::Block(other) => acc = block_add_rows(acc, other)
              _ => ()
            }
          }
          return block_matrix_expr(acc)
        }
      }
      if collapsed.any(fn(item) {
          match item {
            MatrixExpr::Block(_) => true
            _ => false
          }
        }) {
        let normalized : Array[MatrixExpr] = []
        let mut changed = false
        let mut template : MatrixExpr? = None
        for item in collapsed {
          match item {
            MatrixExpr::Block(_) => {
              template = Some(item)
              break
            }
            _ => ()
          }
        }
        match template {
          Some(MatrixExpr::Block(rows)) =>
            for item in collapsed {
              match item {
                MatrixExpr::Identity(_) =>
                  match block_identity_like(rows) {
                    Some(identity) => {
                      normalized.push(identity)
                      changed = true
                    }
                    None => normalized.push(item)
                  }
                _ => normalized.push(item)
              }
            }
          _ => ()
        }
        if changed {
          return mat_add(normalized).block_collapse()
        }
      }
      matrix_expr_add(collapsed)
    }
    MatrixExpr::Mul(items) => {
      let collapsed = items.map(fn(item) { item.block_collapse() })
      guard !collapsed.is_empty() else { return MatrixExpr::Identity(0) }
      let mut acc = collapsed[0]
      for item in collapsed[1:] {
        acc = match (acc, item) {
          (MatrixExpr::BlockDiag(lhs), MatrixExpr::BlockDiag(rhs)) =>
            if lhs.length() == rhs.length() {
              block_diag_matrix_expr(
                Array::makei(lhs.length(), fn(i) { mat_mul([lhs[i], rhs[i]]) }),
              )
            } else {
              mat_mul([MatrixExpr::BlockDiag(lhs), MatrixExpr::BlockDiag(rhs)])
            }
          (MatrixExpr::Block(lhs), MatrixExpr::Block(rhs)) =>
            if block_mul_compatible(lhs, rhs) {
              block_matrix_expr(block_mul_rows(lhs, rhs))
            } else {
              mat_mul([MatrixExpr::Block(lhs), MatrixExpr::Block(rhs)])
            }
          (MatrixExpr::BlockDiag(lhs), MatrixExpr::Block(rhs)) =>
            if lhs.length() == rhs.length() {
              let rows : Array[Array[MatrixExpr]] = []
              for i in 0..
            if !lhs.is_empty() && lhs[0].length() == rhs.length() {
              let rows : Array[Array[MatrixExpr]] = []
              for i in 0.. mat_mul([left, right])
        }
      }
      acc
    }
    MatrixExpr::BlockDiag(items) =>
      matrix_expr_block_diag(items.map(fn(item) { item.block_collapse() }))
    MatrixExpr::Block(rows) => {
      let collapsed = rows.map(fn(row) {
        row.map(fn(item) { item.block_collapse() })
      })
      matrix_expr_block(collapsed)
    }
    MatrixExpr::Transpose(inner) =>
      matrix_expr_transpose(inner.block_collapse())
    MatrixExpr::Inverse(inner) => matrix_expr_inverse(inner.block_collapse())
    MatrixExpr::Pow(inner, exp) => matrix_expr_pow(inner.block_collapse(), exp)
    _ => self
  }
}