Chris@5: 
Chris@94: module yetilab.matrix.matrix;
Chris@94: 
Chris@214: // A matrix is an array of vectors.
Chris@94: 
Chris@101: // A matrix can be stored in either column-major (the default) or
Chris@101: // row-major format. Storage order is an efficiency concern only:
Chris@101: // every API function operating on matrix objects will return the same
Chris@101: // result regardless of storage order.  (The transpose function just
Chris@101: // switches the row/column order without moving the elements.)
Chris@18: 
Chris@214: //!!! check that we are not unnecessarily copying in the transform functions
Chris@214: 
Chris@214: vec = load yetilab.block.vector;
Chris@99: bf = load yetilab.block.blockfuncs;
Chris@9: 
Chris@214: load yetilab.block.vectortype;
Chris@195: load yetilab.matrix.matrixtype;
Chris@195: 
Chris@208: size m =
Chris@208:     case m of
Chris@208:     RowM r:
Chris@208:         major = length r;
Chris@208:         { 
Chris@208:             rows = major, 
Chris@214:             columns = if major > 0 then vec.length r[0] else 0 fi,
Chris@208:         };
Chris@208:     ColM c:
Chris@208:         major = length c;
Chris@208:         { 
Chris@214:             rows = if major > 0 then vec.length c[0] else 0 fi,
Chris@208:             columns = major, 
Chris@208:         };
Chris@208:     esac;
Chris@208: 
Chris@208: width m = (size m).columns;
Chris@208: height m = (size m).rows;
Chris@208: 
Chris@208: getAt row col m =
Chris@208:     case m of
Chris@214:     RowM rows: r = rows[row]; vec.at col r;
Chris@214:     ColM cols: c = cols[col]; vec.at row c;
Chris@208:     esac;
Chris@208: 
Chris@208: getColumn j m =
Chris@208:     case m of
Chris@214:     RowM rows: vec.fromList (map do i: getAt i j m done [0..length rows-1]);
Chris@208:     ColM cols: cols[j];
Chris@208:     esac;
Chris@208: 
Chris@208: getRow i m =
Chris@208:     case m of
Chris@208:     RowM rows: rows[i];
Chris@214:     ColM cols: vec.fromList (map do j: getAt i j m done [0..length cols-1]);
Chris@208:     esac;
Chris@208: 
Chris@214: /*
Chris@208: setAt row col n m = //!!! dangerous, could modify copies -- should it be allowed?
Chris@208:     case m of
Chris@214:     RowM rows: r = rows[row]; (vec.data r)[col] := n;
Chris@214:     ColM cols: c = cols[col]; (vec.data c)[row] := n;
Chris@208:     esac;
Chris@214: */
Chris@208: 
Chris@208: isRowMajor? m =
Chris@208:     case m of
Chris@208:     RowM _: true;
Chris@208:     ColM _: false;
Chris@208:     esac;
Chris@94: 
Chris@201: newColMajorStorage { rows, columns } = 
Chris@97:     if rows < 1 then array []
Chris@214:     else array (map \(vec.zeros rows) [1..columns])
Chris@97:     fi;
Chris@94: 
Chris@98: zeroMatrix { rows, columns } = 
Chris@208:     ColM (newColMajorStorage { rows, columns });
Chris@201: 
Chris@201: zeroMatrixWithTypeOf m { rows, columns } = 
Chris@208:     if isRowMajor? m then
Chris@208:         RowM (newColMajorStorage { rows = columns, columns = rows });
Chris@201:     else
Chris@208:         ColM (newColMajorStorage { rows, columns });
Chris@201:     fi;
Chris@5: 
Chris@214: zeroSizeMatrix () = zeroMatrix { rows = 0, columns = 0 };
Chris@214: 
Chris@98: generate f { rows, columns } =
Chris@214:     if rows < 1 or columns < 1 then zeroSizeMatrix ()
Chris@214:     else
Chris@214:         m = array (map \(new double[rows]) [1..columns]);
Chris@214:         for [0..columns-1] do col:
Chris@214:             for [0..rows-1] do row:
Chris@214:                 m[col][row] := f row col;
Chris@214:             done;
Chris@5:         done;
Chris@214:         ColM (array (map vec.vector m))
Chris@214:     fi;
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@100: transposed m =
Chris@208:     case m of
Chris@208:     RowM d: ColM d;
Chris@208:     ColM d: RowM d;
Chris@208:     esac;
Chris@100: 
Chris@100: flipped m =
Chris@208:     if isRowMajor? m then
Chris@208:         generate do row col: getAt row col m done (size m);
Chris@100:     else
Chris@100:         transposed
Chris@208:            (generate do row col: getAt col row m done
Chris@208:             { rows = (width m), columns = (height m) });
Chris@100:     fi;
Chris@100: 
Chris@161: toRowMajor m =
Chris@208:     if isRowMajor? m then m else flipped m fi;
Chris@161: 
Chris@161: toColumnMajor m =
Chris@208:     if not isRowMajor? m then m else flipped m fi;
Chris@161: 
Chris@100: // Matrices with different storage order but the same contents are
Chris@100: // equal (but comparing them is slow)
Chris@97: equal m1 m2 =
Chris@208:     if size m1 != size m2 then false
Chris@208:     elif isRowMajor? m1 != isRowMajor? m2 then equal (flipped m1) m2;
Chris@100:     else
Chris@214:         compare d1 d2 = all id (map2 vec.equal d1 d2);
Chris@208:         case m1 of
Chris@208:         RowM d1: case m2 of RowM d2: compare d1 d2; _: false; esac;
Chris@208:         ColM d1: case m2 of ColM d2: compare d1 d2; _: false; esac;
Chris@100:         esac
Chris@100:     fi;
Chris@97: 
Chris@214: /*!!! not needed now it's immutable?
Chris@95: copyOf m =
Chris@214:    (copyOfData d = (array (map vec.copyOf d));
Chris@208:     case m of
Chris@208:     RowM d: RowM (copyOfData d);
Chris@208:     ColM d: ColM (copyOfData d);
Chris@208:     esac);
Chris@214: */
Chris@6: 
Chris@163: newMatrix type data = //!!! NB does not copy data
Chris@96:    (tagger = case type of RowMajor (): RowM; ColumnMajor (): ColM esac;
Chris@214:     if empty? data or vec.empty? (head data)
Chris@201:     then zeroSizeMatrix ()
Chris@208:     else tagger (array data)
Chris@96:     fi);
Chris@96: 
Chris@163: newRowVector data = //!!! NB does not copy data
Chris@208:     RowM (array [data]);
Chris@96: 
Chris@163: newColumnVector data = //!!! NB does not copy data
Chris@208:     ColM (array [data]);
Chris@8: 
Chris@195: scaled factor m = //!!! v inefficient
Chris@208:     generate do row col: factor * (getAt row col m) done (size m);
Chris@98: 
Chris@98: sum' m1 m2 =
Chris@208:     if (size m1) != (size m2)
Chris@208:     then failWith "Matrices are not the same size: \(size m1), \(size m2)";
Chris@98:     else
Chris@208:         generate do row col: getAt row col m1 + getAt row col m2 done (size m1);
Chris@98:     fi;
Chris@98: 
Chris@98: product m1 m2 =
Chris@208:     if (size m1).columns != (size m2).rows
Chris@208:     then failWith "Matrix dimensions incompatible: \(size m1), \(size m2) (\((size m1).columns != (size m2).rows)";
Chris@98:     else
Chris@98:         generate do row col:
Chris@208:             bf.sum (bf.multiply (getRow row m1) (getColumn col m2))
Chris@208:         done { rows = (size m1).rows, columns = (size m2).columns }
Chris@98:     fi;
Chris@98: 
Chris@161: asRows m =
Chris@208:     map do i: getRow i m done [0 .. (height m) - 1];
Chris@161: 
Chris@161: asColumns m =
Chris@208:     map do i: getColumn i m done [0 .. (width m) - 1];
Chris@161: 
Chris@178: concatAgainstGrain tagger getter counter mm =
Chris@208:    (n = counter (size (head mm));
Chris@208:     tagger (array
Chris@177:        (map do i:
Chris@214:            vec.concat (map (getter i) mm)
Chris@208:            done [0..n-1])));
Chris@177: 
Chris@178: concatWithGrain tagger getter counter mm =
Chris@208:     tagger (array
Chris@177:        (concat
Chris@177:            (map do m:
Chris@208:                n = counter (size m);
Chris@208:                map do i: getter i m done [0..n-1]
Chris@208:                done mm)));
Chris@177: 
Chris@178: checkDimensionsFor direction first mm =
Chris@178:    (counter = if direction == Horizontal () then (.rows) else (.columns) fi;
Chris@208:     n = counter (size first);
Chris@208:     if not (all id (map do m: counter (size m) == n done mm)) then
Chris@208:         failWith "Matrix dimensions incompatible for concat (found \(map do m: counter (size m) done mm) not all of which are \(n))";
Chris@178:     fi);
Chris@178: 
Chris@187: concat direction mm = //!!! doc: storage order is taken from first matrix in sequence
Chris@187:     //!!! would this be better as separate concatHorizontal/concatVertical functions?
Chris@190:     case mm of
Chris@190:     first::rest: 
Chris@178:         checkDimensionsFor direction first mm;
Chris@208:         row = isRowMajor? first;
Chris@178:         // horizontal, row-major: against grain with rows
Chris@178:         // horizontal, col-major: with grain with cols
Chris@178:         // vertical, row-major: with grain with rows
Chris@178:         // vertical, col-major: against grain with cols
Chris@178:         case direction of
Chris@178:         Horizontal ():
Chris@208:             if row then concatAgainstGrain RowM getRow (.rows) mm;
Chris@208:             else concatWithGrain ColM getColumn (.columns) mm;
Chris@178:             fi;
Chris@178:         Vertical ():
Chris@208:             if row then concatWithGrain RowM getRow (.rows) mm;
Chris@208:             else concatAgainstGrain ColM getColumn (.columns) mm;
Chris@178:             fi;
Chris@178:         esac;
Chris@190:     [single]: single;
Chris@190:     _: zeroSizeMatrix ();
Chris@190:     esac;
Chris@177: 
Chris@187: rowSlice start count m = //!!! doc: storage order same as input
Chris@208:     if isRowMajor? m then
Chris@208:         RowM (array (map ((flip getRow) m) [start .. start + count - 1]))
Chris@187:     else 
Chris@214:         ColM (array (map (vec.rangeOf start count) (asColumns m)))
Chris@187:     fi;
Chris@187: 
Chris@187: columnSlice start count m = //!!! doc: storage order same as input
Chris@208:     if not isRowMajor? m then
Chris@208:         ColM (array (map ((flip getColumn) m) [start .. start + count - 1]))
Chris@187:     else 
Chris@214:         RowM (array (map (vec.rangeOf start count) (asRows m)))
Chris@187:     fi;
Chris@187: 
Chris@201: resizedTo newsize m =
Chris@208:    (if newsize == (size m) then
Chris@201:         m
Chris@208:     elif (height m) == 0 or (width m) == 0 then
Chris@202:         zeroMatrixWithTypeOf m newsize;
Chris@201:     else
Chris@208:         growrows = newsize.rows - (height m);
Chris@208:         growcols = newsize.columns - (width m);
Chris@208:         rowm = isRowMajor? m;
Chris@201:         resizedTo newsize
Chris@201:             if rowm and growrows < 0 then
Chris@201:                 rowSlice 0 newsize.rows m
Chris@201:             elif (not rowm) and growcols < 0 then 
Chris@201:                 columnSlice 0 newsize.columns m
Chris@201:             elif growrows < 0 then 
Chris@201:                 rowSlice 0 newsize.rows m
Chris@201:             elif growcols < 0 then 
Chris@201:                 columnSlice 0 newsize.columns m
Chris@201:             else
Chris@201:                 if growrows > 0 then
Chris@201:                     concat (Vertical ())
Chris@208:                        [m, zeroMatrixWithTypeOf m ((size m) with { rows = growrows })]
Chris@201:                 else
Chris@201:                     concat (Horizontal ())
Chris@208:                        [m, zeroMatrixWithTypeOf m ((size m) with { columns = growcols })]
Chris@201:                 fi
Chris@201:             fi
Chris@202:     fi);
Chris@201: 
Chris@5: {
Chris@208:     size,
Chris@208:     width,
Chris@208:     height,
Chris@208:     getAt,
Chris@208:     getColumn,
Chris@208:     getRow,
Chris@214: //    setAt,
Chris@208:     isRowMajor?,
Chris@208:     generate,
Chris@208:     constMatrix,
Chris@208:     randomMatrix,
Chris@208:     zeroMatrix,
Chris@208:     identityMatrix,
Chris@208:     zeroSizeMatrix,
Chris@208:     equal,
Chris@214: //    copyOf,
Chris@208:     transposed,
Chris@208:     flipped,
Chris@208:     toRowMajor,
Chris@208:     toColumnMajor,
Chris@208:     scaled,
Chris@208:     resizedTo,
Chris@208:     asRows,
Chris@208:     asColumns,
Chris@208:     sum = sum',
Chris@208:     product,
Chris@208:     concat,
Chris@208:     rowSlice,
Chris@208:     columnSlice,
Chris@208:     newMatrix,
Chris@208:     newRowVector,
Chris@208:     newColumnVector,
Chris@208: }
Chris@208: as
Chris@208: {
Chris@208: //!!! check whether these are right to be .selector rather than just selector
Chris@208: 
Chris@208:     size is matrix -> { .rows is number, .columns is number },
Chris@208:     width is matrix -> number,
Chris@208:     height is matrix -> number,
Chris@208:     getAt is number -> number -> matrix -> number,
Chris@214:     getColumn is number -> matrix -> vector,
Chris@214:     getRow is number -> matrix -> vector,
Chris@214: //    setAt is number -> number -> number -> matrix -> (), //!!! lose?
Chris@208:     isRowMajor? is matrix -> boolean,
Chris@195:     generate is (number -> number -> number) -> { .rows is number, .columns is number } -> matrix,
Chris@195:     constMatrix is number -> { .rows is number, .columns is number } -> matrix,
Chris@195:     randomMatrix is { .rows is number, .columns is number } -> matrix,
Chris@195:     zeroMatrix is { .rows is number, .columns is number } -> matrix, 
Chris@195:     identityMatrix is { .rows is number, .columns is number } -> matrix, 
Chris@195:     zeroSizeMatrix is () -> matrix,
Chris@195:     equal is matrix -> matrix -> boolean,
Chris@214: //    copyOf is matrix -> matrix,
Chris@195:     transposed is matrix -> matrix,
Chris@195:     flipped is matrix -> matrix, 
Chris@195:     toRowMajor is matrix -> matrix, 
Chris@195:     toColumnMajor is matrix -> matrix,
Chris@195:     scaled is number -> matrix -> matrix,
Chris@195:     resizedTo is { .rows is number, .columns is number } -> matrix -> matrix,
Chris@214:     asRows is matrix -> list<vector>, 
Chris@214:     asColumns is matrix -> list<vector>,
Chris@208:     sum is matrix -> matrix -> matrix,
Chris@195:     product is matrix -> matrix -> matrix,
Chris@195:     concat is (Horizontal () | Vertical ()) -> list<matrix> -> matrix,
Chris@195:     rowSlice is number -> number -> matrix -> matrix, 
Chris@195:     columnSlice is number -> number -> matrix -> matrix,
Chris@214:     newMatrix is (ColumnMajor () | RowMajor ()) -> list<vector> -> matrix, 
Chris@214:     newRowVector is vector -> matrix, 
Chris@214:     newColumnVector is vector -> matrix,
Chris@5: }
Chris@5: