Mercurial > hg > may

changeset 257:f00ab8baa6d7

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make sum, difference, scaled and abs use sparse operations and return sparse matrices if all inputs are sparse
author Chris Cannam
date Tue, 21 May 2013 22:37:28 +0100
parents 97e97a2e6a4e
children f3b7b5d20f88
files yetilab/matrix/matrix.yeti yetilab/matrix/test/speedtest.yeti yetilab/matrix/test/test_matrix.yeti
diffstat 3 files changed, 162 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- 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
--- a/yetilab/matrix/test/speedtest.yeti	Tue May 21 22:36:39 2013 +0100
+++ b/yetilab/matrix/test/speedtest.yeti	Tue May 21 22:37:28 2013 +0100
@@ -83,7 +83,7 @@
    (print "                                 ";
     println "isSparse: \(mat.isSparse? m), density \(mat.density m)");
 
-println "\nM * M multiplies:\n";
+println "\nM * M multiplies (and a few sums):\n";
 
 sz = 500;
 
@@ -131,6 +131,34 @@
 
 println "";
 
+print "CMD + CMD... ";
+reportOn (time \(mat.sum cmd cmd));
+
+print "CMD + RMD... ";
+reportOn (time \(mat.sum cmd rmd));
+
+print "RMD + CMD... ";
+reportOn (time \(mat.sum rmd cmd));
+
+print "RMD + RMD... ";
+reportOn (time \(mat.sum rmd rmd));
+
+println "";
+
+print "CMS + CMS... ";
+reportOn (time \(mat.sum cms cms));
+
+print "CMS + RMS... ";
+reportOn (time \(mat.sum cms rms));
+
+print "RMS + CMS... ";
+reportOn (time \(mat.sum rms cms));
+
+print "RMS + RMS... ";
+reportOn (time \(mat.sum rms rms));
+
+println "";
+
 print "CMD * CMD... ";
 reportOn (time \(mat.product cmd cmd));
 
@@ -143,7 +171,7 @@
 print "RMD * RMD... ";
 reportOn (time \(mat.product rmd rmd));
 
-println "\nLarge sparse M * M multiplies:\n";
+println "\nLarge sparse M * M multiplies and adds:\n";
 
 sz = 5000000;
 nnz = 10000;
@@ -178,5 +206,18 @@
 print "RMS * RMS... ";
 reportOn (time \(mat.product rms rms));
 
+println "";
+
+print "CMS + CMS... ";
+reportOn (time \(mat.sum cms cms));
+
+print "CMS + RMS... ";
+reportOn (time \(mat.sum cms rms));
+
+print "RMS + CMS... ";
+reportOn (time \(mat.sum rms cms));
+
+print "RMS + RMS... ";
+reportOn (time \(mat.sum rms rms));
 
 ();
--- a/yetilab/matrix/test/test_matrix.yeti	Tue May 21 22:36:39 2013 +0100
+++ b/yetilab/matrix/test/test_matrix.yeti	Tue May 21 22:37:28 2013 +0100
@@ -215,6 +215,32 @@
     yrt
 ),
 
+"sparseSum-\(name)": \(
+    s = mat.newSparseMatrix (ColumnMajor ()) { rows = 2, columns = 3 } [
+        { i = 0, j = 0, v = 1 },
+        { i = 0, j = 2, v = 2 },
+        { i = 1, j = 1, v = 4 },
+    ];
+    t = mat.newSparseMatrix (ColumnMajor ()) { rows = 2, columns = 3 } [
+        { i = 0, j = 1, v = 7 },
+        { i = 1, j = 0, v = 5 },
+        { i = 1, j = 1, v = 6 },
+    ];
+    tot = mat.sum s t;
+    mat.isSparse? tot and
+        compareMatrices tot (mat.sum (mat.toDense s) t) and
+        compareMatrices tot (mat.sum (mat.toDense s) (mat.toDense t)) and
+        compareMatrices tot (mat.sum s (mat.toDense t)) and
+        compareMatrices tot 
+           (mat.newSparseMatrix (RowMajor ()) { rows = 2, columns = 3 } [
+               { i = 0, j = 0, v = 1 },
+               { i = 0, j = 1, v = 7 },
+               { i = 0, j = 2, v = 2 },
+               { i = 1, j = 0, v = 5 },
+               { i = 1, j = 1, v = 10 },
+            ])
+),
+
 "difference-\(name)": \(
     compareMatrices
        (mat.difference (constMatrix 2 { rows = 3, columns = 4 })
@@ -232,6 +258,32 @@
     yrt
 ),
 
+"sparseDifference-\(name)": \(
+    s = mat.newSparseMatrix (ColumnMajor ()) { rows = 2, columns = 3 } [
+        { i = 0, j = 0, v = 1 },
+        { i = 0, j = 2, v = 2 },
+        { i = 1, j = 1, v = 4 },
+    ];
+    t = mat.newSparseMatrix (ColumnMajor ()) { rows = 2, columns = 3 } [
+        { i = 0, j = 1, v = 7 },
+        { i = 1, j = 0, v = 5 },
+        { i = 1, j = 1, v = 6 },
+    ];
+    diff = mat.difference s t;
+    mat.isSparse? diff and
+        compareMatrices diff (mat.difference (mat.toDense s) t) and
+        compareMatrices diff (mat.difference (mat.toDense s) (mat.toDense t)) and
+        compareMatrices diff (mat.difference s (mat.toDense t)) and
+        compareMatrices diff 
+           (mat.newSparseMatrix (RowMajor ()) { rows = 2, columns = 3 } [
+               { i = 0, j = 0, v = 1 },
+               { i = 0, j = 1, v = -7 },
+               { i = 0, j = 2, v = 2 },
+               { i = 1, j = 0, v = -5 },
+               { i = 1, j = 1, v = -2 },
+            ])
+),
+
 "abs-\(name)": \(
     compareMatrices
        (mat.abs (newMatrix (ColumnMajor ()) [[-1,4],[2,-5],[-3,0]]))