diff yetilab/matrix/matrix.yeti @ 257:f00ab8baa6d7

Make sum, difference, scaled and abs use sparse operations and return sparse matrices if all inputs are sparse

--- a/yetilab/matrix/matrix.yeti	Tue May 21 22:36:39 2013 +0100
+++ b/yetilab/matrix/matrix.yeti	Tue May 21 22:37:28 2013 +0100
@@ -121,6 +121,12 @@
     _: vec.fromList (map do j: getAt i j m done [0..width m - 1]);
     esac;
 
+asRows m =
+    map do i: getRow i m done [0 .. (height m) - 1];
+
+asColumns m =
+    map do i: getColumn i m done [0 .. (width m) - 1];
+
 isRowMajor? m =
     case m of
     DenseRows _: true;
@@ -169,6 +175,7 @@
 swapij =
     map do { i, j, v }: { i = j, j = i, v } done;
 
+//!!! should use { row = , column = , value = } instead of i, j, v?
 enumerateSparse m =
    (enumerate { values, indices, pointers } =
         concat
@@ -349,30 +356,83 @@
 newColumnVector data = //!!! NB does not copy data
     DenseCols (array [data]);
 
-scaled factor m = //!!! v inefficient
-    generate do row col: factor * (getAt row col m) done (size m);
-
 thresholded threshold m = //!!! v inefficient; and should take a threshold function?
     generate do row col:
         v = getAt row col m; if (abs v) > threshold then v else 0 fi
     done (size m);
 
+denseLinearOp op m1 m2 =
+    if isRowMajor? m1 then
+        newMatrix (RowMajor ()) 
+           (map2 do c1 c2: op c1 c2 done (asRows m1) (asRows m2));
+    else
+        newMatrix (ColumnMajor ()) 
+           (map2 do c1 c2: op c1 c2 done (asColumns m1) (asColumns m2));
+    fi;
+
+sparseSumOrDifference op m1 m2 =
+   (h = [:];
+    for (enumerate m1) do { i, j, v }:
+        if not (i in h) then h[i] := [:] fi;
+        h[i][j] := v;
+    done;
+    for (enumerate m2) do { i, j, v }:
+        if not (i in h) then h[i] := [:] fi;
+        if j in h[i] then h[i][j] := op h[i][j] v;
+        else h[i][j] := op 0 v;
+        fi;
+    done;
+    entries = concat
+       (map do i:
+            kk = keys h[i];
+            map2 do j v: { i, j, v } done kk (map (at h[i]) kk)
+            done (keys h));
+    makeSparse
+        if isRowMajor? m1 then (RowMajor ()) else (ColumnMajor ()) fi
+        (size m1) 
+        entries);
+
 sum' m1 m2 =
     if (size m1) != (size m2)
     then failWith "Matrices are not the same size: \(size m1), \(size m2)";
+    elif isSparse? m1 and isSparse? m2 then
+        sparseSumOrDifference (+) m1 m2;
     else
-        generate do row col: getAt row col m1 + getAt row col m2 done (size m1);
+        denseLinearOp bf.add m1 m2;
     fi;
 
-difference m1 m2 = //!!! doc: m1 - m2, not m2 - m1
+difference m1 m2 =
     if (size m1) != (size m2)
     then failWith "Matrices are not the same size: \(size m1), \(size m2)";
+    elif isSparse? m1 and isSparse? m2 then
+        sparseSumOrDifference (-) m1 m2;
     else
-        generate do row col: getAt row col m1 - getAt row col m2 done (size m1);
+        denseLinearOp bf.subtract m1 m2;
+    fi;
+
+scaled factor m =
+    if isSparse? m then
+        makeSparse
+            if isRowMajor? m then (RowMajor ()) else (ColumnMajor ()) fi
+            (size m)
+            (map do { i, j, v }: { i, j, v = factor * v } done (enumerate m))
+    elif isRowMajor? m then
+        newMatrix (RowMajor ()) (map (bf.scaled factor) (asRows m));
+    else
+        newMatrix (ColumnMajor ()) (map (bf.scaled factor) (asColumns m));
     fi;
 
 abs' m =
-    generate do row col: abs (getAt row col m) done (size m);
+    if isSparse? m then
+        makeSparse
+            if isRowMajor? m then (RowMajor ()) else (ColumnMajor ()) fi
+            (size m)
+            (map do { i, j, v }: { i, j, v = abs v } done (enumerate m))
+    elif isRowMajor? m then
+        newMatrix (RowMajor ()) (map bf.abs (asRows m));
+    else
+        newMatrix (ColumnMajor ()) (map bf.abs (asColumns m));
+    fi;
 
 sparseProductLeft size m1 m2 =
    (e = enumerateSparse m1;
@@ -451,12 +511,6 @@
         fi;
     fi;
 
-asRows m =
-    map do i: getRow i m done [0 .. (height m) - 1];
-
-asColumns m =
-    map do i: getColumn i m done [0 .. (width m) - 1];
-
 concatAgainstGrain tagger getter counter mm =
    (n = counter (size (head mm));
     tagger (array