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@222: vec = load yetilab.vector.vector;
Chris@222: bf = load yetilab.vector.blockfuncs;
Chris@9: 
Chris@222: load yetilab.vector.vectortype;
Chris@195: load yetilab.matrix.matrixtype;
Chris@195: 
Chris@208: size m =
Chris@208:     case m of
Chris@234:     DenseRows 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@234:     DenseCols 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@237:     SparseCSR { values, indices, pointers, extent }:
Chris@234:         {
Chris@236:             rows = (length pointers) - 1,
Chris@237:             columns = extent
Chris@234:         };
Chris@237:     SparseCSC { values, indices, pointers, extent }:
Chris@234:         {
Chris@237:             rows = extent,
Chris@236:             columns = (length pointers) - 1
Chris@234:         };
Chris@208:     esac;
Chris@208: 
Chris@208: width m = (size m).columns;
Chris@208: height m = (size m).rows;
Chris@208: 
Chris@249: nonZeroValues m =
Chris@249:    (nz d =
Chris@242:         sum
Chris@242:            (map do v:
Chris@242:                 sum (map do n: if n == 0 then 0 else 1 fi done (vec.list v))
Chris@242:                 done d);
Chris@242:     case m of 
Chris@249:     DenseRows d: nz d;
Chris@249:     DenseCols d: nz d;
Chris@249:     SparseCSR d: vec.length d.values;
Chris@249:     SparseCSC d: vec.length d.values;
Chris@242:     esac);
Chris@242: 
Chris@249: density m =
Chris@249:    ({ rows, columns } = size m;
Chris@249:     cells = rows * columns;
Chris@249:     (nonZeroValues m) / cells);
Chris@249: 
Chris@236: sparseSlice n d =
Chris@236:    (start = d.pointers[n];
Chris@236:     end = d.pointers[n+1];
Chris@236:     { 
Chris@236:         values = vec.slice d.values start end,
Chris@239:         indices = slice d.indices start end,
Chris@236:     });
Chris@236: 
Chris@252: nonEmptySlices d =
Chris@252:    (ne = array [];
Chris@252:     for [0..length d.pointers - 2] do i:
Chris@252:         if d.pointers[i] != d.pointers[i+1] then
Chris@252:             push ne i
Chris@252:         fi
Chris@252:     done;
Chris@252:     ne);
Chris@252: 
Chris@236: fromSlice n m d =
Chris@236:    (slice = sparseSlice n d;
Chris@236:     var v = 0;
Chris@236:     for [0..length slice.indices - 1] do i:
Chris@236:         if slice.indices[i] == m then
Chris@260:             v := vec.at slice.values i;
Chris@236:         fi
Chris@236:     done;
Chris@236:     v);
Chris@236: 
Chris@236: filledSlice n d =
Chris@236:    (slice = sparseSlice n d;
Chris@236:     dslice = new double[d.extent];
Chris@239:     for [0..length slice.indices - 1] do i:
Chris@260:         dslice[slice.indices[i]] := vec.at slice.values i;
Chris@239:     done;
Chris@236:     vec.vector dslice);
Chris@236: 
Chris@260: at' m row col =
Chris@208:     case m of
Chris@260:     DenseRows rows: r = rows[row]; vec.at r col;
Chris@260:     DenseCols cols: c = cols[col]; vec.at c row;
Chris@236:     SparseCSR data: fromSlice row col data;
Chris@236:     SparseCSC data: fromSlice col row data;
Chris@208:     esac;
Chris@208: 
Chris@208: getColumn j m =
Chris@208:     case m of
Chris@234:     DenseCols cols: cols[j];
Chris@236:     SparseCSC data: filledSlice j data;
Chris@260:     _: vec.fromList (map do i: at' m i j done [0..height m - 1]);
Chris@208:     esac;
Chris@208: 
Chris@208: getRow i m =
Chris@208:     case m of
Chris@234:     DenseRows rows: rows[i];
Chris@236:     SparseCSR data: filledSlice i data; 
Chris@260:     _: vec.fromList (map do j: at' m i j done [0..width m - 1]);
Chris@208:     esac;
Chris@208: 
Chris@257: asRows m =
Chris@257:     map do i: getRow i m done [0 .. (height m) - 1];
Chris@257: 
Chris@257: asColumns m =
Chris@257:     map do i: getColumn i m done [0 .. (width m) - 1];
Chris@257: 
Chris@208: isRowMajor? m =
Chris@208:     case m of
Chris@234:     DenseRows _: true;
Chris@234:     DenseCols _: false;
Chris@234:     SparseCSR _: true;
Chris@234:     SparseCSC _: false;
Chris@234:     esac;
Chris@234: 
Chris@234: isSparse? m =
Chris@234:     case m of
Chris@234:     DenseRows _: false;
Chris@234:     DenseCols _: false;
Chris@234:     SparseCSR _: true;
Chris@234:     SparseCSC _: true;
Chris@208:     esac;
Chris@94: 
Chris@261: typeOf m =
Chris@261:     if isRowMajor? m then RowMajor ()
Chris@261:     else ColumnMajor ()
Chris@261:     fi;
Chris@261: 
Chris@261: flippedTypeOf m =
Chris@261:     if isRowMajor? m then ColumnMajor ()
Chris@261:     else RowMajor ()
Chris@261:     fi;
Chris@261: 
Chris@244: newColumnMajorStorage { 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@244:     DenseCols (newColumnMajorStorage { rows, columns });
Chris@201: 
Chris@201: zeroMatrixWithTypeOf m { rows, columns } = 
Chris@208:     if isRowMajor? m then
Chris@244:         DenseRows (newColumnMajorStorage { rows = columns, columns = rows });
Chris@201:     else
Chris@244:         DenseCols (newColumnMajorStorage { 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@234:         DenseCols (array (map vec.vector m))
Chris@214:     fi;
Chris@5: 
Chris@255: swapij =
Chris@255:     map do { i, j, v }: { i = j, j = i, v } done;
Chris@255: 
Chris@257: //!!! should use { row = , column = , value = } instead of i, j, v?
Chris@235: enumerateSparse m =
Chris@240:    (enumerate { values, indices, pointers } =
Chris@240:         concat
Chris@240:            (map do i:
Chris@240:                 start = pointers[i];
Chris@240:                 end = pointers[i+1];
Chris@240:                 map2 do j v: { i, j, v } done 
Chris@240:                     (slice indices start end)
Chris@240:                     (vec.list (vec.slice values start end))
Chris@240:                 done [0..length pointers - 2]);
Chris@235:     case m of
Chris@255:     SparseCSC d: swapij (enumerate d);
Chris@255:     SparseCSR d: enumerate d;
Chris@237:      _: [];
Chris@235:     esac);
Chris@235: 
Chris@255: enumerateDense m =
Chris@255:    (enumerate d =
Chris@255:         concat
Chris@255:            (map do i:
Chris@255:                 vv = d[i];
Chris@255:                 map2 do j v: { i, j, v } done
Chris@255:                     [0..vec.length vv - 1]
Chris@255:                     (vec.list vv);
Chris@255:                 done [0..length d - 1]);
Chris@255:     case m of
Chris@255:     DenseCols c: swapij (enumerate c);
Chris@255:     DenseRows r: enumerate r;
Chris@255:      _: [];
Chris@255:     esac);
Chris@255: 
Chris@255: enumerate m =
Chris@255:     if isSparse? m then enumerateSparse m else enumerateDense m fi;
Chris@255: 
Chris@261: // Make a sparse matrix from entries whose i, j values are known to be
Chris@261: // within range
Chris@238: makeSparse type size data =
Chris@244:    (isRow = case type of RowMajor (): true; ColumnMajor (): false esac;
Chris@238:     ordered = 
Chris@238:         sortBy do a b:
Chris@238:             if a.maj == b.maj then a.min < b.min else a.maj < b.maj fi
Chris@238:         done
Chris@238:            (map
Chris@238:                 if isRow then
Chris@238:                     do { i, j, v }: { maj = i, min = j, v } done;
Chris@238:                 else
Chris@238:                     do { i, j, v }: { maj = j, min = i, v } done;
Chris@238:                 fi
Chris@255:                (filter do d: d.v != 0 done data));
Chris@238:     tagger = if isRow then SparseCSR else SparseCSC fi;
Chris@238:     majorSize = if isRow then size.rows else size.columns fi;
Chris@238:     minorSize = if isRow then size.columns else size.rows fi;
Chris@251:     pointers = array [0];
Chris@251:     setArrayCapacity pointers (size.rows + 1);
Chris@251:     fillPointers n i data =
Chris@251:         if n < majorSize then
Chris@251:             case data of
Chris@238:             d::rest:
Chris@251:                (for [n..d-1] \(push pointers i);
Chris@251:                 fillPointers d (i+1) rest);
Chris@251:              _:
Chris@251:                 for [n..majorSize-1] \(push pointers i);
Chris@238:             esac;
Chris@238:         fi;
Chris@251:     fillPointers 0 0 (map (.maj) ordered);
Chris@238:     tagger {
Chris@238:         values = vec.fromList (map (.v) ordered),
Chris@238:         indices = array (map (.min) ordered),
Chris@251:         pointers,
Chris@238:         extent = minorSize,
Chris@238:     });
Chris@238: 
Chris@261: // Make a sparse matrix from entries that may contain out-of-range
Chris@261: // cells which need to be filtered out. This is the public API for
Chris@261: // makeSparse and is also used to discard out-of-range cells from
Chris@261: // resizedTo.
Chris@261: newSparseMatrix type size data =
Chris@261:     makeSparse type size
Chris@261:        (filter
Chris@261:             do { i, j, v }:
Chris@261:                 i == int i and i >= 0 and i < size.rows and 
Chris@261:                 j == int j and j >= 0 and j < size.columns
Chris@261:             done data);
Chris@261: 
Chris@243: toSparse m =
Chris@238:     if isSparse? m then m
Chris@238:     else
Chris@261:         makeSparse (typeOf m) (size m) (enumerateDense m);
Chris@238:     fi;
Chris@238: 
Chris@238: toDense m =
Chris@238:     if not (isSparse? m) then m
Chris@238:     elif isRowMajor? m then
Chris@238:         DenseRows (array (map do row: getRow row m done [0..height m - 1]));
Chris@238:     else
Chris@238:         DenseCols (array (map do col: getColumn col m done [0..width m - 1]));
Chris@238:     fi;
Chris@235: 
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@234:     DenseRows d: DenseCols d;
Chris@234:     DenseCols d: DenseRows d;
Chris@236:     SparseCSR d: SparseCSC d;
Chris@236:     SparseCSC d: SparseCSR d;
Chris@208:     esac;
Chris@100: 
Chris@100: flipped m =
Chris@235:     if isSparse? m then
Chris@261:         makeSparse (flippedTypeOf m) (size m) (enumerateSparse m)
Chris@100:     else
Chris@235:         if isRowMajor? m then
Chris@260:             generate do row col: at' m row col done (size m);
Chris@235:         else
Chris@235:             transposed
Chris@260:                (generate do row col: at' m col row done
Chris@235:                 { rows = (width m), columns = (height m) });
Chris@235:         fi
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@238: equal'' comparator vecComparator m1 m2 =
Chris@238:     // Prerequisite: m1 and m2 have same sparse-p and storage order
Chris@241:    (compareVecLists vv1 vv2 = all id (map2 vecComparator vv1 vv2);
Chris@238:     compareSparse d1 d2 =
Chris@238:         d1.extent == d2.extent and
Chris@238:         vecComparator d1.values d2.values and
Chris@241:         d1.indices == d2.indices and
Chris@241:         d1.pointers == d2.pointers;
Chris@238:     case m1 of
Chris@238:     DenseRows d1:
Chris@238:         case m2 of DenseRows d2: compareVecLists d1 d2; _: false; esac;
Chris@238:     DenseCols d1:
Chris@238:         case m2 of DenseCols d2: compareVecLists d1 d2; _: false; esac;
Chris@238:     SparseCSR d1:
Chris@238:         case m2 of SparseCSR d2: compareSparse d1 d2; _: false; esac;
Chris@238:     SparseCSC d1:
Chris@238:         case m2 of SparseCSC d2: compareSparse d1 d2; _: false; esac;
Chris@238:     esac);
Chris@238: 
Chris@238: equal' comparator vecComparator m1 m2 =
Chris@226:     if size m1 != size m2 then 
Chris@226:         false
Chris@226:     elif isRowMajor? m1 != isRowMajor? m2 then
Chris@238:         equal' comparator vecComparator (flipped m1) m2;
Chris@238:     elif isSparse? m1 != isSparse? m2 then
Chris@238:         if isSparse? m1 then
Chris@243:             equal' comparator vecComparator m1 (toSparse m2)
Chris@238:         else
Chris@243:             equal' comparator vecComparator (toSparse m1) m2
Chris@238:         fi
Chris@100:     else
Chris@238:         equal'' comparator vecComparator m1 m2
Chris@100:     fi;
Chris@97: 
Chris@228: // Compare matrices using the given comparator for individual cells.
Chris@228: // Note that matrices with different storage order but the same
Chris@228: // contents are equal, although comparing them is slow.
Chris@249: //!!! Document the fact that sparse matrices can only be equal if they
Chris@249: // have the same set of non-zero cells (regardless of comparator used)
Chris@228: equalUnder comparator =
Chris@238:     equal' comparator (vec.equalUnder comparator);
Chris@228: 
Chris@228: equal =
Chris@238:     equal' (==) vec.equal;
Chris@226: 
Chris@163: newMatrix type data = //!!! NB does not copy data
Chris@238:    (tagger = case type of RowMajor (): DenseRows; ColumnMajor (): DenseCols 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@234:     DenseRows (array [data]);
Chris@96: 
Chris@163: newColumnVector data = //!!! NB does not copy data
Chris@234:     DenseCols (array [data]);
Chris@8: 
Chris@257: denseLinearOp op m1 m2 =
Chris@257:     if isRowMajor? m1 then
Chris@261:         newMatrix (typeOf m1) 
Chris@257:            (map2 do c1 c2: op c1 c2 done (asRows m1) (asRows m2));
Chris@257:     else
Chris@261:         newMatrix (typeOf m1) 
Chris@257:            (map2 do c1 c2: op c1 c2 done (asColumns m1) (asColumns m2));
Chris@257:     fi;
Chris@257: 
Chris@257: sparseSumOrDifference op m1 m2 =
Chris@257:    (h = [:];
Chris@257:     for (enumerate m1) do { i, j, v }:
Chris@257:         if not (i in h) then h[i] := [:] fi;
Chris@257:         h[i][j] := v;
Chris@257:     done;
Chris@257:     for (enumerate m2) do { i, j, v }:
Chris@257:         if not (i in h) then h[i] := [:] fi;
Chris@257:         if j in h[i] then h[i][j] := op h[i][j] v;
Chris@257:         else h[i][j] := op 0 v;
Chris@257:         fi;
Chris@257:     done;
Chris@257:     entries = concat
Chris@257:        (map do i:
Chris@257:             kk = keys h[i];
Chris@257:             map2 do j v: { i, j, v } done kk (map (at h[i]) kk)
Chris@257:             done (keys h));
Chris@261:     makeSparse (typeOf m1) (size m1) entries);
Chris@257: 
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@257:     elif isSparse? m1 and isSparse? m2 then
Chris@257:         sparseSumOrDifference (+) m1 m2;
Chris@98:     else
Chris@257:         denseLinearOp bf.add m1 m2;
Chris@98:     fi;
Chris@98: 
Chris@257: difference m1 m2 =
Chris@229:     if (size m1) != (size m2)
Chris@229:     then failWith "Matrices are not the same size: \(size m1), \(size m2)";
Chris@257:     elif isSparse? m1 and isSparse? m2 then
Chris@257:         sparseSumOrDifference (-) m1 m2;
Chris@229:     else
Chris@257:         denseLinearOp bf.subtract m1 m2;
Chris@257:     fi;
Chris@257: 
Chris@257: scaled factor m =
Chris@257:     if isSparse? m then
Chris@261:         makeSparse (typeOf m) (size m)
Chris@261:            (map do { i, j, v }: { i, j, v = factor * v } done (enumerate m))
Chris@257:     elif isRowMajor? m then
Chris@261:         newMatrix (typeOf m) (map (bf.scaled factor) (asRows m));
Chris@257:     else
Chris@261:         newMatrix (typeOf m) (map (bf.scaled factor) (asColumns m));
Chris@229:     fi;
Chris@229: 
Chris@229: abs' m =
Chris@257:     if isSparse? m then
Chris@261:         makeSparse (typeOf m) (size m)
Chris@261:            (map do { i, j, v }: { i, j, v = abs v } done (enumerate m))
Chris@257:     elif isRowMajor? m then
Chris@261:         newMatrix (typeOf m) (map bf.abs (asRows m));
Chris@257:     else
Chris@261:         newMatrix (typeOf m) (map bf.abs (asColumns m));
Chris@257:     fi;
Chris@229: 
Chris@258: filter f m =
Chris@258:     if isSparse? m then
Chris@261:         makeSparse (typeOf m) (size m)
Chris@261:            (map do { i, j, v }: { i, j, v = if f v then v else 0 fi } done
Chris@261:                (enumerate m))
Chris@258:     else
Chris@258:         vfilter = vec.fromList . (map do i: if f i then i else 0 fi done) . vec.list;
Chris@258:         if isRowMajor? m then
Chris@261:             newMatrix (typeOf m) (map vfilter (asRows m));
Chris@258:         else
Chris@261:             newMatrix (typeOf m) (map vfilter (asColumns m));
Chris@258:         fi;
Chris@258:     fi;
Chris@258: 
Chris@249: sparseProductLeft size m1 m2 =
Chris@271:    (d1 = case m1 of
Chris@271:          SparseCSR d: d; 
Chris@271:          SparseCSC d: d;
Chris@271:          _: failWith "sparseProductLeft called for non-sparse m1";
Chris@271:          esac;
Chris@271:     swap = not isRowMajor? m1;
Chris@249:     data = array (map \(new double[size.rows]) [1..size.columns]);
Chris@271:     vv = d1.values;
Chris@271:     ii = d1.indices;
Chris@271:     pp = d1.pointers;
Chris@249:     for [0..size.columns - 1] do j':
Chris@249:         c = getColumn j' m2;
Chris@271:         var i = 0;
Chris@271:         for [0..length ii - 1] do ix:
Chris@271:             j = ii[ix];
Chris@271:             ix == pp[i+1] loop (i := i + 1);
Chris@271:             if swap then
Chris@271:                 data[j'][j] := data[j'][j] + (vec.at vv ix) * (vec.at c i);
Chris@271:             else
Chris@271:                 data[j'][i] := data[j'][i] + (vec.at vv ix) * (vec.at c j);
Chris@271:             fi;
Chris@249:         done;
Chris@249:     done;
Chris@249:     DenseCols (array (map vec.vector (list data))));
Chris@249: 
Chris@249: sparseProductRight size m1 m2 =
Chris@249:    (e = enumerateSparse m2;
Chris@249:     data = array (map \(new double[size.columns]) [1..size.rows]);
Chris@249:     for [0..size.rows - 1] do i':
Chris@249:         r = getRow i' m1;
Chris@249:         for e do { v, i, j }:
Chris@260:             data[i'][j] := data[i'][j] + v * (vec.at r i);
Chris@249:         done;
Chris@249:     done;
Chris@249:     DenseRows (array (map vec.vector (list data))));
Chris@249: 
Chris@251: sparseProduct size m1 m2 =
Chris@251:     case m2 of
Chris@251:     SparseCSC d:
Chris@251:        (e = enumerateSparse m1;
Chris@251:         entries =
Chris@251:            (map do j':
Chris@251:                 cs = sparseSlice j' d;
Chris@252:                 hin = mapIntoHash
Chris@260:                    (at cs.indices) (vec.at cs.values)
Chris@252:                    [0..length cs.indices - 1];
Chris@252:                 hout = [:];
Chris@252:                 for e do { v, i, j }:
Chris@252:                     if j in hin then
Chris@252:                         p = v * hin[j];
Chris@252:                         hout[i] := p + (if i in hout then hout[i] else 0 fi);
Chris@252:                     fi
Chris@252:                 done;
Chris@252:                 map do i:
Chris@252:                     { i, j = j', v = hout[i] }
Chris@252:                 done (keys hout);
Chris@252:             done (nonEmptySlices d));
Chris@251:         makeSparse (ColumnMajor ()) size (concat entries));
Chris@251:     SparseCSR _:
Chris@251:         sparseProduct size m1 (flipped m2);
Chris@251:      _: failWith "sparseProduct called for non-sparse matrices";
Chris@251:     esac;
Chris@251: 
Chris@249: denseProduct size m1 m2 =
Chris@249:    (data = array (map \(new double[size.rows]) [1..size.columns]);
Chris@249:     for [0..size.rows - 1] do i:
Chris@249:         row = getRow i m1;
Chris@249:         for [0..size.columns - 1] do j:
Chris@249:             data[j][i] := bf.sum (bf.multiply row (getColumn j m2));
Chris@249:         done;
Chris@249:     done;
Chris@249:     DenseCols (array (map vec.vector (list data))));
Chris@249: 
Chris@98: product m1 m2 =
Chris@208:     if (size m1).columns != (size m2).rows
Chris@246:     then failWith "Matrix dimensions incompatible: \(size m1), \(size m2) (\((size m1).columns) != \((size m2).rows))";
Chris@249:     else 
Chris@249:         size = { rows = (size m1).rows, columns = (size m2).columns };
Chris@249:         if isSparse? m1 then
Chris@251:             if isSparse? m2 then
Chris@251:                 sparseProduct size m1 m2
Chris@251:             else
Chris@251:                 sparseProductLeft size m1 m2
Chris@251:             fi
Chris@249:         elif isSparse? m2 then
Chris@249:             sparseProductRight size m1 m2
Chris@249:         else
Chris@249:             denseProduct size m1 m2
Chris@249:         fi;
Chris@98:     fi;
Chris@98: 
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@259: sparseConcat direction first mm =
Chris@259:    (dimension d f = if direction == d then sum (map f mm) else f first fi;
Chris@259:     rows = dimension (Vertical ()) height;
Chris@259:     columns = dimension (Horizontal ()) width;
Chris@259:     entries ioff joff ui uj mm =
Chris@259:         case mm of 
Chris@259:         m::rest:
Chris@259:            (map do { i, j, v }: { i = i + ioff, j = j + joff, v }
Chris@259:                 done (enumerate m)) ++
Chris@259:            (entries
Chris@259:                (ioff + ui * height m)
Chris@259:                (joff + uj * width m)
Chris@259:                 ui uj rest);
Chris@259:          _: []
Chris@259:         esac;
Chris@261:     makeSparse (typeOf first) { rows, columns }
Chris@259:         if direction == Vertical () then entries 0 0 1 0 mm
Chris@259:         else entries 0 0 0 1 mm fi);
Chris@259: 
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@259:     case length mm of
Chris@259:     0: zeroSizeMatrix ();
Chris@259:     1: head mm;
Chris@259:     _:
Chris@259:         first = head mm;
Chris@178:         checkDimensionsFor direction first mm;
Chris@259:         if all isSparse? mm then
Chris@259:             sparseConcat direction first mm
Chris@259:         else
Chris@259:             row = isRowMajor? first;
Chris@259:             // horizontal, row-major: against grain with rows
Chris@259:             // horizontal, col-major: with grain with cols
Chris@259:             // vertical, row-major: with grain with rows
Chris@259:             // vertical, col-major: against grain with cols
Chris@259:             case direction of
Chris@259:             Horizontal ():
Chris@259:                 if row then concatAgainstGrain DenseRows getRow (.rows) mm;
Chris@259:                 else concatWithGrain DenseCols getColumn (.columns) mm;
Chris@259:                 fi;
Chris@259:             Vertical ():
Chris@259:                 if row then concatWithGrain DenseRows getRow (.rows) mm;
Chris@259:                 else concatAgainstGrain DenseCols getColumn (.columns) mm;
Chris@259:                 fi;
Chris@259:             esac;
Chris@259:         fi;
Chris@190:     esac;
Chris@177: 
Chris@260: //!!! doc note: argument order chosen for consistency with std module slice
Chris@260: rowSlice m start end = //!!! doc: storage order same as input
Chris@208:     if isRowMajor? m then
Chris@260:         DenseRows (array (map ((flip getRow) m) [start .. end - 1]))
Chris@187:     else 
Chris@260:         DenseCols (array (map do v: vec.slice v start end done (asColumns m)))
Chris@187:     fi;
Chris@187: 
Chris@260: //!!! doc note: argument order chosen for consistency with std module slice
Chris@260: columnSlice m start end = //!!! doc: storage order same as input
Chris@208:     if not isRowMajor? m then
Chris@260:         DenseCols (array (map ((flip getColumn) m) [start .. end - 1]))
Chris@187:     else 
Chris@260:         DenseRows (array (map do v: vec.slice v start end done (asRows m)))
Chris@187:     fi;
Chris@187: 
Chris@201: resizedTo newsize m =
Chris@208:    (if newsize == (size m) then
Chris@201:         m
Chris@261:     elif isSparse? m then
Chris@261:         // don't call makeSparse directly: want to discard
Chris@261:         // out-of-range cells
Chris@261:         newSparseMatrix (typeOf m) newsize (enumerateSparse 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@260:                 rowSlice m 0 newsize.rows
Chris@201:             elif (not rowm) and growcols < 0 then 
Chris@260:                 columnSlice m 0 newsize.columns
Chris@201:             elif growrows < 0 then 
Chris@260:                 rowSlice m 0 newsize.rows
Chris@201:             elif growcols < 0 then 
Chris@260:                 columnSlice m 0 newsize.columns
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@246:     density,
Chris@249:     nonZeroValues,
Chris@260:     at = at',
Chris@208:     getColumn,
Chris@208:     getRow,
Chris@208:     isRowMajor?,
Chris@234:     isSparse?,
Chris@208:     generate,
Chris@208:     constMatrix,
Chris@208:     randomMatrix,
Chris@208:     zeroMatrix,
Chris@208:     identityMatrix,
Chris@208:     zeroSizeMatrix,
Chris@208:     equal,
Chris@226:     equalUnder,
Chris@208:     transposed,
Chris@208:     flipped,
Chris@208:     toRowMajor,
Chris@208:     toColumnMajor,
Chris@238:     toSparse,
Chris@238:     toDense,
Chris@208:     scaled,
Chris@208:     resizedTo,
Chris@208:     asRows,
Chris@208:     asColumns,
Chris@208:     sum = sum',
Chris@229:     difference,
Chris@229:     abs = abs',
Chris@258:     filter,
Chris@208:     product,
Chris@208:     concat,
Chris@208:     rowSlice,
Chris@208:     columnSlice,
Chris@208:     newMatrix,
Chris@208:     newRowVector,
Chris@208:     newColumnVector,
Chris@261:     newSparseMatrix,
Chris@255:     enumerate
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@246:     density is matrix -> number,
Chris@249:     nonZeroValues is matrix -> number,
Chris@260:     at is matrix -> number -> number -> number,
Chris@214:     getColumn is number -> matrix -> vector,
Chris@214:     getRow is number -> matrix -> vector,
Chris@208:     isRowMajor? is matrix -> boolean,
Chris@234:     isSparse? 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@226:     equalUnder is (number -> number -> boolean) -> matrix -> matrix -> boolean,
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@243:     toSparse is matrix -> matrix,
Chris@238:     toDense is matrix -> matrix,
Chris@195:     scaled is number -> matrix -> matrix,
Chris@243:     thresholded 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@229:     difference is matrix -> matrix -> matrix,
Chris@229:     abs is matrix -> matrix,
Chris@258:     filter is (number -> boolean) -> matrix -> matrix,
Chris@195:     product is matrix -> matrix -> matrix,
Chris@195:     concat is (Horizontal () | Vertical ()) -> list<matrix> -> matrix,
Chris@260:     rowSlice is matrix -> number -> number -> matrix, 
Chris@260:     columnSlice is matrix -> number -> number -> matrix,
Chris@214:     newMatrix is (ColumnMajor () | RowMajor ()) -> list<vector> -> matrix, 
Chris@214:     newRowVector is vector -> matrix, 
Chris@214:     newColumnVector is vector -> matrix,
Chris@255:     newSparseMatrix is (ColumnMajor () | RowMajor ()) -> { .rows is number, .columns is number } -> list<{ .i is number, .j is number, .v is number }> -> matrix,
Chris@255:     enumerate is matrix -> list<{ .i is number, .j is number, .v is number }>
Chris@5: }
Chris@5: