changeset 250:6a141098c059

Merge from sparse branch

--- a/yetilab/matrix/matrix.yeti	Sun May 19 18:42:10 2013 +0100
+++ b/yetilab/matrix/matrix.yeti	Mon May 20 22:17:35 2013 +0100
@@ -19,68 +19,128 @@
 
 size m =
     case m of
-    RowM r:
+    DenseRows r:
         major = length r;
         { 
             rows = major, 
             columns = if major > 0 then vec.length r[0] else 0 fi,
         };
-    ColM c:
+    DenseCols c:
         major = length c;
         { 
             rows = if major > 0 then vec.length c[0] else 0 fi,
             columns = major, 
         };
+    SparseCSR { values, indices, pointers, extent }:
+        {
+            rows = (length pointers) - 1,
+            columns = extent
+        };
+    SparseCSC { values, indices, pointers, extent }:
+        {
+            rows = extent,
+            columns = (length pointers) - 1
+        };
     esac;
 
 width m = (size m).columns;
 height m = (size m).rows;
 
+nonZeroValues m =
+   (nz d =
+        sum
+           (map do v:
+                sum (map do n: if n == 0 then 0 else 1 fi done (vec.list v))
+                done d);
+    case m of 
+    DenseRows d: nz d;
+    DenseCols d: nz d;
+    SparseCSR d: vec.length d.values;
+    SparseCSC d: vec.length d.values;
+    esac);
+
+density m =
+   ({ rows, columns } = size m;
+    cells = rows * columns;
+    (nonZeroValues m) / cells);
+
+sparseSlice n d =
+   (start = d.pointers[n];
+    end = d.pointers[n+1];
+    { 
+        values = vec.slice d.values start end,
+        indices = slice d.indices start end,
+    });
+
+fromSlice n m d =
+   (slice = sparseSlice n d;
+    var v = 0;
+    for [0..length slice.indices - 1] do i:
+        if slice.indices[i] == m then
+            v := vec.at i slice.values;
+        fi
+    done;
+    v);
+
+filledSlice n d =
+   (slice = sparseSlice n d;
+    dslice = new double[d.extent];
+    for [0..length slice.indices - 1] do i:
+        dslice[slice.indices[i]] := vec.at i slice.values;
+    done;
+    vec.vector dslice);
+
 getAt row col m =
     case m of
-    RowM rows: r = rows[row]; vec.at col r;
-    ColM cols: c = cols[col]; vec.at row c;
+    DenseRows rows: r = rows[row]; vec.at col r;
+    DenseCols cols: c = cols[col]; vec.at row c;
+    SparseCSR data: fromSlice row col data;
+    SparseCSC data: fromSlice col row data;
     esac;
 
 getColumn j m =
     case m of
-    RowM rows: vec.fromList (map do i: getAt i j m done [0..length rows-1]);
-    ColM cols: cols[j];
+    DenseCols cols: cols[j];
+    SparseCSC data: filledSlice j data;
+    _: vec.fromList (map do i: getAt i j m done [0..height m - 1]);
     esac;
 
 getRow i m =
     case m of
-    RowM rows: rows[i];
-    ColM cols: vec.fromList (map do j: getAt i j m done [0..length cols-1]);
+    DenseRows rows: rows[i];
+    SparseCSR data: filledSlice i data; 
+    _: vec.fromList (map do j: getAt i j m done [0..width m - 1]);
     esac;
 
-/*
-setAt row col n m = //!!! dangerous, could modify copies -- should it be allowed?
-    case m of
-    RowM rows: r = rows[row]; (vec.data r)[col] := n;
-    ColM cols: c = cols[col]; (vec.data c)[row] := n;
-    esac;
-*/
-
 isRowMajor? m =
     case m of
-    RowM _: true;
-    ColM _: false;
+    DenseRows _: true;
+    DenseCols _: false;
+    SparseCSR _: true;
+    SparseCSC _: false;
     esac;
 
-newColMajorStorage { rows, columns } = 
+isSparse? m =
+    case m of
+    DenseRows _: false;
+    DenseCols _: false;
+    SparseCSR _: true;
+    SparseCSC _: true;
+    esac;
+
+newColumnMajorStorage { rows, columns } = 
     if rows < 1 then array []
     else array (map \(vec.zeros rows) [1..columns])
     fi;
 
 zeroMatrix { rows, columns } = 
-    ColM (newColMajorStorage { rows, columns });
+    DenseCols (newColumnMajorStorage { rows, columns });
 
 zeroMatrixWithTypeOf m { rows, columns } = 
     if isRowMajor? m then
-        RowM (newColMajorStorage { rows = columns, columns = rows });
+        DenseRows (newColumnMajorStorage { rows = columns, columns = rows });
     else
-        ColM (newColMajorStorage { rows, columns });
+        DenseCols (newColumnMajorStorage { rows, columns });
     fi;
 
 zeroSizeMatrix () = zeroMatrix { rows = 0, columns = 0 };
@@ -94,7 +154,90 @@
                 m[col][row] := f row col;
             done;
         done;
-        ColM (array (map vec.vector m))
+        DenseCols (array (map vec.vector m))
+    fi;
+
+enumerateSparse m =
+   (enumerate { values, indices, pointers } =
+        concat
+           (map do i:
+                start = pointers[i];
+                end = pointers[i+1];
+                map2 do j v: { i, j, v } done 
+                    (slice indices start end)
+                    (vec.list (vec.slice values start end))
+                done [0..length pointers - 2]);
+    case m of
+    SparseCSC d: 
+        map do { i, j, v }: { i = j, j = i, v } done (enumerate d);
+    SparseCSR d:
+        enumerate d;
+     _: [];
+    esac);
+
+makeSparse type size data =
+   (isRow = case type of RowMajor (): true; ColumnMajor (): false esac;
+    ordered = 
+        sortBy do a b:
+            if a.maj == b.maj then a.min < b.min else a.maj < b.maj fi
+        done
+           (map
+                if isRow then
+                    do { i, j, v }: { maj = i, min = j, v } done;
+                else
+                    do { i, j, v }: { maj = j, min = i, v } done;
+                fi
+                data);
+    tagger = if isRow then SparseCSR else SparseCSC fi;
+    majorSize = if isRow then size.rows else size.columns fi;
+    minorSize = if isRow then size.columns else size.rows fi;
+    majorPointers acc nn n i data =
+        if n < nn then 
+            case data of 
+            d::rest:
+                majorPointers (acc ++ (map \(i) [n..d-1])) nn d (i+1) rest;
+             _: 
+                majorPointers (acc ++ [i]) nn (n+1) i [];
+            esac;
+        else
+            acc
+        fi;
+    tagger {
+        values = vec.fromList (map (.v) ordered),
+        indices = array (map (.min) ordered),
+        pointers = array (majorPointers [0] majorSize 0 0 (map (.maj) ordered)),
+        extent = minorSize,
+    });
+
+toSparse m =
+    if isSparse? m then m
+    else
+        { rows, columns } = size m;
+        enumerate m ii jj =
+            case ii of
+            i::irest:
+                case jj of
+                j::rest:
+                    v = getAt i j m;
+                    if v != 0 then { i, j, v } :. \(enumerate m ii rest)
+                    else enumerate m ii rest
+                    fi;
+                 _: enumerate m irest [0..columns-1];
+                esac;
+             _: [];
+            esac;
+        makeSparse 
+            if isRowMajor? m then RowMajor () else ColumnMajor () fi
+               (size m)
+               (enumerate m [0..rows-1] [0..columns-1]);
+    fi;
+
+toDense m =
+    if not (isSparse? m) then m
+    elif isRowMajor? m then
+        DenseRows (array (map do row: getRow row m done [0..height m - 1]));
+    else
+        DenseCols (array (map do col: getColumn col m done [0..width m - 1]));
     fi;
 
 constMatrix n = generate do row col: n done;
@@ -103,17 +246,27 @@
 
 transposed m =
     case m of
-    RowM d: ColM d;
-    ColM d: RowM d;
+    DenseRows d: DenseCols d;
+    DenseCols d: DenseRows d;
+    SparseCSR d: SparseCSC d;
+    SparseCSC d: SparseCSR d;
     esac;
 
 flipped m =
-    if isRowMajor? m then
-        generate do row col: getAt row col m done (size m);
+    if isSparse? m then
+        if isRowMajor? m then
+            makeSparse (ColumnMajor ()) (size m) (enumerateSparse m)
+        else
+            makeSparse (RowMajor ()) (size m) (enumerateSparse m)
+        fi
     else
-        transposed
-           (generate do row col: getAt col row m done
-            { rows = (width m), columns = (height m) });
+        if isRowMajor? m then
+            generate do row col: getAt row col m done (size m);
+        else
+            transposed
+               (generate do row col: getAt col row m done
+                { rows = (width m), columns = (height m) });
+        fi
     fi;
 
 toRowMajor m =
@@ -122,53 +275,72 @@
 toColumnMajor m =
     if not isRowMajor? m then m else flipped m fi;
 
-equal' vecComparator m1 m2 =
+equal'' comparator vecComparator m1 m2 =
+    // Prerequisite: m1 and m2 have same sparse-p and storage order
+   (compareVecLists vv1 vv2 = all id (map2 vecComparator vv1 vv2);
+    compareSparse d1 d2 =
+        d1.extent == d2.extent and
+        vecComparator d1.values d2.values and
+        d1.indices == d2.indices and
+        d1.pointers == d2.pointers;
+    case m1 of
+    DenseRows d1:
+        case m2 of DenseRows d2: compareVecLists d1 d2; _: false; esac;
+    DenseCols d1:
+        case m2 of DenseCols d2: compareVecLists d1 d2; _: false; esac;
+    SparseCSR d1:
+        case m2 of SparseCSR d2: compareSparse d1 d2; _: false; esac;
+    SparseCSC d1:
+        case m2 of SparseCSC d2: compareSparse d1 d2; _: false; esac;
+    esac);
+
+equal' comparator vecComparator m1 m2 =
     if size m1 != size m2 then 
         false
     elif isRowMajor? m1 != isRowMajor? m2 then
-        equal' vecComparator (flipped m1) m2;
+        equal' comparator vecComparator (flipped m1) m2;
+    elif isSparse? m1 != isSparse? m2 then
+        if isSparse? m1 then
+            equal' comparator vecComparator m1 (toSparse m2)
+        else
+            equal' comparator vecComparator (toSparse m1) m2
+        fi
     else
-        compare d1 d2 = all id (map2 vecComparator d1 d2);
-        case m1 of
-        RowM d1: case m2 of RowM d2: compare d1 d2; _: false; esac;
-        ColM d1: case m2 of ColM d2: compare d1 d2; _: false; esac;
-        esac
+        equal'' comparator vecComparator m1 m2
     fi;
 
 // Compare matrices using the given comparator for individual cells.
 // Note that matrices with different storage order but the same
 // contents are equal, although comparing them is slow.
+//!!! Document the fact that sparse matrices can only be equal if they
+// have the same set of non-zero cells (regardless of comparator used)
 equalUnder comparator =
-    equal' (vec.equalUnder comparator);
+    equal' comparator (vec.equalUnder comparator);
 
 equal =
-    equal' vec.equal;
-
-/*!!! not needed now it's immutable?
-copyOf m =
-   (copyOfData d = (array (map vec.copyOf d));
-    case m of
-    RowM d: RowM (copyOfData d);
-    ColM d: ColM (copyOfData d);
-    esac);
-*/
+    equal' (==) vec.equal;
 
 newMatrix type data = //!!! NB does not copy data
-   (tagger = case type of RowMajor (): RowM; ColumnMajor (): ColM esac;
+   (tagger = case type of RowMajor (): DenseRows; ColumnMajor (): DenseCols esac;
     if empty? data or vec.empty? (head data)
     then zeroSizeMatrix ()
     else tagger (array data)
     fi);
 
 newRowVector data = //!!! NB does not copy data
-    RowM (array [data]);
+    DenseRows (array [data]);
 
 newColumnVector data = //!!! NB does not copy data
-    ColM (array [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);
+
 sum' m1 m2 =
     if (size m1) != (size m2)
     then failWith "Matrices are not the same size: \(size m1), \(size m2)";
@@ -186,13 +358,50 @@
 abs' m =
     generate do row col: abs (getAt row col m) done (size m);
 
+sparseProductLeft size m1 m2 =
+   (e = enumerateSparse m1;
+    data = array (map \(new double[size.rows]) [1..size.columns]);
+    for [0..size.columns - 1] do j':
+        c = getColumn j' m2;
+        for e do { v, i, j }:
+            data[j'][i] := data[j'][i] + v * (vec.at j c);
+        done;
+    done;
+    DenseCols (array (map vec.vector (list data))));
+
+sparseProductRight size m1 m2 =
+   (e = enumerateSparse m2;
+    data = array (map \(new double[size.columns]) [1..size.rows]);
+    for [0..size.rows - 1] do i':
+        r = getRow i' m1;
+        for e do { v, i, j }:
+            data[i'][j] := data[i'][j] + v * (vec.at i r);
+        done;
+    done;
+    DenseRows (array (map vec.vector (list data))));
+
+denseProduct size m1 m2 =
+   (data = array (map \(new double[size.rows]) [1..size.columns]);
+    for [0..size.rows - 1] do i:
+        row = getRow i m1;
+        for [0..size.columns - 1] do j:
+            data[j][i] := bf.sum (bf.multiply row (getColumn j m2));
+        done;
+    done;
+    DenseCols (array (map vec.vector (list data))));
+
 product m1 m2 =
     if (size m1).columns != (size m2).rows
-    then failWith "Matrix dimensions incompatible: \(size m1), \(size m2) (\((size m1).columns != (size m2).rows)";
-    else
-        generate do row col:
-            bf.sum (bf.multiply (getRow row m1) (getColumn col m2))
-        done { rows = (size m1).rows, columns = (size m2).columns }
+    then failWith "Matrix dimensions incompatible: \(size m1), \(size m2) (\((size m1).columns) != \((size m2).rows))";
+    else 
+        size = { rows = (size m1).rows, columns = (size m2).columns };
+        if isSparse? m1 then
+            sparseProductLeft size m1 m2
+        elif isSparse? m2 then
+            sparseProductRight size m1 m2
+        else
+            denseProduct size m1 m2
+        fi;
     fi;
 
 asRows m =
@@ -235,30 +444,31 @@
         // vertical, col-major: against grain with cols
         case direction of
         Horizontal ():
-            if row then concatAgainstGrain RowM getRow (.rows) mm;
-            else concatWithGrain ColM getColumn (.columns) mm;
+            if row then concatAgainstGrain DenseRows getRow (.rows) mm;
+            else concatWithGrain DenseCols getColumn (.columns) mm;
             fi;
         Vertical ():
-            if row then concatWithGrain RowM getRow (.rows) mm;
-            else concatAgainstGrain ColM getColumn (.columns) mm;
+            if row then concatWithGrain DenseRows getRow (.rows) mm;
+            else concatAgainstGrain DenseCols getColumn (.columns) mm;
             fi;
         esac;
     [single]: single;
     _: zeroSizeMatrix ();
     esac;
 
+//!!! inconsistent with std.slice which has start..end not start+count (see also vec.slice/rangeOf)
 rowSlice start count m = //!!! doc: storage order same as input
     if isRowMajor? m then
-        RowM (array (map ((flip getRow) m) [start .. start + count - 1]))
+        DenseRows (array (map ((flip getRow) m) [start .. start + count - 1]))
     else 
-        ColM (array (map (vec.rangeOf start count) (asColumns m)))
+        DenseCols (array (map (vec.rangeOf start count) (asColumns m)))
     fi;
 
 columnSlice start count m = //!!! doc: storage order same as input
     if not isRowMajor? m then
-        ColM (array (map ((flip getColumn) m) [start .. start + count - 1]))
+        DenseCols (array (map ((flip getColumn) m) [start .. start + count - 1]))
     else 
-        RowM (array (map (vec.rangeOf start count) (asRows m)))
+        DenseRows (array (map (vec.rangeOf start count) (asRows m)))
     fi;
 
 resizedTo newsize m =
@@ -294,11 +504,13 @@
     size,
     width,
     height,
+    density,
+    nonZeroValues,
     getAt,
     getColumn,
     getRow,
-//    setAt,
     isRowMajor?,
+    isSparse?,
     generate,
     constMatrix,
     randomMatrix,
@@ -307,12 +519,14 @@
     zeroSizeMatrix,
     equal,
     equalUnder,
-//    copyOf,
     transposed,
     flipped,
     toRowMajor,
     toColumnMajor,
+    toSparse,
+    toDense,
     scaled,
+    thresholded,
     resizedTo,
     asRows,
     asColumns,
@@ -326,6 +540,7 @@
     newMatrix,
     newRowVector,
     newColumnVector,
+    newSparseMatrix = makeSparse
 }
 as
 {
@@ -334,11 +549,13 @@
     size is matrix -> { .rows is number, .columns is number },
     width is matrix -> number,
     height is matrix -> number,
+    density is matrix -> number,
+    nonZeroValues is matrix -> number,
     getAt is number -> number -> matrix -> number,
     getColumn is number -> matrix -> vector,
     getRow is number -> matrix -> vector,
-//    setAt is number -> number -> number -> matrix -> (), //!!! lose?
     isRowMajor? is matrix -> boolean,
+    isSparse? is matrix -> boolean,
     generate is (number -> number -> number) -> { .rows is number, .columns is number } -> matrix,
     constMatrix is number -> { .rows is number, .columns is number } -> matrix,
     randomMatrix is { .rows is number, .columns is number } -> matrix,
@@ -347,12 +564,14 @@
     zeroSizeMatrix is () -> matrix,
     equal is matrix -> matrix -> boolean,
     equalUnder is (number -> number -> boolean) -> matrix -> matrix -> boolean,
-//    copyOf is matrix -> matrix,
     transposed is matrix -> matrix,
     flipped is matrix -> matrix, 
     toRowMajor is matrix -> matrix, 
     toColumnMajor is matrix -> matrix,
+    toSparse is matrix -> matrix,
+    toDense is matrix -> matrix,
     scaled is number -> matrix -> matrix,
+    thresholded is number -> matrix -> matrix,
     resizedTo is { .rows is number, .columns is number } -> matrix -> matrix,
     asRows is matrix -> list<vector>, 
     asColumns is matrix -> list<vector>,
@@ -366,5 +585,6 @@
     newMatrix is (ColumnMajor () | RowMajor ()) -> list<vector> -> matrix, 
     newRowVector is vector -> matrix, 
     newColumnVector is vector -> matrix,
+    newSparseMatrix is (ColumnMajor () | RowMajor ()) -> { .rows is number, .columns is number } -> list<{ .i is number, .j is number, .v is number }> -> matrix
 }
 

--- a/yetilab/matrix/matrixtype.yeti	Sun May 19 18:42:10 2013 +0100
+++ b/yetilab/matrix/matrixtype.yeti	Mon May 20 22:17:35 2013 +0100
@@ -3,7 +3,21 @@
 
 load yetilab.vector.vectortype;
 
-typedef opaque matrix = RowM. array<vector> | ColM. array<vector>;
+typedef opaque matrix =
+    DenseRows. array<vector> | // array of rows
+    DenseCols. array<vector> | // array of columns
+    SparseCSR. {
+        .values is vector,
+        .indices is array<number>, // column index of each value
+        .pointers is array<number>, // offset of first value in each row
+        .extent is number // max possible index + 1, i.e. number of columns
+        } |
+    SparseCSC. {
+        .values is vector,
+        .indices is array<number>, // row index of each value
+        .pointers is array<number>, // offset of first value in each column
+        .extent is number // max pointers index + 1, i.e. number of rows
+        };
 
 ();
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/yetilab/matrix/test/speedtest.yeti	Mon May 20 22:17:35 2013 +0100
@@ -0,0 +1,142 @@
+
+program yetilab.matrix.test.speedtest;
+
+mat = load yetilab.matrix.matrix;
+vec = load yetilab.vector.vector;
+
+{ compare, compareUsing } = load yetilab.test.test;
+
+compareMatrices = compareUsing mat.equal;
+
+time f =
+   (start = System#currentTimeMillis();
+    result = f ();
+    end = System#currentTimeMillis();
+    println " \(end-start)ms";
+    result);
+
+makeMatrices sz sparsity =
+   (print "Making \(sz) * \(sz) random matrix...";
+    m = time \(mat.randomMatrix { rows = sz, columns = sz });
+    makeSparse () = 
+       (print "Making \(sparsity * 100)% sparse version (as dense matrix)...";
+        t = time \(mat.thresholded sparsity m);
+        println "Reported density: \(mat.density t) (non-zero values: \(mat.nonZeroValues t))";
+        print "Converting to sparse matrix...";
+        s = time \(mat.toSparse t);
+        println "Reported density: \(mat.density s) (non-zero values: \(mat.nonZeroValues s))";
+        s);
+    s = makeSparse ();
+    println "Making types:";
+    print "Col-major dense...";
+    cmd = time \(mat.toColumnMajor m);
+    print "Row-major dense...";
+    rmd = time \(mat.toRowMajor m);
+    print "Col-major sparse...";
+    cms = time \(mat.toColumnMajor s);
+    print "Row-major sparse...";
+    rms = time \(mat.toRowMajor s);
+    println "";
+    { cmd, rmd, cms, rms });
+
+println "\nR * M multiplies:\n";
+
+sz = 2000;
+sparsity = 0.98;
+
+{ cmd, rmd, cms, rms } = makeMatrices sz sparsity;
+
+row = mat.newRowVector (vec.fromList (map \(Math#random()) [1..sz]));
+col = mat.newColumnVector (vec.fromList (map \(Math#random()) [1..sz]));
+
+print "R * CMD... ";
+a = (time \(mat.product row cmd));
+
+print "R * RMD... ";
+b = (time \(mat.product row rmd));
+
+print "R * CMS... ";
+c = (time \(mat.product row cms));
+
+print "R * RMS... ";
+d = (time \(mat.product row rms));
+
+println "\nChecking results: \(compareMatrices a b) \(compareMatrices c d)";
+
+println "\nM * C multiplies:\n";
+
+print "CMD * C... ";
+a = (time \(mat.product cmd col));
+
+print "RMD * C... ";
+b = (time \(mat.product rmd col));
+
+print "CMS * C... ";
+c = (time \(mat.product cms col));
+
+print "RMS * C... ";
+d = (time \(mat.product rms col));
+
+println "\nChecking results: \(compareMatrices a b) \(compareMatrices c d)";
+
+println "\nM * M multiplies:\n";
+
+sz = 500;
+
+{ cmd, rmd, cms, rms } = makeMatrices sz sparsity;
+
+print "CMS * CMD... ";
+\() (time \(mat.product cms cmd));
+
+print "CMS * RMD... ";
+\() (time \(mat.product cms rmd));
+
+print "RMS * CMD... ";
+\() (time \(mat.product rms cmd));
+
+print "RMS * RMD... ";
+\() (time \(mat.product rms rmd));
+
+println "";
+
+print "CMD * CMS... ";
+\() (time \(mat.product cmd cms));
+
+print "CMD * RMS... ";
+\() (time \(mat.product cmd rms));
+
+print "RMD * CMS... ";
+\() (time \(mat.product rmd cms));
+
+print "RMD * RMS... ";
+\() (time \(mat.product rmd rms));
+
+println "";
+
+print "CMS * CMS... ";
+\() (time \(mat.product cms cms));
+
+print "CMS * RMS... ";
+\() (time \(mat.product cms rms));
+
+print "RMS * CMS... ";
+\() (time \(mat.product rms cms));
+
+print "RMS * RMS... ";
+\() (time \(mat.product rms rms));
+
+println "";
+
+print "CMD * CMD... ";
+\() (time \(mat.product cmd cmd));
+
+print "CMD * RMD... ";
+\() (time \(mat.product cmd rmd));
+
+print "RMD * CMD... ";
+\() (time \(mat.product rmd cmd));
+
+print "RMD * RMD... ";
+\() (time \(mat.product rmd rmd));
+
+();

--- a/yetilab/matrix/test/test_matrix.yeti	Sun May 19 18:42:10 2013 +0100
+++ b/yetilab/matrix/test/test_matrix.yeti	Mon May 20 22:17:35 2013 +0100
@@ -90,11 +90,11 @@
 ),
 
 "equal-\(name)": \(
-    m = newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]];
+    m = newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]];
     n = m;
-    p = newMatrix (RowMajor ()) [[1,2,3],[4,5,6]];
-    q = newMatrix (ColumnMajor ()) [[1,2,3],[4,5,6]];
-    r = newMatrix (ColumnMajor ()) [[1,4],[2,5]];
+    p = newMatrix (RowMajor ()) [[1,0,3],[4,5,6]];
+    q = newMatrix (ColumnMajor ()) [[1,0,3],[4,5,6]];
+    r = newMatrix (ColumnMajor ()) [[1,4],[0,5]];
     compareMatrices m n and
         compareMatrices m p and
         compareMatrices n p and
@@ -106,8 +106,8 @@
     p = newMatrix (ColumnMajor ()) [[1,2,3],[4,5,6]];
     q = newMatrix (ColumnMajor ()) [[1,2,3],[4,5,6]];
     r = newMatrix (ColumnMajor ()) [[4,3,1],[3,1,2]];
-    s = newMatrix (ColumnMajor ()) [[1,2,5],[6,7,8]];
-    t = newMatrix (ColumnMajor ()) [[1,2,5],[6,7,9]];
+    s = newMatrix (ColumnMajor ()) [[1,4,5],[6,7,8]];
+    t = newMatrix (ColumnMajor ()) [[1,4,5],[6,7,9]];
     mat.equalUnder (==) p p and
         mat.equalUnder (==) p q and
         mat.equalUnder (!=) p r and
@@ -123,30 +123,7 @@
            (map do col: mat.getAt row col m == generator row col done [0..2])
             done [0..1])
 ),
-/*!!!
-"setAt-\(name)": \(
-    generator row col = row * 10 + col;
-    m = generate generator { rows = 2, columns = 3 };
-    mat.setAt 1 2 16 m;
-    compare (mat.getAt 1 2 m) 16 and
-        compare (mat.getAt 1 1 m) 11 and
-        compare (mat.getAt 0 2 m) 2
-),
 
-"copyOfEqual-\(name)": \(
-    m = constMatrix 2 { rows = 3, columns = 4 };
-    m'' = mat.copyOf m;
-    compareMatrices m'' m
-),
-
-"copyOfAlias-\(name)": \(
-    m = constMatrix 2 { rows = 3, columns = 4 };
-    m' = m;
-    m'' = mat.copyOf m;
-    mat.setAt 0 0 6 m;
-    compareMatrices m' m and not mat.equal m m'';
-),
-*/
 "transposedEmpty-\(name)": \(
     compare (mat.size (mat.transposed (constMatrix 2 { rows = 0, columns = 0 }))) { columns = 0, rows = 0 } and
         compare (mat.size (mat.transposed (constMatrix 2 { rows = 0, columns = 4 }))) { columns = 0, rows = 0 } and
@@ -175,14 +152,14 @@
 ),
 
 "flipped-\(name)": \(
-    m = newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]];
+    m = newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]];
     m' = mat.flipped m;
-    m'' = newMatrix (RowMajor ()) [[1,2,3],[4,5,6]];
+    m'' = newMatrix (RowMajor ()) [[1,0,3],[4,5,6]];
     compareMatrices m m' and compareMatrices m m'' and compareMatrices m' m'';
 ),
 
 "flipped-back-\(name)": \(
-    m = newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]];
+    m = newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]];
     compareMatrices m (mat.flipped (mat.flipped m));
 ),
 
@@ -192,18 +169,18 @@
 ),
 
 "toRowMajor-\(name)": \(
-    m = newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]];
+    m = newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]];
     m' = mat.toRowMajor m;
-    m'' = newMatrix (RowMajor ()) [[1,2,3],[4,5,6]];
+    m'' = newMatrix (RowMajor ()) [[1,0,3],[4,5,6]];
     m''' = mat.toRowMajor m'';
     compareMatrices m m' and compareMatrices m m'' and compareMatrices m' m''
         and compareMatrices m m''';
 ),
 
 "toColumnMajor-\(name)": \(
-    m = newMatrix (RowMajor ()) [[1,4],[2,5],[3,6]];
+    m = newMatrix (RowMajor ()) [[1,4],[0,5],[3,6]];
     m' = mat.toColumnMajor m;
-    m'' = newMatrix (ColumnMajor ()) [[1,2,3],[4,5,6]];
+    m'' = newMatrix (ColumnMajor ()) [[1,0,3],[4,5,6]];
     m''' = mat.toColumnMajor m'';
     compareMatrices m m' and compareMatrices m m'' and compareMatrices m' m''
         and compareMatrices m m''';
@@ -310,29 +287,29 @@
 "asRows-\(name)": \(
     compare
        (map vec.list
-           (mat.asRows (newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]])))
-        [[1,2,3],[4,5,6]];
+           (mat.asRows (newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]])))
+        [[1,0,3],[4,5,6]];
 ),
 
 "asColumns-\(name)": \(
     compare
        (map vec.list
-           (mat.asColumns (newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]])))
-        [[1,4],[2,5],[3,6]];
+           (mat.asColumns (newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]])))
+        [[1,4],[0,5],[3,6]];
 ),
 
 "concat-horiz-\(name)": \(
     compareMatrices
        (mat.concat (Horizontal ()) 
-          [(newMatrix (ColumnMajor ()) [[1,4],[2,5]]),
+          [(newMatrix (ColumnMajor ()) [[1,4],[0,5]]),
            (newMatrix (RowMajor ()) [[3],[6]])])
-       (newMatrix (ColumnMajor ()) [[1,4],[2,5],[3,6]])
+       (newMatrix (ColumnMajor ()) [[1,4],[0,5],[3,6]])
 ),
 
 "concatFail-horiz-\(name)": \(
     try
         \() (mat.concat (Horizontal ()) 
-          [(newMatrix (ColumnMajor ()) [[1,4],[2,5]]),
+          [(newMatrix (ColumnMajor ()) [[1,4],[0,5]]),
            (newMatrix (ColumnMajor ()) [[3],[6]])]);
         false
     catch FailureException e:
@@ -343,15 +320,15 @@
 "concat-vert-\(name)": \(
     compareMatrices
        (mat.concat (Vertical ()) 
-          [(newMatrix (ColumnMajor ()) [[1,4],[2,5]]),
+          [(newMatrix (ColumnMajor ()) [[1,4],[0,5]]),
            (newMatrix (RowMajor ()) [[3,6]])])
-       (newMatrix (ColumnMajor ()) [[1,4,3],[2,5,6]])
+       (newMatrix (ColumnMajor ()) [[1,4,3],[0,5,6]])
 ),
 
 "concatFail-vert-\(name)": \(
     try
         \() (mat.concat (Vertical ()) 
-          [(newMatrix (ColumnMajor ()) [[1,4],[2,5]]),
+          [(newMatrix (ColumnMajor ()) [[1,4],[0,5]]),
            (newMatrix (RowMajor ()) [[3],[6]])]);
         false
     catch FailureException e:
@@ -361,24 +338,64 @@
 
 "rowSlice-\(name)": \(
     compareMatrices
-       (mat.rowSlice 1 2 (newMatrix (RowMajor ()) [[1,2],[3,4],[5,6],[7,8]]))
-       (newMatrix (RowMajor ()) [[3,4],[5,6]])
+       (mat.rowSlice 1 2 (newMatrix (RowMajor ()) [[1,0],[3,4],[0,6],[7,8]]))
+       (newMatrix (RowMajor ()) [[3,4],[0,6]])
 ),
 
 "columnSlice-\(name)": \(
     compareMatrices
-       (mat.columnSlice 1 2 (newMatrix (RowMajor ()) [[1,2,3,4],[5,6,7,8]]))
-       (newMatrix (RowMajor ()) [[2,3],[6,7]])
+       (mat.columnSlice 1 2 (newMatrix (RowMajor ()) [[1,0,3,4],[0,6,7,8]]))
+       (newMatrix (RowMajor ()) [[0,3],[6,7]])
+),
+
+"density-\(name)": \(
+    compare (mat.density (newMatrix (ColumnMajor ()) [[1,2,0],[0,5,0]])) (3/6) and
+        compare (mat.density (newMatrix (ColumnMajor ()) [[1,2,3],[4,5,6]])) (6/6) and
+        compare (mat.density (newMatrix (ColumnMajor ()) [[0,0,0],[0,0,0]])) 0
+),
+
+"nonZeroValues-\(name)": \(
+    compare (mat.nonZeroValues (newMatrix (ColumnMajor ()) [[1,2,0],[0,5,0]])) 3 and
+        compare (mat.nonZeroValues (newMatrix (ColumnMajor ()) [[1,2,3],[4,5,6]])) 6 and
+        compare (mat.nonZeroValues (newMatrix (ColumnMajor ()) [[0,0,0],[0,0,0]])) 0
+),
+
+"toSparse-\(name)": \(
+    m = newMatrix (ColumnMajor ()) [[1,2,0],[-1,-4,6],[0,0,3]];
+    compareMatrices (mat.toSparse m) m and
+        compareMatrices (mat.toDense (mat.toSparse m)) m and
+        compare (mat.density (mat.toSparse m)) (6/9)
+),
+
+"toDense-\(name)": \(
+    m = newMatrix (ColumnMajor ()) [[1,2,0],[-1,-4,6],[0,0,3]];
+    compareMatrices (mat.toDense m) m and
+        compareMatrices (mat.toSparse (mat.toDense m)) m
+),
+
+"thresholded-\(name)": \(
+    m = newMatrix (ColumnMajor ()) [[1,2,0],[-1,-4,6],[0,0,3]];
+    compareMatrices
+       (mat.thresholded 2 m)
+       (newMatrix (ColumnMajor ()) [[0,0,0],[0,-4,6],[0,0,3]]) and
+        compare (mat.density (mat.thresholded 2 m)) (3/9)
 ),
 
 ]);
 
-colhash = makeTests "column-major" id;
-rowhash = makeTests "row-major" mat.flipped;
+colhash = makeTests "column-dense" id;
+rowhash = makeTests "row-dense" mat.flipped;
+sparsecolhash = makeTests "column-sparse" mat.toSparse;
+
+// there are two possible orders for constructing a sparse row-major
+// matrix from a dense col-major one, so test them both:
+sparserowhash1 = makeTests "row-sparse-a" (mat.toSparse . mat.flipped);
+sparserowhash2 = makeTests "row-sparse-b" (mat.flipped . mat.toSparse);
 
 all = [:];
-for (keys colhash) do k: all[k] := colhash[k] done;
-for (keys rowhash) do k: all[k] := rowhash[k] done;
+for [ colhash, rowhash, sparsecolhash, sparserowhash1, sparserowhash2 ] do h:
+    for (keys h) do k: all[k] := h[k] done;
+done;
 
 all is hash<string, () -> boolean>;
 

--- a/yetilab/vector/vector.yeti	Sun May 19 18:42:10 2013 +0100
+++ b/yetilab/vector/vector.yeti	Mon May 20 22:17:35 2013 +0100
@@ -72,6 +72,9 @@
 rangeOf start len v is number -> number -> ~double[] -> ~double[] =
     Arrays#copyOfRange(v, start, start + len);
 
+slice v start end is ~double[] -> number -> number -> ~double[] =
+    rangeOf start (end - start) v;
+
 resizedTo n v is number -> ~double[] -> ~double[] =
     Arrays#copyOf(v, n);
 
@@ -103,6 +106,7 @@
     equal,
     equalUnder,
     rangeOf,
+    slice,
     resizedTo,
     concat,
 } as {
@@ -122,6 +126,7 @@
     equal is vector -> vector -> boolean,
     equalUnder is (number -> number -> boolean) -> vector -> vector -> boolean,
     rangeOf is number -> number -> vector -> vector, //!!! not well-named now vector arg is at the end
+    slice is vector -> number -> number -> vector, //!!! duplication with rangeOf (std module function on arrays is called slice though)
     resizedTo is number -> vector -> vector,
     concat is list?<vector> -> vector,
 }