-- Copyright 1993-1998, by the Cecil Project
-- Department of Computer Science and Engineering, University of Washington
-- See the LICENSE file for license information.

(--DOC
  Matrices support querying the number of rows and columns of a
  matrix, extracting an element, a row, or a column of the matrix, and
  converting the matrix into a vector of row vectors. The `do' and
  `indices_do' methods support iterating over the matrix. Conformable
  matrices can be added and multiplied. A matrix of comparable values
  is itself comparable, pointwise. Mutable matrices support changing a
  given matrix element.

  One concrete implementation of matrices exists, based on
  representing matrices by a vector of vectors. The
  `new_vector_matrix_init' functions supports creating a new
  `vector_matrix' of a given size with its elements initialized as
  specified by the `init_fn' closure.
--)


abstract object matrix[T];
  extend type matrix[`T] subtypes matrix[`S >= T];

  -- query dimensions of matrix
  signature num_rows(matrix[`T]):int;
  signature num_cols(matrix[`T]):int;

  -- fetch element of matrix
  signature fetch(matrix[`T], row:int, col:int):T;

  -- fetch row and/or column slice(s) of matrix
  method row(m@:matrix[`T], row:int):indexed[T] {
      new_i_vector_init[T](m.num_cols, &(col:int){ m.fetch(row, col) }) }
  method col(m@:matrix[`T], col:int):indexed[T] {
      new_i_vector_init[T](m.num_rows, &(row:int){ m.fetch(row, col) }) }
  method rows(m@:matrix[`T]):indexed[indexed[T]] {
      new_i_vector_init[indexed[T]](m.num_rows, &(row:int){ m.row(row) }) }

  -- iterators
  method indices_do(m@:matrix[`T], c:&(int,int):void):void (** inline **) {
      m.num_rows.do(&(row:int){
	  m.num_cols.do(&(col:int){
	      eval(c, row, col);
	  });
      });
  }
  method do(m@:matrix[`T], c:&(int,int,T):void):void (** inline **) {
      m.indices_do(&(row:int,col:int){
	  eval(c, row, col, m.fetch(row,col));
      });
  }

  -- add conformant matrices
  method +(m1@:matrix[`T <= num], m2@:matrix[T]):matrix[T] {
      assert(m1.num_rows = m2.num_rows & { m1.num_cols = m2.num_cols },
	     "cannot add matrices that aren't conformant");
      m1.copy_init(m1.num_rows, m1.num_cols, &(row:int, col:int){
	  m1.fetch(row,col) + m2.fetch(row,col) }) }

  -- multiply conformant matrices
  method *(m1@:matrix[`T <= num], m2@:matrix[T]):matrix[T|int] {
      assert(m1.num_cols = m2.num_rows,
	     "cannot multiply matrices that aren't conformant");
      m1.copy_init[T|int](m1.num_rows, m2.num_cols, &(row:int, col:int){
	  (-- really should write this as a dot-product of two vectors:
	      m1.row(row) * m2.col(col)
	   --)
	  let var sum:T|int := 0;
	  m1.num_cols.do(&(i:int){
	      sum := sum + m1.fetch(row, i) * m2.fetch(i, col);
	  });
	  sum }) }

  -- multiply conformant matrices, using an ugly triply-nested-loop-based alg.
  method *_ugly(m1@:matrix[`T <= num], m2@:matrix[T])
     	       :matrix[T|int] {
      assert(m1.num_cols = m2.num_rows,
	     "cannot multiply matrices that aren't conformant");
      let result:m_matrix[T|int] :=
        m1.copy_mutable_init[T|int](m1.num_rows, m2.num_cols,
					&(row:int, col:int){ 0 });
      m1.num_rows.do(&(row:int){
	  m2.num_cols.do(&(col:int){
	      let var sum:T|int := 0;
	      m1.num_cols.do(&(i:int){
		  sum := sum + m1.fetch(row, i) * m2.fetch(i, col);
	      });
	      result.store(row, col, sum);
	  });
      });
      result }

  extend matrix[`T <= comparable[T]] isa comparable[matrix[T]];
  method =(m1@:matrix[`T <= comparable[T]], m2@:matrix[`T]):bool {
      m1 == m2 | { m1.rows = m2.rows } }

  method print_string(m@:matrix[`T]):string {
      let var s:string := "matrix{";
      let var first:bool := true;
      m.num_rows.do(&(row:int){
	  if(first, { first := false; }, { s := s || ",\n       "; });
	  s := s || "{" || m.row(row).elems_print_string || "}";
      });
      s || "}" }

  method copy_init(m@:matrix[`T], num_rows:int, num_cols:int,
		   init:&(int,int):T):matrix[T] {
      copy_init[T](m, num_rows, num_cols, init) }
  signature copy_init[T](matrix[T], num_rows:int, num_cols:int,
			 init:&(int,int):T):matrix[T];

  method copy_mutable_init(m@:matrix[`T], num_rows:int, num_cols:int,
			   init:&(int,int):T):m_matrix[T] {
      copy_mutable_init[T](m, num_rows, num_cols, init) }
  signature copy_mutable_init[T](matrix[T], num_rows:int, num_cols:int,
				 init:&(int,int):T):m_matrix[T];


-- m_matrix is a mutable matrix

abstract object m_matrix[T] isa matrix[T];
  signature store(matrix[`T], row:int, col:int, value:T):void;


-- vector_matrix is an implementation of mutable matrices
-- represented using nested vectors

template representation vector_matrix[T] isa m_matrix[T];
  field rows(@:vector_matrix[`T]):m_vector[m_vector[T]];

  method num_rows(m@:vector_matrix[`T]):int { m.rows.length }
  method num_cols(m@:vector_matrix[`T]):int { m.rows.fetch(0, { ^0 }).length }

  method fetch(m@:vector_matrix[`T], row:int, col:int):T { (m.rows!row)!col }
  method store(m@:vector_matrix[`T], row:int, col:int, value:T):void {
      (m.rows!row)!col := value;
  }
  method row(m@:vector_matrix[`T], row:int):indexed[T] { m.rows!row }

  method new_vector_matrix_init[T](num_rows:int, num_cols:int,
				   init:&(int,int):T):m_matrix[T] {
      concrete object isa vector_matrix[T] {
	  rows := new_m_vector_init[m_vector[T]](num_rows, &(row:int){
	      new_m_vector_init[T](num_cols, &(col:int){
		  init.eval(row, col)
	      })
	  })
      } }

  method copy_init[T](m@:vector_matrix[T], num_rows:int, num_cols:int,
		      init:&(int,int):T):matrix[T] {
      copy_mutable_init[T](m, num_rows, num_cols, init) }

  method copy_mutable_init[T](m@:vector_matrix[T], num_rows:int, num_cols:int,
			      init:&(int,int):T):m_matrix[T] {
      new_vector_matrix_init[T](num_rows, num_cols, init) }