May

/**

 * 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.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 =

    jsOfList (map do r: jsOfList (map jsOfNum (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,

}