view src/may/matrix.yeti @ 591:eb27901664cd

Check for consistent lengths in fromRows/fromColumns; test for same in mapRows/mapColumns

/**
 * Matrices. A matrix is a two-dimensional (NxM) container of
 * double-precision floating point values.
 *
 * A matrix may be dense or sparse.
 * 
 * A dense matrix (the default) is just a series of vectors, making up
 * the matrix "grid". The values may be stored in either column-major
 * order, in which case the series consists of one vector for each
 * column in the matrix, or row-major order, in which case the series
 * consists of one vector for each row. The default is column-major.
 * 
 * A sparse matrix has a more complex representation in which only the
 * non-zero values are stored. This is typically used for matrices
 * containing sparse data, that is, data in which most of the values
 * are zero: using a sparse representation is more efficient than a
 * dense one (in both time and memory) if the matrix is very large but
 * contains a relatively low proportion of non-zero values. Like dense
 * matrices, sparse ones may be column-major or row-major.
 * 
 * The choice of dense or sparse, row- or column-major is a question
 * of efficiency alone. All functions in this module should return the
 * same results regardless of how the matrices they operate on are
 * represented. However, differences in performance can be very large
 * and it is often worth converting matrices to a different storage
 * format if you know they can be more efficiently manipulated that
 * way. For example, multiplying two matrices is fastest if the first
 * is in column-major and the second in row-major order.
 * 
 * Use the isRowMajor? and isSparse? functions to query the storage
 * format of a matrix; use the flipped function to convert between
 * column-major and row-major storage; and use toSparse and toDense to
 * convert between sparse and dense storage.
 *
 * Note that the matrix size is preserved even if at least one
 * dimension is zero. That is, it is legal to have matrices of size
 * 0x0, 0x4, 1x0 etc, and they are distinct from each other.
 */

module may.matrix;

{ ceil, floor, random } = load may.mathmisc;

vec = load may.vector;

load yeti.experimental.json;

typedef opaque matrix_t = {
    size is { rows is number, columns is number },
    data is
        DenseRows array<vec.vector_t> | // array of rows
        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
            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
            extent is number // max pointers index + 1, i.e. number of rows
        }
};

size m = m.size;
width m = m.size.columns;
height m = m.size.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.data 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,
    });

nonEmptySlices d =
   (ne = array [];
    for [0..length d.pointers - 2] do i:
        if d.pointers[i] != d.pointers[i+1] then
            push ne i
        fi
    done;
    ne);

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 slice.values i;
        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 slice.values i;
    done;
    vec.vector dslice);

at' m row col =
    case m.data of
    DenseRows rows: r = rows[row]; vec.at r col;
    DenseCols cols: c = cols[col]; vec.at c row;
    SparseCSR data: fromSlice row col data;
    SparseCSC data: fromSlice col row data;
    esac;

//!!! better as getXx or just xx?

