# HG changeset patch
# User Chris Cannam
# Date 1395418951 0
# Node ID bdee7ac8cbcdf4630ce643fefe3553b9b129be41
# Parent  acc244cab1b7b83e5d47647e6416894d54d41699
Provisionally add ivector, use in sparse matrices

diff -r acc244cab1b7 -r bdee7ac8cbcd src/may/matrix.yeti
--- a/src/may/matrix.yeti	Fri Mar 21 10:33:04 2014 +0000
+++ b/src/may/matrix.yeti	Fri Mar 21 16:22:31 2014 +0000
@@ -41,6 +41,7 @@
 module may.matrix;
 
 vec = load may.vector;
+ivec = load may.vector.ivector;
 
 load yeti.experimental.json;
 
@@ -51,14 +52,14 @@
         DenseCols array<vec.vector_t> | // array of columns
         SparseCSR {
             values is vec.vector_t,
-            indices is array<number>, // column index of each value
-            pointers is array<number>, // offset of first value in each row
+            indices is ivec.ivector_t, // column index of each value
+            pointers is ivec.ivector_t, // offset of first value in each row
             extent is number // max possible index + 1, i.e. number of columns
         } |
         SparseCSC {
             values is vec.vector_t,
-            indices is array<number>, // row index of each value
-            pointers is array<number>, // offset of first value in each column
+            indices is ivec.ivector_t, // row index of each value
+            pointers is ivec.ivector_t, // offset of first value in each column
             extent is number // max pointers index + 1, i.e. number of rows
         }
 };
@@ -86,17 +87,17 @@
     (nonZeroValues m) / cells);
 
 sparseSlice n d =
-   (start = d.pointers[n];
-    end = d.pointers[n+1];
+   (start = ivec.at d.pointers n;
+    end = ivec.at d.pointers (n+1);
     { 
         values = vec.slice d.values start end,
-        indices = slice d.indices start end,
+        indices = ivec.slice d.indices start end,
     });
 
 nonEmptySlices d =
    (ne = array [];
-    for [0..length d.pointers - 2] do i:
-        if d.pointers[i] != d.pointers[i+1] then
+    for [0..ivec.length d.pointers - 2] do i:
+        if ivec.at d.pointers i != ivec.at d.pointers (i+1) then
             push ne i
         fi
     done;
@@ -105,8 +106,8 @@
 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
+    for [0..ivec.length slice.indices - 1] do i:
+        if ivec.at slice.indices i == m then
             v := vec.at slice.values i;
         fi
     done;
@@ -115,8 +116,8 @@
 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 slice.values i;
+    for [0..ivec.length slice.indices - 1] do i:
+        dslice[ivec.at slice.indices i] := vec.at slice.values i;
     done;
     vec.vector dslice);
 
@@ -254,12 +255,12 @@
    (enumerate { values, indices, pointers } =
         concat
            (map do i:
-                start = pointers[i];
-                end = pointers[i+1];
+                start = ivec.at pointers i;
+                end = ivec.at pointers (i+1);
                 map2 do j v: { i, j, v } done 
-                    (slice indices start end)
+                    (ivec.list (ivec.slice indices start end))
                     (vec.list (vec.slice values start end))
-                done [0..length pointers - 2]);
+                done [0..ivec.length pointers - 2]);
     case m.data of
     SparseCSC d: swapij (enumerate d);
     SparseCSR d: enumerate d;
@@ -320,8 +321,8 @@
         size,
         data = tagger {
             values = vec.fromList (map (.v) ordered),
-            indices = array (map (.min) ordered),
-            pointers,
+            indices = ivec.fromList (map (.min) ordered),
+            pointers = ivec.fromList (list pointers),
             extent = minorSize,
         }
     });
@@ -401,8 +402,8 @@
     compareSparse d1 d2 =
         d1.extent == d2.extent and
         vecComparator d1.values d2.values and
-        d1.indices == d2.indices and
-        d1.pointers == d2.pointers;
+        ivec.equal d1.indices d2.indices and
+        ivec.equal d1.pointers d2.pointers;
     case m1.data of
     DenseRows d1:
         case m2.data of DenseRows d2: compareVecLists d1 d2; _: false; esac;
@@ -577,10 +578,10 @@
     for [0..size.columns - 1] do j':
         c = getColumn j' m2;
         var p = 0;
-        for [0..length indices - 1] do ix:
-            ix == pointers[p+1] loop (p := p + 1);
-            i = if rows then p else indices[ix] fi;
-            j = if rows then indices[ix] else p fi;
+        for [0..ivec.length indices - 1] do ix:
+            ix == ivec.at pointers (p+1) loop (p := p + 1);
+            i = if rows then p else ivec.at indices ix fi;
+            j = if rows then ivec.at indices ix else p fi;
             data[j'][i] := data[j'][i] + (vec.at values ix) * (vec.at c j);
         done;
     done;
@@ -591,17 +592,17 @@
         case m2.data of
         SparseCSR d: d;
         SparseCSC d: d;
-        _: failWith "sparseProductLeft called for non-sparse m1";
+        _: failWith "sparseProductRight called for non-sparse m2";
         esac;
     rows = isRowMajor? m2;
     data = array (map \(new double[size.columns]) [1..size.rows]);
     for [0..size.rows - 1] do i':
         r = getRow i' m1;
         var p = 0;
-        for [0..length indices - 1] do ix:
-            ix == pointers[p+1] loop (p := p + 1);
-            i = if rows then p else indices[ix] fi;
-            j = if rows then indices[ix] else p fi;
+        for [0..ivec.length indices - 1] do ix:
+            ix == ivec.at pointers (p+1) loop (p := p + 1);
+            i = if rows then p else ivec.at indices ix fi;
+            j = if rows then ivec.at indices ix else p fi;
             data[i'][j] := data[i'][j] + (vec.at values ix) * (vec.at r i);
         done;
     done;
@@ -618,21 +619,21 @@
             esac;
         rows = isRowMajor? m1;
         var p = 0;
-        pindices = new int[length indices];
-        for [0..length indices - 1] do ix:
-            ix == pointers[p+1] loop (p := p + 1);
+        pindices = new int[ivec.length indices];
+        for [0..ivec.length indices - 1] do ix:
+            ix == ivec.at pointers (p+1) loop (p := p + 1);
             pindices[ix] := p;
         done;
         entries =
            (map do j':
                 cs = sparseSlice j' d;
                 hin = mapIntoHash
-                   (at cs.indices) (vec.at cs.values)
-                   [0..length cs.indices - 1];
+                   (ivec.at cs.indices) (vec.at cs.values)
+                   [0..ivec.length cs.indices - 1];
                 hout = [:];
-                for [0..length indices - 1] do ix:
-                    i = if rows then pindices[ix] else indices[ix] fi;
-                    j = if rows then indices[ix] else pindices[ix] fi;
+                for [0..ivec.length indices - 1] do ix:
+                    i = if rows then pindices[ix] else ivec.at indices ix fi;
+                    j = if rows then ivec.at indices ix else pindices[ix] fi;
                     if j in hin then
                         p = (vec.at values ix) * hin[j];
                         hout[i] := p + (if i in hout then hout[i] else 0 fi);
diff -r acc244cab1b7 -r bdee7ac8cbcd src/may/vector.yeti
--- a/src/may/vector.yeti	Fri Mar 21 10:33:04 2014 +0000
+++ b/src/may/vector.yeti	Fri Mar 21 16:22:31 2014 +0000
@@ -31,7 +31,10 @@
     a);
 
 /// Return a vector of length n, containing all ones.
-ones n = consts 1.0 n;
+ones n is number -> ~double[] =
+   (a = zeros n;
+    Arrays#fill(a, 1.0);
+    a);
 
 /// Return a vector of the values in the given list.
 fromList l is list?<number> -> ~double[] =
diff -r acc244cab1b7 -r bdee7ac8cbcd src/may/vector/ivector.yeti
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/may/vector/ivector.yeti	Fri Mar 21 16:22:31 2014 +0000
@@ -0,0 +1,220 @@
+
+/**
+ * Vector of int. A typesafe, immutable wrapper around a Java primitive
+ * array of int. Much more limited than vector of double.
+ */
+
+module may.vector.ivector;
+
+import java.util: Arrays;
+
+/// Return a vector of n zeros.
+zeros n is number -> ~int[] =
+    new int[n];
+
+/// Return a vector of length n, containing all m.
+consts m n is number -> number -> ~int[] =
+   (a = zeros n;
+    Arrays#fill(a, m);
+    a);
+
+/// Return a vector of length n, containing all ones.
+ones n is number -> ~int[] =
+   (a = zeros n;
+    Arrays#fill(a, 1);
+    a);
+
+/// Return a vector of the values in the given list.
+fromList l is list?<number> -> ~int[] =
+    l as ~int[];
+
+/// Return the given vector as a list.
+list' a is ~int[] -> list<number> =
+    list a;
+
+/// Return the given vector as a Yeti array.
+array' a is ~int[] -> array<number> =
+    array a;
+
+/// Return the length of the given vector.
+length' v is ~int[] -> number =
+    length v;
+
+/// Return true if the given vector is empty (has length 0).
+empty?' =
+    empty? . list';
+
+/// Return element n in the given vector v. (The function name and
+/// argument order are chosen for symmetry with the similar standard
+/// library array function.)
+at' v n is ~int[] -> number -> number =
+    v[n];
+
+/// Return true if the given vectors are equal, using the standard ==
+/// comparator on their elements.
+equal v1 v2 is ~int[] -> ~int[] -> boolean =
+    list v1 == list v2;
+
+/// Return true if the given vectors are equal, when applying the
+/// given numerical comparator to each element.
+equalUnder comparator v1 v2 =
+    length' v1 == length' v2 and
+        all id (map2 comparator (list v1) (list v2));
+
+/// Return another copy of the given vector.
+copyOf v is ~int[] -> ~int[] =
+    Arrays#copyOf(v, length v);
+
+/// Return the given vector as a primitive array. Modifying the
+/// array contents will not affect the original vector.
+primitive = copyOf;
+
+/// Return a new vector containing a subset of the elements of the
+/// given vector, from index start (inclusive) to index finish
+/// (exclusive). (The function name and argument order are chosen for
+/// symmetry with the standard library slice and strSlice functions.)
+slice v start finish is ~int[] -> number -> number -> ~int[] =
+   (len = length v;
+    if start == 0 and finish == len then v
+    elif start < 0 then slice v 0 finish
+    elif start > len then slice v len finish
+    elif finish < start then slice v start start
+    elif finish > len then slice v start len
+    else
+        Arrays#copyOfRange(v, start, finish);
+    fi);
+
+/// Return a vector of length n, containing the contents of the given
+/// vector v. If v is longer than n, the contents will be truncated;
+/// if shorter, they will be padded with zeros.
+resizedTo n v is number -> ~int[] -> ~int[] =
+    if n == length v then v;
+    else Arrays#copyOf(v, n);
+    fi;
+
+/// Return a vector that is the reverse of the given vector.  Name
+/// chosen (in preference to passive "reversed") for symmetry with the
+/// standard library list reverse function.
+reverse v is ~int[] -> ~int[] =
+   (len = length v;
+    a = new int[len];
+    for [0..len-1] do i:
+        a[len-i-1] := v[i];
+    done;
+    a);
+
+/// Return a single new vector that contains the contents of all the
+/// given vectors, in order. (Unlike the standard module list concat
+/// function, this one cannot be lazy.)
+concat' vv is list?<~int[]> -> ~int[] =
+    case vv of
+        [v]: v;
+        v0::rest:
+           (len = sum (map length' vv);
+            vout = Arrays#copyOf(v0 is ~int[], len);
+            var base = length' v0;
+            for rest do v: 
+                vlen = length' v;
+                System#arraycopy(v, 0, vout, base, vlen);
+                base := base + vlen;
+            done;
+            vout);
+        _: zeros 0;
+    esac;
+
+/// Return a single new vector that contains the contents of the given
+/// vector, repeated n times. The vector will therefore have length n
+/// times the length of v.
+repeated v n is ~int[] -> number -> ~int[] =
+    concat' (map \(v) [1..n]);
+
+max' v is ~int[] -> number = 
+   (var mx = 0;
+    for [0..length v - 1] do i:
+        if i == 0 or v[i] > mx then
+            mx := v[i];
+        fi
+    done;
+    mx);
+
+maxindex v is ~int[] -> number =
+   (var mx = 0;
+    var mi = -1;
+    for [0..length v - 1] do i:
+        if i == 0 or v[i] > mx then
+            mx := v[i];
+            mi := i;
+        fi
+    done;
+    mi);
+
+min' v is ~int[] -> number = 
+   (var mn = 0;
+    for [0..length v - 1] do i:
+        if i == 0 or v[i] < mn then
+            mn := v[i];
+        fi
+    done;
+    mn);
+
+minindex v is ~int[] -> number =
+   (var mn = 0;
+    var mi = -1;
+    for [0..length v - 1] do i:
+        if i == 0 or v[i] < mn then
+            mn := v[i];
+            mi := i;
+        fi
+    done;
+    mi);
+
+typedef opaque ivector_t = ~int[];
+
+{
+    zeros,
+    consts,
+    ones,
+    vector v = v,
+    primitive,
+    fromList,
+    list = list',
+    array = array',
+    length = length',
+    empty? = empty?',
+    at = at',
+    equal,
+    equalUnder,
+    slice,
+    resizedTo,
+    reverse,
+    repeated,
+    concat = concat',
+    max = max',
+    min = min',
+    maxindex,
+    minindex,
+} as {
+    zeros is number -> ivector_t,
+    consts is number -> number -> ivector_t,
+    ones is number -> ivector_t,
+    vector is ~int[] -> ivector_t,
+    primitive is ivector_t -> ~int[],
+    fromList is list?<number> -> ivector_t,
+    list is ivector_t -> list<number>,
+    array is ivector_t -> array<number>,
+    length is ivector_t -> number,
+    empty? is ivector_t -> boolean,
+    at is ivector_t -> number -> number,
+    equal is ivector_t -> ivector_t -> boolean,
+    equalUnder is (number -> number -> boolean) -> ivector_t -> ivector_t -> boolean,
+    slice is ivector_t -> number -> number -> ivector_t,
+    resizedTo is number -> ivector_t -> ivector_t,
+    reverse is ivector_t -> ivector_t,
+    repeated is ivector_t -> number -> ivector_t,
+    concat is list?<ivector_t> -> ivector_t,
+    max is ivector_t -> number,
+    min is ivector_t -> number,
+    maxindex is ivector_t -> number,
+    minindex is ivector_t -> number,
+}
+