view yetilab/matrix/matrix.yeti @ 99:9832210dc42c

More tests (some failing) and a bit more implementation

module yetilab.matrix.matrix;

// A matrix is an array of fvectors (i.e. primitive double[]s).

// A matrix can be either RowMajor, akin to a C multidimensional array
// in which each row is a separate fvector, or ColumnMajor, akin to a
// FORTAN multidimensional array in which each column is a separate
// fvector. The default is ColumnMajor. Storage order is an efficiency
// concern only, all operations behave identically regardless.  (The
// transpose function just switches the row/column order without
// moving the elements.)

vec = load yetilab.block.fvector;
block = load yetilab.block.block;
bf = load yetilab.block.blockfuncs;

make d = {
    get data () = d,
    get size () =
        case d of
        RowM r:
            major = length r;
            { 
                rows = major, 
                columns = if major > 0 then vec.length r[0] else 0 fi,
            };
        ColM c:
            major = length c;
            { 
                rows = if major > 0 then vec.length c[0] else 0 fi,
                columns = major, 
            };
        esac,
    getColumn j =
        case d of
        RowM rows: block.fromList (map do i: getAt i j done [0..length rows-1]);
        ColM cols: block.block cols[j];
        esac,
    getRow i =
        case d of
        RowM rows: block.block rows[i];
        ColM cols: block.fromList (map do j: getAt i j done [0..length cols-1]);
        esac,
    getAt row col =
        case d of
        RowM rows: r = rows[row]; (r is ~double[])[col];
        ColM cols: c = cols[col]; (c is ~double[])[row];
        esac,
    setAt row col n =
        case d of
        RowM rows: r = rows[row]; (r is ~double[])[col] := n;
        ColM cols: c = cols[col]; (c is ~double[])[row] := n;
        esac,
    get isRowMajor? () =
        case d of
        RowM _: true;
        ColM _: false;
        esac,
    };

newStorage { rows, columns } = 
    if rows < 1 then array []
    else array (map \(vec.zeros rows) [1..columns])
    fi;

zeroMatrix { rows, columns } = 
    make (ColM (newStorage { rows, columns }));

generate f { rows, columns } =
   (m = newStorage { rows, columns };
    for [0..columns-1] do col:
        for [0..rows-1] do row:
            m[col][row] := f row col;
        done;
    done;
    make (ColM m));

constMatrix n = generate do row col: n done;
randomMatrix = generate do row col: Math#random() done;
identityMatrix = constMatrix 1;

width m = m.size.columns;
height m = m.size.rows;

//!!! should matrices with the same content but different storage order be equal?
equal m1 m2 =
   (compare d1 d2 =
        all id (map2 vec.equal d1 d2);
    case m1.data of
    RowM d1:
        case m2.data of RowM d2: compare d1 d2; _: false; esac;
    ColM d1:
        case m2.data of ColM d2: compare d1 d2; _: false; esac;
    esac);

copyOf m =
   (copyOfData d = (array (map vec.copyOf d));
    make
       (case m.data of
        RowM d: RowM (copyOfData d);
        ColM d: ColM (copyOfData d);
        esac));

transposed m =
    make
       (case m.data of
        RowM d: ColM d;
        ColM d: RowM d;
        esac);

//!!! Change storage from column to row major order (or back
//again). Is there a word for this?
flipped m =
    if m.isRowMajor? then id else transposed fi
       (generate do row col: m.getAt col row done
           { rows = m.size.columns, columns = m.size.rows });

newMatrix type data is RowMajor () | ColumnMajor () -> list?<list?<number>> -> 'a =
   (tagger = case type of RowMajor (): RowM; ColumnMajor (): ColM esac;
    if empty? data or empty? (head data)
    then zeroMatrix { rows = 0, columns = 0 }
    else make (tagger (array (map vec.vector data)))
    fi);

newRowVector data = 
    newMatrix (RowMajor ()) [data];

newColumnVector data = 
    newMatrix (ColumnMajor ()) [data];

scaled factor m =
    generate do row col: factor * m.getAt row col done m.size;

sum' m1 m2 =
    if m1.size != m2.size
    then failWith "Matrices are not the same size: \(m1.size), \(m2.size)";
    else
        generate do row col: m1.getAt row col + m2.getAt row col done m1.size;
    fi;

product m1 m2 =
    if m1.size.columns != m2.size.rows
    then failWith "Matrix dimensions incompatible: \(m1.size), \(m2.size) (\(m1.size.columns != m2.size.rows)";
    else
    //!!! super-slow!
        generate do row col:
            bf.sum (bf.multiply (m1.getRow row) (m2.getColumn col))
        done { rows = m1.size.rows, columns = m2.size.columns }
    fi;

{
constMatrix, randomMatrix, zeroMatrix, identityMatrix,
generate,
width, height,
equal,
copyOf,
transposed,
flipped,
scaled,
sum = sum', product,
newMatrix, newRowVector, newColumnVector,
}