//!!! arguably getRow, getColumn, getDiagonal should have m as first arg for symmetry with at
getColumn j m =
    case m.data of
    DenseCols cols: cols[j];
    SparseCSC data: filledSlice j data;
    _: vec.fromList (map do i: at' m i j done [0..height m - 1]);
    esac;

getRow i m =
    case m.data of
    DenseRows rows: rows[i];
    SparseCSR data: filledSlice i data; 
    _: vec.fromList (map do j: at' m i j done [0..width m - 1]);
    esac;

getDiagonal k m =
   (ioff = if k < 0 then -k else 0 fi;
    joff = if k > 0 then  k else 0 fi;
    n = min (width m - joff) (height m - ioff);
    vec.fromList (map do i: at' m (i + ioff) (i + joff) done [0..n - 1]));

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.data of
    DenseRows _: true;
    DenseCols _: false;
    SparseCSR _: true;
    SparseCSC _: false;
    esac;

isSparse? m =
    case m.data of
    DenseRows _: false;
    DenseCols _: false;
    SparseCSR _: true;
    SparseCSC _: true;
    esac;

taggerForTypeOf m =
    if isRowMajor? m then Rows
    else Columns
    fi;

taggerForFlippedTypeOf m =
    if isRowMajor? m then Columns
    else Rows
    fi;

flippedSize { rows, columns } = { rows = columns, columns = rows };

newColumnMajorStorage { rows, columns } = 
    array (map \(vec.zeros rows) [1..columns]);

zeroMatrix size = 
    {
        size,
        data = DenseCols (newColumnMajorStorage size)
    };

zeroMatrixWithTypeOf m size = 
    {
        size,
        data =
            if isRowMajor? m then
                DenseRows (newColumnMajorStorage (flippedSize size));
            else
                DenseCols (newColumnMajorStorage size);
            fi
    };

zeroSizeMatrix () = zeroMatrix { rows = 0, columns = 0 };

newMatrix size d =
    case d of
    Rows rr: 
        if (length rr) != size.rows then
            failWith "Wrong number of rows in row-major newMatrix (\(length rr), size calls for \(size.rows))";
        elif not (all do r: vec.length r == size.columns done rr) then
            failWith "Wrong or inconsistent number of columns in rows in row-major newMatrix (\(map vec.length rr)), size calls for \(size.columns))";
        else
            {
                size,
                data = DenseRows (array rr)
            }
        fi;
    Columns cc: 
        if (length cc) != size.columns then
            failWith "Wrong number of columns in column-major newMatrix (\(length cc), size calls for \(size.columns))";
        elif not (all do c: vec.length c == size.rows done cc) then
            failWith "Wrong or inconsistent number of rows in in columns in column-major newMatrix (\(map vec.length cc)), size calls for \(size.rows))";
        else
            {
                size,
                data = DenseCols (array cc)
            }
        fi;
    esac;

newMatrixMatching m d = newMatrix (size m) ((taggerForTypeOf m) d);

generate f { rows, columns } =
   (m = array (map \(new double[rows]) [1..columns]);
    for [0..columns-1] do col:
        for [0..rows-1] do row:
            m[col][row] := f row col;
        done;
    done;
    {
        size = { rows, columns },
        data = DenseCols (array (map vec.vector m))
    });

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
           (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.data of
    SparseCSC d: swapij (enumerate d);
    SparseCSR d: enumerate d;
     _: [];
    esac);

enumerateDense m =
   (enumerate d =
        concat
           (map do i:
                vv = d[i];
                map2 do j v: { i, j, v } done
                    [0..vec.length vv - 1]
                    (vec.list vv);
                done [0..length d - 1]);
    case m.data of
    DenseCols c: swapij (enumerate c);
    DenseRows r: enumerate r;
     _: [];
    esac);

enumerate m =
    if isSparse? m then enumerateSparse m else enumerateDense m fi;

// Make a sparse matrix from entries whose i, j values are known to be
// within range
newSparse size d =
   (isRow = case d of Rows _: true; Columns _: false esac;
    data = case d of Rows rr: rr; Columns cc: cc 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
               (filter do d: d.v != 0 done 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;
    pointers = array [0];
    setArrayCapacity pointers (size.rows + 1);
    fillPointers n i data =
        if n < majorSize then
            case data of
            d::rest:
               (for [n..d-1] \(push pointers i);
                fillPointers d (i+1) rest);
             _:
                for [n..majorSize-1] \(push pointers i);
            esac;
        fi;
    fillPointers 0 0 (map (.maj) ordered);
    {
        size,
        data = tagger {
            values = vec.fromList (map (.v) ordered),
            indices = array (map (.min) ordered),
            pointers,
            extent = minorSize,
        }
    });

newSparseMatching m d = newSparse (size m) ((taggerForTypeOf m) d);

// Make a sparse matrix from entries that may contain out-of-range
// cells which need to be filtered out. This is the public API for
// newSparse and is also used to discard out-of-range cells from
// resizedTo.
//!!! doc: i is row number, j is column number (throughout, for sparse stuff). Would calling them row/column be better?
//!!! doc: Rows/Columns determines the storage order, the input data are treated the same either way (perhaps this does mean row/column would be better than i/j)
newSparseMatrix size d =
   (tagger = case d of Rows _: Rows; Columns _: Columns esac;
    data = case d of Rows rr: rr; Columns cc: cc esac;
    data = filter
        do { i, j, v }:
            i == int i and i >= 0 and i < size.rows and 
            j == int j and j >= 0 and j < size.columns
        done data;
    newSparse size (tagger data));

toSparse m =
    if isSparse? m then m
    else newSparseMatching m (enumerateDense m);
    fi;

toDense m =
    {
        size = (size m),
        data = 
            if not (isSparse? m) then m.data
            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;
randomMatrix = generate do row col: random () done;
identityMatrix = constMatrix 1;

transposed m =
    {
        size = flippedSize (size m),
        data = 
            case m.data of
            DenseRows d: DenseCols d;
            DenseCols d: DenseRows d;
            SparseCSR d: SparseCSC d;
            SparseCSC d: SparseCSR d;
            esac
    };

flipped m =
    if isSparse? m then
        newSparse (size m) ((taggerForFlippedTypeOf m) (enumerateSparse m))
    else
        if isRowMajor? m then
            generate do row col: at' m row col done (size m);
        else
            transposed
               (generate do row col: at' m col row done (flippedSize (size m)));
        fi
    fi;

toRowMajor m =
    if isRowMajor? m then m else flipped m fi;

toColumnMajor m =
    if not isRowMajor? m then m else flipped m fi;

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.data of
    DenseRows d1:
        case m2.data of DenseRows d2: compareVecLists d1 d2; _: false; esac;
    DenseCols d1:
        case m2.data of DenseCols d2: compareVecLists d1 d2; _: false; esac;
    SparseCSR d1:
        case m2.data of SparseCSR d2: compareSparse d1 d2; _: false; esac;
    SparseCSC d1:
        case m2.data 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' 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
        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' comparator (vec.equalUnder comparator);

equal =
    equal' (==) vec.equal;

fromRows rows =
   (if any do r: vec.length r != vec.length (head rows) done rows then
        failWith "Inconsistent row lengths in fromRows (\(map vec.length rows))";
    fi;
    {
        size = { 
            rows = length rows, 
            columns = 
                if empty? rows then 0
                else vec.length (head rows) 
                fi,
        },
        data = DenseRows (array rows)
    });

fromColumns cols =
   (if any do c: vec.length c != vec.length (head cols) done cols then
        failWith "Inconsistent column lengths in fromColumns (\(map vec.length cols))";
    fi;
    {
        size = { 
            columns = length cols, 
            rows = 
                if empty? cols then 0
                else vec.length (head cols) 
                fi,
        },
        data = DenseCols (array cols)
    });

fromLists data = 
    case data of
    Rows rr: fromRows (map vec.fromList rr);
    Columns cc: fromColumns (map vec.fromList cc);
    esac;

newRowVector data = //!!! NB does not copy data
    fromRows (array [data]);

newColumnVector data = //!!! NB does not copy data
    fromColumns (array [data]);

denseLinearOp op m1 m2 =
    if isRowMajor? m1 then
        newMatrixMatching m1
           (map2 do c1 c2: op c1 c2 done (asRows m1) (asRows m2));
    else
        newMatrixMatching m1
           (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));
    newSparseMatching m1 entries);

sum' mm =
    case mm of
    m1::m2::rest:
        sum' 
           (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
                add2 v1 v2 = vec.add [v1,v2];
                denseLinearOp add2 m1 m2;
            fi :: rest);
    [m1]: m1;
    _: failWith "Empty argument list";
    esac;
    
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
        denseLinearOp vec.subtract m1 m2;
    fi;

scaled factor m =
    if isSparse? m then
        newSparseMatching m
           (map do { i, j, v }: { i, j, v = factor * v } done (enumerate m))
    elif isRowMajor? m then
        newMatrixMatching m (map (vec.scaled factor) (asRows m));
    else
        newMatrixMatching m (map (vec.scaled factor) (asColumns m));
    fi;

abs' m =
    if isSparse? m then
        newSparseMatching m
           (map do { i, j, v }: { i, j, v = abs v } done (enumerate m))
    elif isRowMajor? m then
        newMatrixMatching m (map vec.abs (asRows m));
    else
        newMatrixMatching m (map vec.abs (asColumns m));
    fi;

negative m =
    if isSparse? m then
        newSparseMatching m
           (map do { i, j, v }: { i, j, v = (-v) } done (enumerate m))
    elif isRowMajor? m then
        newMatrixMatching m (map vec.negative (asRows m));
    else
        newMatrixMatching m (map vec.negative (asColumns m));
    fi;

//!!! doc: filter by predicate, always returns sparse matrix
filter' f m =
    newSparseMatching m
       (map do { i, j, v }: { i, j, v = if f v then v else 0 fi } done
           (enumerate m));

any' f m =
    any f (map (.v) (enumerate m));

all' f m =
    all f (map (.v) (enumerate m));

sparseProductLeft size m1 m2 =
   ({ values, indices, pointers } = 
        case m1.data of
        SparseCSR d: d;
        SparseCSC d: d;
        _: failWith "sparseProductLeft called for non-sparse m1";
        esac;
    rows = isRowMajor? m1;
    data = array (map \(new double[size.rows]) [1..size.columns]);
    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;
            data[j'][i] := data[j'][i] + (vec.at values ix) * (vec.at c j);
        done;
    done;
    newMatrix size (Columns (array (map vec.vector (list data)))));

sparseProductRight size m1 m2 =
   ({ values, indices, pointers } = 
        case m2.data of
        SparseCSR d: d;
        SparseCSC d: d;
        _: failWith "sparseProductLeft called for non-sparse m1";
        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;
            data[i'][j] := data[i'][j] + (vec.at values ix) * (vec.at r i);
        done;
    done;
    newMatrix size (Rows (array (map vec.vector (list data)))));

sparseProduct size m1 m2 =
    case m2.data of
    SparseCSC d:
       ({ values, indices, pointers } =
            case m1.data of
            SparseCSR d1: d1;
            SparseCSC d1: d1;
            _: failWith "sparseProduct called for non-sparse matrices";
            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[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];
                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;
                    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);
                    fi;
                done;
                map do i:
                    { i, j = j', v = hout[i] }
                done (keys hout);
            done (nonEmptySlices d));
        newSparse size (Columns (concat entries)));
    SparseCSR _:
        sparseProduct size m1 (flipped m2);
     _: failWith "sparseProduct called for non-sparse matrices";
    esac;

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] := vec.sum (vec.multiply [row, getColumn j m2]);
        done;
    done;
    newMatrix size (Columns (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 
        size = { rows = (size m1).rows, columns = (size m2).columns };
        if isSparse? m1 then
            if isSparse? m2 then
                sparseProduct size m1 m2
            else
                sparseProductLeft size m1 m2
            fi
        elif isSparse? m2 then
            sparseProductRight size m1 m2
        else
            denseProduct size m1 m2
        fi;
    fi;

entryWiseProduct mm =
    case mm of
    m1::m2::rest:
        entryWiseProduct
           (if (size m1) != (size m2)
            then failWith "Matrices are not the same size: \(size m1), \(size m2)";
            else 
                if isSparse? m1 then
                    newSparse (size m1)
                       ((taggerForTypeOf m1)
                           (map do { i, j, v }: { i, j, v = v * (at' m2 i j) } done
                               (enumerateSparse m1)))
                elif isSparse? m2 then
                    entryWiseProduct (m2::m1::rest)
                else
                    if isRowMajor? m1 then
                        fromRows (array (map2 do v1 v2: vec.multiply [v1,v2] done
                           (asRows m1) (asRows m2)));
                    else
                        fromColumns (array (map2 do v1 v2: vec.multiply [v1,v2] done
                           (asColumns m1) (asColumns m2)));
                    fi
                fi
            fi :: rest);
    [m1]: m1;
    _: failWith "Empty argument list";
    esac;

entryWiseDivide m1 m2 =
    if (size m1) != (size m2)
    then failWith "Matrices are not the same size: \(size m1), \(size m2)";
    else 
        if isSparse? m1 then
            newSparse (size m1)
               ((taggerForTypeOf m1)
                   (map do { i, j, v }: { i, j, v = v / (at' m2 i j) } done
                       (enumerateSparse m1)))
        // For m2 to be sparse makes no sense (divide by zero all over
        // the shop).
        else
            if isRowMajor? m1 then
                fromRows (array (map2 vec.divide (asRows m1) (asRows m2)));
            else
                fromColumns (array (map2 vec.divide (asColumns m1) (asColumns m2)));
            fi
        fi
    fi;

concatAgainstGrain tagger getter counter mm =
   (n = counter (size (head mm));
    tagger (array
       (map do i:
           vec.concat (map (getter i) mm)
           done [0..n-1])));

concatWithGrain tagger getter counter mm =
    tagger (array
       (concatMap do m:
           n = counter (size m);
           map do i: getter i m done [0..n-1]
        done mm));

sparseConcat direction first mm =
   (dimension d f = if direction == d then sum (map f mm) else f first fi;
    rows = dimension (Vertical ()) height;
    columns = dimension (Horizontal ()) width;
    entries ioff joff ui uj mm acc =
        case mm of 
        m::rest:
            entries
               (ioff + ui * height m)
               (joff + uj * width m)
                ui uj rest
               ((map do { i, j, v }: { i = i + ioff, j = j + joff, v }
                 done (enumerate m)) ++ acc);
         _: acc;
        esac;
    newSparse { rows, columns }
       ((taggerForTypeOf first)
           (if direction == Vertical () then entries 0 0 1 0 mm []
            else entries 0 0 0 1 mm [] fi)));

sumDimensions sumCounter checkCounter mm =
   (check = checkCounter (size (head mm));
    sum
       (map do m:
            s = size m;
            if (checkCounter s) != check then
                failWith "Matrix dimensions incompatible for concat (found \(map do m: checkCounter (size m) done mm) not all of which are \(check))";
            else
                sumCounter s;
            fi
        done mm));

concatHorizontal mm = //!!! doc: storage order is taken from first matrix in sequence; concat is obviously not lazy (unlike std module)
    case mm of
    [m]: m;
    first::rest:
       (w = sumDimensions (.columns) (.rows) mm;
        if all isSparse? mm then
            sparseConcat (Horizontal ()) first mm
        else
            row = isRowMajor? first;
            {
                size = { rows = height first, columns = w },
                data =
                    // horizontal, row-major: against grain with rows
                    // horizontal, col-major: with grain with cols
                    if row then concatAgainstGrain DenseRows getRow (.rows) mm;
                    else concatWithGrain DenseCols getColumn (.columns) mm;
                    fi
            };
        fi);
     _: zeroSizeMatrix ();
    esac;

concatVertical mm = //!!! doc: storage order is taken from first matrix in sequence; concat is obviously not lazy (unlike std module)
    case mm of
    [m]: m;
    first::rest:
       (h = sumDimensions (.rows) (.columns) mm;
        if all isSparse? mm then
            sparseConcat (Vertical ()) first mm
        else
            row = isRowMajor? first;
            {
                size = { rows = h, columns = width first },
                data = 
                    // vertical, row-major: with grain with rows
                    // vertical, col-major: against grain with cols
                    if row then concatWithGrain DenseRows getRow (.rows) mm;
                    else concatAgainstGrain DenseCols getColumn (.columns) mm;
                    fi,
            };
        fi);
     _: zeroSizeMatrix ();
    esac;

//!!! next two v. clumsy

//!!! doc note: argument order chosen for consistency with std module slice
//!!! NB always returns dense matrix, should have sparse version
rowSlice m start end = //!!! doc: storage order same as input
    if start < 0 then rowSlice m 0 end
    elif start > height m then rowSlice m (height m) end
    else
        if end < start then rowSlice m start start
        elif end > height m then rowSlice m start (height m)
        else
            if isRowMajor? m then
                newMatrix { rows = end - start, columns = width m }
                   (Rows
                       (array (map ((flip getRow) m) [start .. end - 1])))
            else 
                newMatrix { rows = end - start, columns = width m }
                   (Columns
                       (array (map do v: vec.slice v start end done (asColumns m))))
            fi;
        fi;
    fi;

//!!! doc note: argument order chosen for consistency with std module slice
//!!! NB always returns dense matrix, should have sparse version
columnSlice m start end = //!!! doc: storage order same as input
    if start < 0 then columnSlice m 0 end
    elif start > width m then columnSlice m (width m) end
    else
        if end < start then columnSlice m start start
        elif end > width m then columnSlice m start (width m)
        else
            if not isRowMajor? m then
                newMatrix { rows = height m, columns = end - start }
                   (Columns
                       (array (map ((flip getColumn) m) [start .. end - 1])))
            else 
                newMatrix { rows = height m, columns = end - start }
                   (Rows
                       (array (map do v: vec.slice v start end done (asRows m))))
            fi;
        fi;
    fi;

resizedTo newsize m =
   (if newsize == (size m) then
        m
    elif isSparse? m then
        // don't call newSparse directly: want to discard
        // out-of-range cells
        newSparseMatrix newsize ((taggerForTypeOf m) (enumerateSparse m))
    elif (height m) == 0 or (width m) == 0 then
        zeroMatrixWithTypeOf m newsize;
    else
        growrows = newsize.rows - (height m);
        growcols = newsize.columns - (width m);
        rowm = isRowMajor? m;
        resizedTo newsize
            if rowm and growrows < 0 then
                rowSlice m 0 newsize.rows
            elif (not rowm) and growcols < 0 then 
                columnSlice m 0 newsize.columns
            elif growrows < 0 then 
                rowSlice m 0 newsize.rows
            elif growcols < 0 then 
                columnSlice m 0 newsize.columns
            else
                if growrows > 0 then
                    concatVertical
                       [m, zeroMatrixWithTypeOf m ((size m) with { rows = growrows })]
                else
                    concatHorizontal
                       [m, zeroMatrixWithTypeOf m ((size m) with { columns = growcols })]
                fi
            fi
    fi);

//!!! doc: always dense
repeatedHorizontal n m =
   (if n == 1 then m
    else
        cols = asColumns m;
        fromColumns (fold do acc _: acc ++ cols done [] [1..n])
    fi);

//!!! doc: always dense
repeatedVertical n m =
   (if n == 1 then m
    else
        rows = asRows m;
        fromRows (fold do acc _: acc ++ rows done [] [1..n])
    fi);

//!!! doc: always dense
tiledTo newsize m =
    if newsize == size m then
        m
    elif (height m) == 0 or (width m) == 0 then
        zeroMatrixWithTypeOf m newsize;
    else    
        h = ceil (newsize.columns / (width m));
        v = ceil (newsize.rows / (height m));
        if isRowMajor? m then
            resizedTo newsize (repeatedHorizontal h (repeatedVertical v m))
        else
            resizedTo newsize (repeatedVertical v (repeatedHorizontal h m))
        fi
    fi;

minValue m =
    if width m == 0 or height m == 0 then 0
    elif isSparse? m then
        minv ll = fold min (head ll) (tail ll);
        minnz = minv (map (.v) (enumerate m));
        if minnz > 0 and nonZeroValues m < (width m * height m) then 0
        else minnz fi;
    elif isRowMajor? m then
        vec.min (vec.fromList (map vec.min (asRows m)));
    else
        vec.min (vec.fromList (map vec.min (asColumns m)));
    fi;

maxValue m =
    if width m == 0 or height m == 0 then 0
    elif isSparse? m then
        maxv ll = fold max (head ll) (tail ll);
        maxnz = maxv (map (.v) (enumerate m));
        if maxnz < 0 and nonZeroValues m < (width m * height m) then 0
        else maxnz fi;
    elif isRowMajor? m then
        vec.max (vec.fromList (map vec.max (asRows m)));
    else
        vec.max (vec.fromList (map vec.max (asColumns m)));
    fi;

total m = 
    if isSparse? m then
        fold (+) 0 (map (.v) (enumerateSparse m));
    elif isRowMajor? m then
        fold (+) 0 (map vec.sum (asRows m));
    else
        fold (+) 0 (map vec.sum (asColumns m));
    fi;

mapRows rf m =
    fromRows (map rf (asRows m));

mapColumns cf m =
    fromColumns (map cf (asColumns m));

format m =
    strJoin "\n"
       (chunk = 8;
        map do b:
            c0 = b * chunk;
            c1 = b * chunk + chunk - 1;
            c1 = if c1 > width m then width m else c1 fi;
            [ "\nColumns \(c0) to \(c1)\n",
              (map do row:
                   map do v:
                       strPad ' ' 10
                          (if v == 0 then "0.0"
                           elif abs v >= 1000.0 or abs v < 0.01 then
                               String#format("%.2E", [v as ~Double])
                           else
                               String#format("%5f", [v as ~Double])
                           fi);
                   done (vec.list row) |> strJoin "";
               done (asRows (columnSlice m c0 (c1 + 1))) |> strJoin "\n")
            ];
        done [0..floor(width m / chunk)] |> concat);

print' = println . format;
eprint' = eprintln . format;

json m =
    jsonList (map do r: jsonList (map jsonNum (vec.list r)) done (asRows m));

{
    size,
    width,
    height,
    density,
    nonZeroValues,
    at = at',
    getColumn,
    getRow,
    getDiagonal,
    isRowMajor?,
    isSparse?,
    generate,
    constMatrix,
    randomMatrix,
    zeroMatrix,
    identityMatrix,
    equal, //!!! if empty is empty?, why is equal not equal? ?
    equalUnder,
    transposed,
    flipped,
    toRowMajor,
    toColumnMajor,
    toSparse,
    toDense,
    scaled,
    minValue,
    maxValue,
    total,
    asRows,
    asColumns,
    sum = sum',
    difference,
    abs = abs',
    negative,
    filter = filter',
    all = all',
    any = any',
    product,
    entryWiseProduct,
    entryWiseDivide,
    resizedTo,
    tiledTo,
    repeatedHorizontal,
    repeatedVertical,
    concatHorizontal,
    concatVertical,
    rowSlice,
    columnSlice,
    mapRows,
    mapColumns,
    fromRows,
    fromColumns,
    fromLists,
    newMatrix,
    newRowVector,
    newColumnVector,
    newSparseMatrix,
    enumerate,
    format,
    print = print',
    eprint = eprint',
    json,
}
as
{
    size is matrix_t -> { rows is number, columns is number },
    width is matrix_t -> number,
    height is matrix_t -> number,
    density is matrix_t -> number,
    nonZeroValues is matrix_t -> number,
    at is matrix_t -> number -> number -> number,
    getColumn is number -> matrix_t -> vec.vector_t,
    getRow is number -> matrix_t -> vec.vector_t,
    getDiagonal is number -> matrix_t -> vec.vector_t,
    isRowMajor? is matrix_t -> boolean,
    isSparse? is matrix_t -> boolean,
    generate is (number -> number -> number) -> { rows is number, columns is number } -> matrix_t,
    constMatrix is number -> { rows is number, columns is number } -> matrix_t,
    randomMatrix is { rows is number, columns is number } -> matrix_t,
    zeroMatrix is { rows is number, columns is number } -> matrix_t, 
    identityMatrix is { rows is number, columns is number } -> matrix_t, 
    equal is matrix_t -> matrix_t -> boolean,
    equalUnder is (number -> number -> boolean) -> matrix_t -> matrix_t -> boolean,
    transposed is matrix_t -> matrix_t,
    flipped is matrix_t -> matrix_t, 
    toRowMajor is matrix_t -> matrix_t, 
    toColumnMajor is matrix_t -> matrix_t,
    toSparse is matrix_t -> matrix_t,
    toDense is matrix_t -> matrix_t,
    scaled is number -> matrix_t -> matrix_t,
    minValue is matrix_t -> number,
    maxValue is matrix_t -> number,
    total is matrix_t -> number,
    asRows is matrix_t -> list<vec.vector_t>, 
    asColumns is matrix_t -> list<vec.vector_t>,
    sum is list?<matrix_t> -> matrix_t,
    difference is matrix_t -> matrix_t -> matrix_t,
    abs is matrix_t -> matrix_t,
    negative is matrix_t -> matrix_t,
    filter is (number -> boolean) -> matrix_t -> matrix_t,
    all is (number -> boolean) -> matrix_t -> boolean,
    any is (number -> boolean) -> matrix_t -> boolean,
    product is matrix_t -> matrix_t -> matrix_t,
    entryWiseProduct is list?<matrix_t> -> matrix_t,
    entryWiseDivide is matrix_t -> matrix_t -> matrix_t,
    resizedTo is { rows is number, columns is number } -> matrix_t -> matrix_t,
    tiledTo is { rows is number, columns is number } -> matrix_t -> matrix_t,
    repeatedHorizontal is number -> matrix_t -> matrix_t,
    repeatedVertical is number -> matrix_t -> matrix_t,
    concatHorizontal is list<matrix_t> -> matrix_t,
    concatVertical is list<matrix_t> -> matrix_t,
    rowSlice is matrix_t -> number -> number -> matrix_t, 
    columnSlice is matrix_t -> number -> number -> matrix_t,
    mapRows is (vec.vector_t -> vec.vector_t) -> matrix_t -> matrix_t,
    mapColumns is (vec.vector_t -> vec.vector_t) -> matrix_t -> matrix_t,
    fromRows is list<vec.vector_t> -> matrix_t, 
    fromColumns is list<vec.vector_t> -> matrix_t, 
    fromLists is (Rows list<list<number>> | Columns list<list<number>>) -> matrix_t,
    newMatrix is { rows is number, columns is number } -> (Rows list<vec.vector_t> | Columns list<vec.vector_t>) -> matrix_t,
    newRowVector is vec.vector_t -> matrix_t, 
    newColumnVector is vec.vector_t -> matrix_t,
    newSparseMatrix is { rows is number, columns is number } -> (Rows list<{ i is number, j is number, v is number }> | Columns list<{ i is number, j is number, v is number }>) -> matrix_t,
    enumerate is matrix_t -> list<{ i is number, j is number, v is number }>,
    format is matrix_t -> string,
    print is matrix_t -> (),
    eprint is matrix_t -> (),
    json is matrix_t -> json,
}