Chris@5: 
Chris@94: module yetilab.matrix.matrix;
Chris@94: 
Chris@94: // A matrix is an array of fvectors (i.e. primitive double[]s).
Chris@94: 
Chris@94: // A matrix can be either RowMajor, akin to a C multidimensional array
Chris@96: // in which each row is a separate fvector, or ColumnMajor, akin to a
Chris@94: // FORTAN multidimensional array in which each column is a separate
Chris@96: // fvector. The default is ColumnMajor. Storage order is an efficiency
Chris@94: // concern only, all operations behave identically regardless.  (The
Chris@94: // transpose function just switches the row/column order without
Chris@94: // moving the elements.)
Chris@18: 
Chris@93: vec = load yetilab.block.fvector;
Chris@96: block = load yetilab.block.block;
Chris@9: 
Chris@96: make d = {
Chris@96:     get data () = d,
Chris@96:     get size () =
Chris@96:         case d of
Chris@96:         RowM r:
Chris@96:             major = length r;
Chris@96:             { 
Chris@96:                 rows = major, 
Chris@96:                 columns = if major > 0 then vec.length r[0] else 0 fi,
Chris@96:             };
Chris@96:         ColM c:
Chris@96:             major = length c;
Chris@96:             { 
Chris@96:                 rows = if major > 0 then vec.length c[0] else 0 fi,
Chris@96:                 columns = major, 
Chris@96:             };
Chris@94:         esac,
Chris@96:     getColumn j =
Chris@96:         case d of
Chris@96:         RowM rows: block.fromList (map do i: getAt i j done [0..length rows-1]);
Chris@96:         ColM cols: block.block cols[j];
Chris@94:         esac,
Chris@96:     getRow i =
Chris@96:         case d of
Chris@96:         RowM rows: block.block rows[i];
Chris@96:         ColM cols: block.fromList (map do j: getAt i j done [0..length cols-1]);
Chris@94:         esac,
Chris@96:     getAt row col =
Chris@96:         case d of
Chris@96:         RowM rows: r = rows[row]; (r is ~double[])[col];
Chris@96:         ColM cols: c = cols[col]; (c is ~double[])[row];
Chris@94:         esac,
Chris@96:     setAt row col n =
Chris@96:         case d of
Chris@96:         RowM rows: r = rows[row]; (r is ~double[])[col] := n;
Chris@96:         ColM cols: c = cols[col]; (c is ~double[])[row] := n;
Chris@95:         esac,
Chris@95:     get isRowMajor? () =
Chris@96:         case d of
Chris@96:         RowM _: true;
Chris@96:         ColM _: false;
Chris@95:         esac,
Chris@94:     };
Chris@94: 
Chris@94: newStorage rows cols = 
Chris@97:     if rows < 1 then array []
Chris@97:     else array (map \(vec.zeros rows) [1..cols])
Chris@97:     fi;
Chris@94: 
Chris@94: zeroMatrix rows cols = 
Chris@96:     make (ColM (newStorage rows cols));
Chris@5: 
Chris@20: generate f rows cols =
Chris@94:    (m = newStorage rows cols;
Chris@94:     for [0..cols-1] do col:
Chris@94:         for [0..rows-1] do row:
Chris@94:             m[col][row] := f row col;
Chris@5:         done;
Chris@5:     done;
Chris@96:     make (ColM m));
Chris@5: 
Chris@20: constMatrix n = generate do row col: n done;
Chris@20: randomMatrix = generate do row col: Math#random() done;
Chris@5: identityMatrix = constMatrix 1;
Chris@5: 
Chris@96: width m = m.size.columns;
Chris@96: height m = m.size.rows;
Chris@6: 
Chris@97: equal m1 m2 =
Chris@97:    (compare d1 d2 =
Chris@97:         all id (map2 vec.equal d1 d2);
Chris@97:     case m1.data of
Chris@97:     RowM d1:
Chris@97:         case m2.data of RowM d2: compare d1 d2; _: false; esac;
Chris@97:     ColM d1:
Chris@97:         case m2.data of ColM d2: compare d1 d2; _: false; esac;
Chris@97:     esac);
Chris@97: 
Chris@95: copyOf m =
Chris@95:    (copyOfData d = (array (map vec.copyOf d));
Chris@96:     make
Chris@95:        (case m.data of
Chris@96:         RowM d: RowM (copyOfData d);
Chris@96:         ColM d: ColM (copyOfData d);
Chris@95:         esac));
Chris@6: 
Chris@95: transposed m =
Chris@96:     make
Chris@95:        (case m.data of
Chris@96:         RowM d: ColM d;
Chris@96:         ColM d: RowM d;
Chris@95:         esac);
Chris@5: 
Chris@96: //!!! Change storage from column to row major order (or back
Chris@96: //again). Is there a word for this?
Chris@96: flipped m =
Chris@96:     if m.isRowMajor? then id else transposed fi
Chris@96:        (generate do row col: m.getAt col row done m.size.columns m.size.rows);
Chris@96: 
Chris@96: newMatrix type data is RowMajor () | ColumnMajor () -> list?<list?<number>> -> 'a =
Chris@96:    (tagger = case type of RowMajor (): RowM; ColumnMajor (): ColM esac;
Chris@97:     if empty? data or empty? (head data)
Chris@96:     then zeroMatrix 0 0
Chris@96:     else make (tagger (array (map vec.vector data)))
Chris@96:     fi);
Chris@96: 
Chris@96: newRowVector data = 
Chris@96:     newMatrix (RowMajor ()) [data];
Chris@96: 
Chris@96: newColumnVector data = 
Chris@96:     newMatrix (ColumnMajor ()) [data];
Chris@8: 
Chris@5: {
Chris@97: constMatrix, randomMatrix, zeroMatrix, identityMatrix,
Chris@96: generate,
Chris@96: width, height,
Chris@97: equal,
Chris@20: copyOf,
Chris@15: transposed,
Chris@96: flipped,
Chris@96: newMatrix, newRowVector, newColumnVector
Chris@5: }
Chris@5: