Chris@5: Chris@94: module yetilab.matrix.matrix; Chris@94: Chris@94: // A matrix is an array of fvectors (i.e. primitive double[]s). Chris@94: Chris@94: // A matrix can be either RowMajor, akin to a C multidimensional array Chris@96: // in which each row is a separate fvector, or ColumnMajor, akin to a Chris@94: // FORTAN multidimensional array in which each column is a separate Chris@96: // fvector. The default is ColumnMajor. Storage order is an efficiency Chris@94: // concern only, all operations behave identically regardless. (The Chris@94: // transpose function just switches the row/column order without Chris@94: // moving the elements.) Chris@18: Chris@93: vec = load yetilab.block.fvector; Chris@96: block = load yetilab.block.block; Chris@99: bf = load yetilab.block.blockfuncs; Chris@9: Chris@96: make d = { Chris@96: get data () = d, Chris@96: get size () = Chris@96: case d of Chris@96: RowM r: Chris@96: major = length r; Chris@96: { Chris@96: rows = major, Chris@96: columns = if major > 0 then vec.length r[0] else 0 fi, Chris@96: }; Chris@96: ColM c: Chris@96: major = length c; Chris@96: { Chris@96: rows = if major > 0 then vec.length c[0] else 0 fi, Chris@96: columns = major, Chris@96: }; Chris@94: esac, Chris@96: getColumn j = Chris@96: case d of Chris@96: RowM rows: block.fromList (map do i: getAt i j done [0..length rows-1]); Chris@96: ColM cols: block.block cols[j]; Chris@94: esac, Chris@96: getRow i = Chris@96: case d of Chris@96: RowM rows: block.block rows[i]; Chris@96: ColM cols: block.fromList (map do j: getAt i j done [0..length cols-1]); Chris@94: esac, Chris@96: getAt row col = Chris@96: case d of Chris@96: RowM rows: r = rows[row]; (r is ~double[])[col]; Chris@96: ColM cols: c = cols[col]; (c is ~double[])[row]; Chris@94: esac, Chris@96: setAt row col n = Chris@96: case d of Chris@96: RowM rows: r = rows[row]; (r is ~double[])[col] := n; Chris@96: ColM cols: c = cols[col]; (c is ~double[])[row] := n; Chris@95: esac, Chris@95: get isRowMajor? () = Chris@96: case d of Chris@96: RowM _: true; Chris@96: ColM _: false; Chris@95: esac, Chris@94: }; Chris@94: Chris@98: newStorage { rows, columns } = Chris@97: if rows < 1 then array [] Chris@98: else array (map \(vec.zeros rows) [1..columns]) Chris@97: fi; Chris@94: Chris@98: zeroMatrix { rows, columns } = Chris@98: make (ColM (newStorage { rows, columns })); Chris@5: Chris@98: generate f { rows, columns } = Chris@98: (m = newStorage { rows, columns }; Chris@98: for [0..columns-1] do col: Chris@94: for [0..rows-1] do row: Chris@94: m[col][row] := f row col; Chris@5: done; Chris@5: done; Chris@96: make (ColM m)); Chris@5: Chris@20: constMatrix n = generate do row col: n done; Chris@20: randomMatrix = generate do row col: Math#random() done; Chris@5: identityMatrix = constMatrix 1; Chris@5: Chris@96: width m = m.size.columns; Chris@96: height m = m.size.rows; Chris@6: Chris@100: transposed m = Chris@100: make Chris@100: (case m.data of Chris@100: RowM d: ColM d; Chris@100: ColM d: RowM d; Chris@100: esac); Chris@100: Chris@100: flipped m = Chris@100: if m.isRowMajor? then Chris@100: generate do row col: m.getAt row col done m.size; Chris@100: else Chris@100: transposed Chris@100: (generate do row col: m.getAt col row done Chris@100: { rows = m.size.columns, columns = m.size.rows }); Chris@100: fi; Chris@100: Chris@100: // Matrices with different storage order but the same contents are Chris@100: // equal (but comparing them is slow) Chris@97: equal m1 m2 = Chris@100: if m1.isRowMajor? != m2.isRowMajor? Chris@100: then equal (flipped m1) m2; Chris@100: else Chris@100: compare d1 d2 = all id (map2 vec.equal d1 d2); Chris@100: case m1.data of Chris@100: RowM d1: case m2.data of RowM d2: compare d1 d2; _: false; esac; Chris@100: ColM d1: case m2.data of ColM d2: compare d1 d2; _: false; esac; Chris@100: esac Chris@100: fi; Chris@97: Chris@95: copyOf m = Chris@95: (copyOfData d = (array (map vec.copyOf d)); Chris@96: make Chris@95: (case m.data of Chris@96: RowM d: RowM (copyOfData d); Chris@96: ColM d: ColM (copyOfData d); Chris@95: esac)); Chris@6: Chris@96: newMatrix type data is RowMajor () | ColumnMajor () -> list?> -> 'a = Chris@96: (tagger = case type of RowMajor (): RowM; ColumnMajor (): ColM esac; Chris@97: if empty? data or empty? (head data) Chris@98: then zeroMatrix { rows = 0, columns = 0 } Chris@96: else make (tagger (array (map vec.vector data))) Chris@96: fi); Chris@96: Chris@96: newRowVector data = Chris@96: newMatrix (RowMajor ()) [data]; Chris@96: Chris@96: newColumnVector data = Chris@96: newMatrix (ColumnMajor ()) [data]; Chris@8: Chris@98: scaled factor m = Chris@98: generate do row col: factor * m.getAt row col done m.size; Chris@98: Chris@98: sum' m1 m2 = Chris@98: if m1.size != m2.size Chris@98: then failWith "Matrices are not the same size: \(m1.size), \(m2.size)"; Chris@98: else Chris@98: generate do row col: m1.getAt row col + m2.getAt row col done m1.size; Chris@98: fi; Chris@98: Chris@98: product m1 m2 = Chris@98: if m1.size.columns != m2.size.rows Chris@98: then failWith "Matrix dimensions incompatible: \(m1.size), \(m2.size) (\(m1.size.columns != m2.size.rows)"; Chris@98: else Chris@98: //!!! super-slow! Chris@98: generate do row col: Chris@99: bf.sum (bf.multiply (m1.getRow row) (m2.getColumn col)) Chris@98: done { rows = m1.size.rows, columns = m2.size.columns } Chris@98: fi; Chris@98: Chris@5: { Chris@97: constMatrix, randomMatrix, zeroMatrix, identityMatrix, Chris@96: generate, Chris@96: width, height, Chris@97: equal, Chris@20: copyOf, Chris@15: transposed, Chris@96: flipped, Chris@98: scaled, Chris@98: sum = sum', product, Chris@98: newMatrix, newRowVector, newColumnVector, Chris@5: } Chris@5: