view yetilab/matrix/matrix.yeti @ 97:d5fc902dcc3f

Initial matrix tests

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;

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 cols = 
    if rows < 1 then array []
    else array (map \(vec.zeros rows) [1..cols])
    fi;

zeroMatrix rows cols = 
    make (ColM (newStorage rows cols));

generate f rows cols =
   (m = newStorage rows cols;
    for [0..cols-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;

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 m.size.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 0 0
    else make (tagger (array (map vec.vector data)))
    fi);

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

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

{
constMatrix, randomMatrix, zeroMatrix, identityMatrix,
generate,
width, height,
equal,
copyOf,
transposed,
flipped,
newMatrix, newRowVector, newColumnVector
}