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@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@93: vec = load yetilab.block.fvector;
Chris@96: block = load yetilab.block.block;
Chris@99: bf = load yetilab.block.blockfuncs;
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@158:     setAt row col n = //!!! dangerous, could modify copies -- should it be allowed?
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@98: newStorage { rows, columns } = 
Chris@97:     if rows < 1 then array []
Chris@98:     else array (map \(vec.zeros rows) [1..columns])
Chris@97:     fi;
Chris@94: 
Chris@98: zeroMatrix { rows, columns } = 
Chris@98:     make (ColM (newStorage { rows, columns }));
Chris@5: 
Chris@98: generate f { rows, columns } =
Chris@98:    (m = newStorage { rows, columns };
Chris@98:     for [0..columns-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@158: zeroSizeMatrix () = zeroMatrix { rows = 0, columns = 0 };
Chris@5: 
Chris@96: width m = m.size.columns;
Chris@96: height m = m.size.rows;
Chris@6: 
Chris@100: transposed m =
Chris@100:     make
Chris@100:        (case m.data of
Chris@100:         RowM d: ColM d;
Chris@100:         ColM d: RowM d;
Chris@100:         esac);
Chris@100: 
Chris@100: flipped m =
Chris@100:     if m.isRowMajor? then
Chris@100:         generate do row col: m.getAt row col done m.size;
Chris@100:     else
Chris@100:         transposed
Chris@100:            (generate do row col: m.getAt col row done
Chris@100:             { rows = m.size.columns, columns = m.size.rows });
Chris@100:     fi;
Chris@100: 
Chris@161: toRowMajor m =
Chris@161:     if m.isRowMajor? then m else flipped m fi;
Chris@161: 
Chris@161: toColumnMajor m =
Chris@161:     if not m.isRowMajor? 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@159:     if m1.size != m2.size then false
Chris@159:     elif m1.isRowMajor? != m2.isRowMajor? then equal (flipped m1) m2;
Chris@100:     else
Chris@100:         compare d1 d2 = all id (map2 vec.equal d1 d2);
Chris@100:         case m1.data of
Chris@100:         RowM d1: case m2.data of RowM d2: compare d1 d2; _: false; esac;
Chris@100:         ColM d1: case m2.data of ColM d2: compare d1 d2; _: false; esac;
Chris@100:         esac
Chris@100:     fi;
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@163: newMatrix type data = //!!! NB does not copy data
Chris@96:    (tagger = case type of RowMajor (): RowM; ColumnMajor (): ColM esac;
Chris@163:     if empty? data or block.empty? (head data)
Chris@98:     then zeroMatrix { rows = 0, columns = 0 }
Chris@163:     else make (tagger (array (map block.data data)))
Chris@96:     fi);
Chris@96: 
Chris@163: newRowVector data = //!!! NB does not copy data
Chris@163:     make (RowM (array [block.data data]));
Chris@96: 
Chris@163: newColumnVector data = //!!! NB does not copy data
Chris@163:     make (ColM (array [block.data data]));
Chris@8: 
Chris@98: scaled factor m =
Chris@98:     generate do row col: factor * m.getAt row col done m.size;
Chris@98: 
Chris@158: resizedTo newsize m =
Chris@158:    (oldsize = m.size;
Chris@158:     if newsize == oldsize then m
Chris@158:     else 
Chris@158:         generate do row col:
Chris@158:             if row < oldsize.rows and col < oldsize.columns
Chris@158:             then m.getAt row col else 0 fi
Chris@158:             done newsize;
Chris@158:     fi);
Chris@158: 
Chris@98: sum' m1 m2 =
Chris@98:     if m1.size != m2.size
Chris@98:     then failWith "Matrices are not the same size: \(m1.size), \(m2.size)";
Chris@98:     else
Chris@98:         generate do row col: m1.getAt row col + m2.getAt row col done m1.size;
Chris@98:     fi;
Chris@98: 
Chris@98: product m1 m2 =
Chris@98:     if m1.size.columns != m2.size.rows
Chris@98:     then failWith "Matrix dimensions incompatible: \(m1.size), \(m2.size) (\(m1.size.columns != m2.size.rows)";
Chris@98:     else
Chris@98:         generate do row col:
Chris@99:             bf.sum (bf.multiply (m1.getRow row) (m2.getColumn col))
Chris@98:         done { rows = m1.size.rows, columns = m2.size.columns }
Chris@98:     fi;
Chris@98: 
Chris@161: asRows m =
Chris@171:     map m.getRow [0 .. m.size.rows - 1];
Chris@161: 
Chris@161: asColumns m =
Chris@171:     map m.getColumn [0 .. m.size.columns - 1];
Chris@161: 
Chris@178: concatAgainstGrain tagger getter counter mm =
Chris@178:    (n = counter (head mm).size;
Chris@177:     make (tagger (array
Chris@177:        (map do i:
Chris@177:            block.data (block.concat (map do m: getter m i done mm))
Chris@178:            done [0..n-1]))));
Chris@177: 
Chris@178: concatWithGrain tagger getter counter mm =
Chris@177:     make (tagger (array
Chris@177:        (concat
Chris@177:            (map do m:
Chris@178:                n = counter m.size;
Chris@178:                map do i: block.data (getter m i) done [0..n-1]
Chris@177:                done mm))));
Chris@177: 
Chris@178: checkDimensionsFor direction first mm =
Chris@178:    (counter = if direction == Horizontal () then (.rows) else (.columns) fi;
Chris@178:     n = counter first.size;
Chris@178:     if not (all id (map do m: counter m.size == n done mm)) then
Chris@178:         failWith "Matrix dimensions incompatible for concat (found \(map do m: counter m.size 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@177:     if empty? mm then zeroSizeMatrix ()
Chris@177:     else
Chris@178:         first = head mm;
Chris@178:         checkDimensionsFor direction first mm;
Chris@178:         row = first.isRowMajor?;
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@178:             if row then concatAgainstGrain RowM (.getRow) (.rows) mm;
Chris@178:             else concatWithGrain ColM (.getColumn) (.columns) mm;
Chris@178:             fi;
Chris@178:         Vertical ():
Chris@178:             if row then concatWithGrain RowM (.getRow) (.rows) mm;
Chris@178:             else concatAgainstGrain ColM (.getColumn) (.columns) mm;
Chris@178:             fi;
Chris@178:         esac;
Chris@177:     fi;
Chris@177: 
Chris@187: rowSlice start count m = //!!! doc: storage order same as input
Chris@187:     if m.isRowMajor? then
Chris@187:         make (RowM (array (map (block.data . m.getRow) [start .. start + count - 1])))
Chris@187:     else 
Chris@187:         make (ColM (array (map (block.data . (block.rangeOf start count)) (asColumns m))))
Chris@187:     fi;
Chris@187: 
Chris@187: columnSlice start count m = //!!! doc: storage order same as input
Chris@187:     if not m.isRowMajor? then
Chris@187:         make (ColM (array (map (block.data . m.getColumn) [start .. start + count - 1])))
Chris@187:     else 
Chris@187:         make (RowM (array (map (block.data . (block.rangeOf start count)) (asRows m))))
Chris@187:     fi;
Chris@187: 
Chris@5: {
Chris@158: constMatrix, randomMatrix, zeroMatrix, identityMatrix, zeroSizeMatrix,
Chris@96: generate,
Chris@96: width, height,
Chris@97: equal,
Chris@20: copyOf,
Chris@15: transposed,
Chris@161: flipped, toRowMajor, toColumnMajor,
Chris@98: scaled,
Chris@158: resizedTo,
Chris@161: asRows, asColumns,
Chris@98: sum = sum', product,
Chris@177: concat,
Chris@187: rowSlice, columnSlice,
Chris@98: newMatrix, newRowVector, newColumnVector,
Chris@5: }
Chris@5: