From d6e5e148aab8257e24f63dc7ca06c51528fa54fe Mon Sep 17 00:00:00 2001
From: Mattia Giambirtone <nocturn9x@nocturn9x.space>
Date: Thu, 22 Dec 2022 14:55:26 +0100
Subject: [PATCH] Many fixes to matrix library and minor changes

---
 src/main.nim                |   2 +-
 src/nn/layer.nim            |  13 +-
 src/nn/util/activations.nim |   6 +-
 src/nn/util/matrix.nim      | 520 +++++++++++++++++++++---------------
 4 files changed, 309 insertions(+), 232 deletions(-)

diff --git a/src/main.nim b/src/main.nim
index fb15852..b1061e4 100644
--- a/src/main.nim
+++ b/src/main.nim
@@ -3,7 +3,7 @@ import nn/util/activations
 import nn/util/losses
 
 
-var net = newNeuralNetwork(@[2, 3, 2], activationFunc=newActivation(sigmoid, proc (x, y: float): float = 0.0), 
+var net = newNeuralNetwork(@[2, 3, 2], activationFunc=newActivation(sigmoid, func (x, y: float): float = 0.0), 
                           lossFunc=newLoss(mse, mse), weightRange=(-1.0, +1.0), learnRate=0.05)
 var prediction = net.predict(newMatrix[float](@[2.7, 3.0]))
 echo prediction
diff --git a/src/nn/layer.nim b/src/nn/layer.nim
index fa3244d..dd9d30a 100644
--- a/src/nn/layer.nim
+++ b/src/nn/layer.nim
@@ -74,13 +74,10 @@ proc compute*(self: Layer, data: Matrix[float]): Matrix[float] =
     ## Computes the output of a given layer with
     ## the given input data and returns it as a
     ## one-dimensional array
-    var sequence = newSeqOfCap[float](self.outputSize)
-    for i, weights in self.weights:
-        # This looks fancy, but it's just abstracting some of the
-        # complexity away to the matrix library and is equivalent
-        # to the nested for-loop approach (although more idiomatic
-        # and probably faster)
-        sequence.add(self.activation.function((weights * data).sum() + self.biases[0, i]))
-    result = newMatrix[float](sequence)
+    result = ((self.weights * data).sum() + self.biases).apply(self.activation.function, axis= -1)
 
+
+proc cost*(self: Layer, x: Matrix[float], Y: Matrix[float]): float =
+    ## Returns the total cost of this layer
+    
     
diff --git a/src/nn/util/activations.nim b/src/nn/util/activations.nim
index 902d871..85b25aa 100644
--- a/src/nn/util/activations.nim
+++ b/src/nn/util/activations.nim
@@ -24,11 +24,11 @@ func htan*(input: float): float =
 
 
 type Activation* = ref object
-    function*: proc (input: float): float
-    derivative*: proc (x, y: float): float
+    function*: proc (input: float): float {.noSideEffect.}
+    derivative*: proc (x, y: float): float {.noSideEffect.}
 
 
-proc newActivation*(function: proc (input: float): float, derivative: proc (x, y: float): float): Activation =
+proc newActivation*(function: proc (input: float): float {.noSideEffect.}, derivative: proc (x, y: float): float {.noSideEffect.}): Activation =
     ## Creates a new activation object
     new(result)
     result.function = function
diff --git a/src/nn/util/matrix.nim b/src/nn/util/matrix.nim
index 31e39f7..eb05b8c 100644
--- a/src/nn/util/matrix.nim
+++ b/src/nn/util/matrix.nim
@@ -14,7 +14,6 @@
 
 from std/strformat import `&`
 from std/sequtils import zip
-from std/strutils import join
 
 
 type
@@ -29,10 +28,6 @@ type
         ## A zero-copy view into a matrix
         m: Matrix[T]    # The matrix that owns the row we point to
         row: int        # The row in the matrix to which we point to
-        # Even though a MatrixView has no
-        # rows, we keep the same field for
-        # consistency with the Matrix type
-        shape*: tuple[rows, cols: int]
 
 
 proc getSize(shape: tuple[rows, cols: int]): int =
@@ -44,6 +39,10 @@ proc getSize(shape: tuple[rows, cols: int]): int =
     return shape.cols * shape.rows
 
 
+proc shape*[T](self: MatrixView[T]): tuple[rows, cols: int] = 
+    return (0, self.m.shape.cols)
+
+
 proc newMatrix*[T](data: seq[T]): Matrix[T] =
     ## Initializes a new matrix from a given
     ## 1D sequence
@@ -94,8 +93,6 @@ proc zeros*[T](shape: tuple[rows, cols: int], order: MatrixOrder = RowMajor): Ma
 func len*[T](self: Matrix[T]): int {.inline.} = self.data[].len()
 func len*[T](self: MatrixView[T]): int {.inline.} = self.shape.cols
 func raw*[T](self: Matrix[T]): ref seq[T] {.inline.} = self.data
-proc dup*[T](self: Matrix[T]): Matrix[T]
-proc copy*[T](self: Matrix[T]): Matrix[T]
 
 
 func getIndex[T](self: Matrix[T], row, col: int): int =
@@ -113,7 +110,7 @@ proc `[]`*[T](self: Matrix[T], row, col: int): T {.raises: [IndexDefect, ValueEr
     ## column into the matrix
     var idx = self.getIndex(row, col)
     when not defined(release):
-        if idx notin 0..<self.data[].len() + 1:
+        if idx notin 0..<self.data[].len():
             raise newException(IndexDefect, &"index ({row}, {col}) is out of range for matrix of shape ({self.shape.rows}, {self.shape.cols})")
     return self.data[idx]
 
@@ -124,12 +121,11 @@ proc `[]`*[T](self: Matrix[T], row: int): MatrixView[T] {.raises: [IndexDefect,
     ## returned
     var idx = self.getIndex(row, 0)
     when not defined(release):
-        if idx notin 0..<self.data[].len() + 1:
+        if idx notin 0..<self.data[].len():
             raise newException(IndexDefect, &"row {row} is out of range for matrix of shape ({self.shape.rows}, {self.shape.cols})")
     new(result)
     result.m = self
     result.row = row
-    result.shape = (0, self.shape.cols)
 
 
 proc `[]`*[T](self: MatrixView[T], col: int): T {.raises: [IndexDefect, ValueError].} =
@@ -137,7 +133,7 @@ proc `[]`*[T](self: MatrixView[T], col: int): T {.raises: [IndexDefect, ValueErr
     ## the matrix view
     var idx = self.m.getIndex(self.row, col)
     when not defined(release):
-        if idx notin 0..<self.m.data[].len() + 1:
+        if idx notin 0..<self.m.data[].len():
             raise newException(IndexDefect, &"column {col} is out of range for view of shape ({self.shape.rows}, {self.shape.cols})")
     result = self.m.data[idx]
 
@@ -147,7 +143,7 @@ proc `[]=`*[T](self: Matrix[T], row, col: int, val: T) {.raises: [IndexDefect, V
     ## column into the matrix to value val
     var idx = self.getIndex(row, col)
     when not defined(release):
-        if idx notin 0..<self.data[].len() + 1:
+        if idx notin 0..<self.data[].len():
             raise newException(IndexDefect, &"index ({row}, {col}) is out of range for matrix of shape ({self.shape.rows}, {self.shape.cols})")
     self.data[idx] = val
 
@@ -158,47 +154,10 @@ proc `[]=`*[T](self: MatrixView[T], col: int, val: T) {.raises: [IndexDefect, Va
     ## val
     var idx = self.m.getIndex(0, col)
     when not defined(release):
-        if idx notin 0..<self.m.data[].len() + 1:
+        if idx notin 0..<self.m.data[].len():
             raise newException(IndexDefect, &"column {col} is out of range for view of shape ({self.shape.rows}, {self.shape.cols})")
     self.m.data[idx] = val
 
-
-proc `$`*[T](self: Matrix[T]): string =
-    ## Stringifies the matrix
-    var col: int
-    if self.shape.cols == 0:
-        col = self.len()
-    else:
-        col = self.shape.cols
-    if self.shape.rows == 0:
-        return &"""[{self.data[].join(", ")}]"""
-    if self.shape.rows > 1:
-        result &= "["
-    for row in 0..<self.shape.rows:
-        result &= "["
-        for column in 0..<col:
-            result &= $self[row, column]
-            if column < col - 1:
-                result &= ", "
-        if row < self.shape.rows - 1:
-            result &= "], \n"
-            result &= " "
-        else:
-            result &= "]"
-    if self.shape.rows > 1:
-        result &= "]"
-
-
-proc `$`*[T](self: MatrixView[T]): string =
-    ## Stringifies the matrix view
-    result = "["
-    var i = 0
-    while i < self.shape.cols:
-        result &= $self[i]
-        if i < self.shape.cols - 1:
-            result &= ", "
-        inc(i)
-    result &= "]"
     
 
 # Shape management
@@ -208,8 +167,7 @@ proc reshape*[T](self: Matrix[T], shape: tuple[rows, cols: int]): Matrix[T] {.ra
         if shape.getSize() != self.data[].len():
             raise newException(ValueError, &"shape ({shape.rows}, {shape.cols}) is invalid for matrix of length {self.len()}")
     result = self.dup()
-    if shape.rows > 1:
-        self.shape = shape
+    result.shape = shape
 
 
 proc reshape*[T](self: Matrix[T], rows, cols: int): Matrix[T] {.raises: [ValueError].} =
@@ -220,11 +178,8 @@ proc reshape*[T](self: Matrix[T], rows, cols: int): Matrix[T] {.raises: [ValueEr
 proc transpose*[T](self: Matrix[T]): Matrix[T] =
     ## Transposes rows and columns in the given 
     ## matrix. No data copies occur
-    result = self.dup()
-    result.data = self.data
-    discard result.reshape(self.shape.cols, self.shape.rows)
-    if result.shape.rows > 0:
-        result.order = if result.order == RowMajor: ColumnMajor else: RowMajor
+    result = self.reshape(self.shape.cols, self.shape.rows)
+    result.order = if result.order == RowMajor: ColumnMajor else: RowMajor
 
 
 proc flatten*[T](self: Matrix[T]): Matrix[T] =
@@ -236,7 +191,7 @@ proc flatten*[T](self: Matrix[T]): Matrix[T] =
 
 
 # Helpers for fast applying of operations along an axis
-proc apply*[T](self: Matrix[T], op: proc (a, b: T): T, b: T, copy: bool = false, axis: int): Matrix[T] =
+proc apply*[T](self: Matrix[T], op: proc (a, b: T): T {.noSideEffect.}, b: T, copy: bool = false, axis: int): Matrix[T] =
     ## Applies a binary operator to every
     ## element in the given axis of the
     ## given matrix (0 = rows, 1 = columns, 
@@ -247,22 +202,23 @@ proc apply*[T](self: Matrix[T], op: proc (a, b: T): T, b: T, copy: bool = false,
         result = self.copy()
     case axis:
         of 0:
-            # Stores the indeces of the values
-            # we'll delete after we're done. This
-            # is because applying along 
-            var indeces: seq[int] = @[]
+            for r in 0..<self.shape.rows:
+                for c in 1..<self.shape.cols:
+                    result[r, 0] = op(result[r, 0], b)
         of 1:
-            discard
+            for r in 0..<self.shape.rows - 1:
+                for c in 0..self.shape.cols - 1:
+                    result[r, c] = op(result[r, c], b)
         of -1:
-            for i, row in self:
+            for i, row in result:
                 for j, item in row:
-                    self[i, j] = op(j, b)
+                    result[i, j] = op(item, b)
         else:
             raise newException(ValueError, &"axis {axis} is invalid for matrix") 
 
 
-proc apply*[T](self: Matrix[T], op: proc (a: T): T, copy: bool = false, axis: int): Matrix[T] =
-    ## Applies a unary operator to every
+proc apply*[T](self: Matrix[T], op: proc (a: T): T {.noSideEffect.}, copy: bool = false, axis: int): Matrix[T] =
+    ## Applies a binary operator to every
     ## element in the given axis of the
     ## given matrix (0 = rows, 1 = columns, 
     ## -1 = both). No copies occur unless 
@@ -270,17 +226,27 @@ proc apply*[T](self: Matrix[T], op: proc (a: T): T, copy: bool = false, axis: in
     result = self
     if copy:
         result = self.copy()
-    for i, row in self:
-        for j, item in row:
-            self[i, j] = op(j)
+    case axis:
+        of 0:
+            for r in 0..<self.shape.rows:
+                for c in 1..<self.shape.cols:
+                    result[r, 0] = op(result[r, 0])
+        of 1:
+            for r in 0..<self.shape.rows - 1:
+                for c in 0..self.shape.cols - 1:
+                    result[r, c] = op(result[r, c])
+        of -1:
+            for i, row in result:
+                for j, item in row:
+                    result[i, j] = op(item)
+        else:
+            raise newException(ValueError, &"axis {axis} is invalid for matrix") 
 
 
-proc apply*[T](self: MatrixView[T], op: proc (a, b: T): T, b: T, copy: bool = false, axis: int): MatrixView[T] =
+proc apply*[T](self: MatrixView[T], op: proc (a, b: T): T {.noSideEffect.}, b: T, copy: bool = false): MatrixView[T] =
     ## Applies a binary operator to every
-    ## element in the given axis of the
-    ## given matrix view (0 = rows, 1 = columns, 
-    ## -1 = both). No copies occur unless 
-    ## copy equals true
+    ## element in the matrix view. No copies
+    ## occur unless copy equals true
     result = self
     if copy:
         result = self.copy()
@@ -288,12 +254,10 @@ proc apply*[T](self: MatrixView[T], op: proc (a, b: T): T, b: T, copy: bool = fa
         self[i] = op(j, b)
 
 
-proc apply*[T](self: MatrixView[T], op: proc (a: T): T, copy: bool = false, axis: int): MatrixView[T] =
+proc apply*[T](self: MatrixView[T], op: proc (a: T): T {.noSideEffect.}, copy: bool = false): MatrixView[T] =
     ## Applies a unary operator to every
-    ## element in the given axis of the
-    ## given matrix view (0 = rows, 1 = columns, 
-    ## -1 = both). No copies occur unless 
-    ## copy equals true
+    ## element in the matrix view. No copies
+    ## occur unless copy equals true
     result = self
     if copy:
         result = self.copy()
@@ -301,39 +265,6 @@ proc apply*[T](self: MatrixView[T], op: proc (a: T): T, copy: bool = false, axis
         self[i] = op(j)
 
 
-# Operations along an axis
-proc sum*[T](self: Matrix[T], axis: int, copy: bool = true): Matrix[T] =
-    ## Performs the sum of all the elements
-    ## on a given axis in-place (unless copy
-    ## equals true). The output matrix is 
-    ## returned
-    var self = self
-    if copy:
-        self = self.copy()
-    result = self
-    var indeces: seq[int] = @[]
-    case axis:
-        of 1:
-            for r in 0..<self.shape.rows:
-                for c in 1..<self.shape.cols:
-                    self[r, 0] = self[r, 0] + self[r, c]
-                    indeces.add(self.getIndex(r, c))
-        of 0:
-            for r in 0..<self.shape.rows - 1:
-                for c in 0..self.shape.cols - 1:
-                    self[r, c] = self[r, c] + self[r + 1, c]
-                    indeces.add(self.getIndex(r + 1, c))
-        else:
-            when not defined(release):
-                raise newException(ValueError, &"axis {axis} is invalid for matrix")
-            else:
-                discard
-    self.shape.cols = self.len()
-    self.shape.rows = 0
-    for i, index in indeces:
-        self.data[].delete(index - i)
-
-
 proc sum*[T](self: Matrix[T]): T =
     ## Returns the sum of all elements
     ## in the matrix
@@ -341,6 +272,48 @@ proc sum*[T](self: Matrix[T]): T =
         result += e
 
 
+# Operations along an axis
+proc sum*[T](self: Matrix[T], axis: int, copy: bool = true): Matrix[T] =
+    ## Performs the sum of all the elements
+    ## on a given axis in-place (unless copy
+    ## equals true). The output matrix is 
+    ## returned
+    when not defined(release):
+        if axis == 1 and self.shape.rows == 0:
+            raise newException(ValueError, &"axis {axis} is invalid for matrix of dimension 1")
+    var self = self
+    if copy:
+        self = self.copy()
+    result = self
+    var added: int  = 0
+    case axis:
+        of 1:
+            for row in result:
+                inc(added)
+                result.data[].add(row.sum())
+        of 0:
+            var row = 0
+            var value: T
+            for col in 0..<result.shape.cols:
+                while row < result.shape.rows:
+                    value += result[row, col]
+                    inc(row)
+                result.data[].add(value)
+                inc(added)
+                value = T.default
+                row = 0
+        else:
+            when not defined(release):
+                raise newException(ValueError, &"axis {axis} is invalid for matrix")
+            else:
+                discard
+    while result.data[].len() > added:
+        result.data[].delete(0)
+    result.shape.rows = 0
+    result.shape.cols = added
+    result.order = RowMajor
+
+
 proc sum*[T](self: MatrixView[T]): T =
     ## Returns the sum of all elements
     ## in the matrix view
@@ -350,40 +323,6 @@ proc sum*[T](self: MatrixView[T]): T =
         inc(i)
 
 
-proc sub*[T](self: Matrix[T], axis: int, copy: bool = true): Matrix[T] =
-    var self = self
-    if copy:
-        self = self.copy()
-    result = self
-    var indeces: seq[int] = @[]
-    case axis:
-        of 1:
-            for r in 0..<self.shape.rows:
-                for c in 1..<self.shape.cols:
-                    self[r, 0] = self[r, 0] - self[r, c]
-                    indeces.add(self.getIndex(r, c))
-        of 0:
-            for r in 0..<self.shape.rows - 1:
-                for c in 0..self.shape.cols - 1:
-                    self[r, c] = self[r, c] - self[r + 1, c]
-                    indeces.add(self.getIndex(r + 1, c))
-        else:
-            when not defined(release):
-                raise newException(ValueError, &"axis {axis} is invalid for matrix")
-            else:
-                discard
-    self.shape.cols = 0
-    self.shape.rows = 1
-    self.order = RowMajor
-    for i, index in indeces:
-        self.data[].delete(index - i)
-
-
-proc sub*[T](self: Matrix[T]): T =
-    for e in self.data[]:
-        result -= e
-
-
 proc copy*[T](self: Matrix[T]): Matrix[T] =
     ## Creates a new copy of the given matrix
     ## (copies the underlying data!)
@@ -428,87 +367,159 @@ proc dup*[T](self: MatrixView[T]): MatrixView[T] =
 
 # Wrappers because builtins are not
 # procvars
-proc add*[T](a, b: T): T = a + b
-proc sub*[T](a, b: T): T = a - b
-proc mul*[T](a, b: T): T = a * b
-proc divide*[T](a, b: T): T = a / b
-proc neg*[T](a: T): T = -a
+func add*[T](a, b: T): T = a + b
+func sub*[T](a, b: T): T = a - b
+func mul*[T](a, b: T): T = a * b
+func divide*[T](a, b: T): T = a / b
+func neg*[T](a: T): T = -a
 
 # Warning: These *all* perform copies of the underlying matrix!
-template `+`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(add, b)
-template `+`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(add, a)
+proc `+`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(add, b, axis= -1)
+proc `+`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(add, a, axis= -1)
 
-template `-`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(sub, b)
-template `-`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(sub, a)
-template `-`*[T](a: Matrix[T]): Matrix[T] = a.copy().apply(neg, a)
+proc `-`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(sub, b, axis= -1)
+proc `-`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(sub, a, axis= -1)
+proc `-`*[T](a: Matrix[T]): Matrix[T] = a.copy().apply(neg, a, axis= -1)
 
-template `*`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(mul, b)
-template `*`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(mul, a)
+proc `*`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(mul, b, axis = -1)
+proc `*`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(mul, a, axis= -1)
 
-template `/`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(divide, b)
-template `/`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(divide, a)
+proc `/`*[T](a: Matrix[T], b: T): Matrix[T] = a.copy().apply(divide, b, axis= -1)
+proc `/`*[T](a: T, b: Matrix[T]): Matrix[T] = b.copy().apply(divide, a, axis= -1)
 
 
+# matrix/matrix operations. They produce a new matrix with the
+# result of the operation
+
+proc `+`*[T](a, b: MatrixView[T]): Matrix[T] =
+    ## Performs the vector sum of the 
+    ## given matrix views and returns a new
+    ## vector with the result
+    when not defined(release):
+        if a.shape.cols != b.shape.cols:  # Basically if their length is different
+            raise newException(ValueError, &"incompatible argument shapes for addition")
+    new(result)
+    new(result.data)
+    result.shape = a.shape
+    result.order = RowMajor
+    result.data[] = newSeqOfCap[T](result.shape.getSize())
+    for i in 0..<a.shape.cols:
+        result.data[].add(a[i] + b[i])
+
+
+proc `+`*[T](a, b: Matrix[T]): Matrix[T] {.raises: [ValueError].} =
+    when not defined(release):
+        if a.shape.rows > 0 and b.shape.rows > 0 and a.shape != b.shape:
+            raise newException(ValueError, &"incompatible argument shapes for addition")
+        elif (a.shape.rows == 0 or b.shape.rows == 0) and a.shape.cols != b.shape.cols:
+            raise newException(ValueError, &"incompatible argument shapes for addition")
+    if a.shape.rows == 0 and b.shape.rows == 0:
+        return a[0] + b[0]
+    new(result)
+    new(result.data)
+    result.data[] = newSeqOfCap[T](result.shape.getSize())
+    result.shape = a.shape
+    result.order = RowMajor
+    if result.shape.rows > 1:
+        for row in 0..<result.shape.rows:
+            for m in a[row] + b[row]:
+                for element in m:
+                    result.data[].add(element)
+    else:
+        result = a[0] + b[0]
+
+
+proc `*`*[T](a, b: MatrixView[T]): Matrix[T] =
+    ## Performs the vector product of the 
+    ## given matrix views and returns a new
+    ## vector with the result
+    when not defined(release):
+        if a.shape.cols != b.shape.cols:  # Basically if their length is different
+            raise newException(ValueError, &"incompatible argument shapes for multiplication")
+    new(result)
+    new(result.data)
+    result.shape = a.shape
+    result.order = RowMajor
+    result.data[] = newSeqOfCap[T](result.shape.getSize())
+    for i in 0..<a.shape.cols:
+        result.data[].add(a[i] * b[i])
+
+
+proc `*`*[T](a, b: Matrix[T]): Matrix[T] {.raises: [ValueError].} =
+    when not defined(release):
+        if a.shape.rows > 0 and b.shape.rows > 0 and a.shape.cols != b.shape.rows:
+            raise newException(ValueError, &"incompatible argument shapes for multiplication")
+        elif (a.shape.rows == 0 or b.shape.rows == 0) and a.shape.cols != b.shape.cols:
+            raise newException(ValueError, &"incompatible argument shapes for multiplication")
+    new(result)
+    new(result.data)
+    result.shape = (a.shape.rows, b.shape.cols)
+    result.order = RowMajor
+    result.data[] = newSeqOfCap[T](result.shape.getSize())
+    if result.shape.rows > 1:
+        if a.shape.rows == b.shape.rows:
+            for row in 0..<result.shape.rows:
+                for m in a[row] * b[row]:
+                    for element in m:
+                        result.data[].add(element)
+        elif b.shape.rows < a.shape.rows:
+            for r1 in b:
+                for r2 in a:
+                    for m in r1 * r2:
+                        for element in m:
+                            result.data[].add(element)
+    else:
+        result = a[0] * b[0]
+
+# Comparison operators. They produce a new matrix of the same
+# shape as the input(s) and containing boolean values (the result of
+# the comparison element-wise). Useful for use in where()
+
 # matrix/scalar comparisons
 proc `==`*[T](a: Matrix[T], b: T): Matrix[bool] = 
     new(result)
     new(result.data)
     result.shape = a.shape
-    result.data[] = newSeqOfCap[bool](result.shape.rows * result.shape.cols)
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
     for e in a.data[]:
         result.data[].add(e == b)
 
 
-# matrix/matrix operations. They produce a new matrix with the
-# result of the operation
-
-proc `+`*[T](a, b: Matrix[T]): Matrix[T] {.raises: [ValueError].} =
-    when not defined(release):
-        if a.shape != b.shape:
-            raise newException(ValueError, "can't add matrices of different shapes")
-        if a.order != b.order:
-            raise newException(ValueError, "can't add matrices with different ordering")
+proc `<`*[T](a: Matrix[T], b: T): Matrix[bool] = 
     new(result)
     new(result.data)
-    result.data[] = newSeqOfCap[T](result.shape.rows * result.shape.cols)
     result.shape = a.shape
-    result.order = RowMajor
-    for (e, k) in zip(a.data[], b.data[]):
-        result.data[].add(e + k)
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
+    for e in a.data[]:
+        result.data[].add(e < b)
 
 
-proc `*`*[T](a: MatrixView[T], b: Matrix[T]): Matrix[T] {.raises: [ValueError].} =
-    when not defined(release):
-        echo a
-        echo b
-        if a.shape.cols != b.shape.rows:
-            raise newException(ValueError, &"incompatible argument shapes for multiplication")
-        if a.m.order != b.order:
-            raise newException(ValueError, "can't multiply matrices with different ordering")
+proc `>`*[T](a: Matrix[T], b: T): Matrix[bool] = 
     new(result)
     new(result.data)
-    result.data[] = newSeqOfCap[T](result.shape.rows * result.shape.cols)
-    for i in 0..<a.shape.cols:
-        result.data[].add(a[i] * b[0, i])
     result.shape = a.shape
-    result.order = RowMajor
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
+    for e in a.data[]:
+        result.data[].add(e > b)
 
 
-proc `*`*[T](a, b: Matrix[T]): Matrix[T] {.raises: [ValueError].} =
-    when not defined(release):
-        if a.shape.cols != b.shape.rows:
-            raise newException(ValueError, &"incompatible argument shapes for multiplication")
-        if a.order != b.order:
-            raise newException(ValueError, "can't multiply matrices with different ordering")
+proc `<=`*[T](a: Matrix[T], b: T): Matrix[bool] = 
     new(result)
     new(result.data)
-    result.shape = (a.shape.rows, b.shape.cols)
-    result.order = RowMajor
-    result.data[] = newSeqOfCap[T](result.shape.rows * result.shape.cols)
-    for (e, k) in zip(a.data[], b.data[]):
-        result.data[].add(e * k)
+    result.shape = a.shape
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
+    for e in a.data[]:
+        result.data[].add(e <= b)
+
+
+proc `>=`*[T](a: Matrix[T], b: T): Matrix[bool] = 
+    new(result)
+    new(result.data)
+    result.shape = a.shape
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
+    for e in a.data[]:
+        result.data[].add(e >= b)
 
-# Comparison operators. They produce a new matrix with boolean values
 
 proc `==`*[T](a, b: Matrix[T]): Matrix[bool] {.raises: [ValueError].} =
     when not defined(release):
@@ -518,7 +529,7 @@ proc `==`*[T](a, b: Matrix[T]): Matrix[bool] {.raises: [ValueError].} =
     new(result.data)
     result.shape = a.shape
     result.order = RowMajor
-    result.data[] = newSeqOfCap[bool](result.shape.rows * result.shape.cols)
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
     for r in 0..<a.shape.rows:
         for c in 0..<a.shape.cols:
             result.data[].add(a[r, c] == b[r, c])
@@ -532,7 +543,7 @@ proc `>`*[T](a, b: Matrix[T]): Matrix[bool] {.raises: [ValueError].} =
     new(result.data)
     result.shape = a.shape
     result.order = RowMajor
-    result.data[] = newSeqOfCap[bool](result.shape.rows * result.shape.cols)
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
     for r in 0..<a.shape.rows:
         for c in 0..<a.shape.cols:
             result.data[].add(a[r, c] > b[r, c])
@@ -546,7 +557,7 @@ proc `>=`*[T](a, b: Matrix[T]): Matrix[bool] {.raises: [ValueError].} =
     new(result.data)
     result.shape = a.shape
     result.order = RowMajor
-    result.data[] = newSeqOfCap[bool](result.shape.rows * result.shape.cols)
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
     for r in 0..<a.shape.rows:
         for c in 0..<a.shape.cols:
             result.data[].add(a[r, c] >= b[r, c])
@@ -560,7 +571,7 @@ proc `<=`*[T](a, b: Matrix[T]): Matrix[bool] {.raises: [ValueError].} =
     new(result.data)
     result.shape = a.shape
     result.order = RowMajor
-    result.data[] = newSeqOfCap[bool](result.shape.rows * result.shape.cols)
+    result.data[] = newSeqOfCap[bool](result.shape.getSize())
     for r in 0..<a.shape.rows:
         for c in 0..<a.shape.cols:
             result.data[].add(a[r, c] <= b[r, c])
@@ -593,9 +604,9 @@ proc toRowMajor*[T](self: Matrix[T]): Matrix[T] =
     self.data[] = @[]
     var idx = 0
     var col = 0
-    while col < result.shape.cols:
-        result.data[].add(orig[idx])
-        idx += result.shape.cols
+    while col < self.shape.cols:
+        self.data[].add(orig[idx])
+        idx += self.shape.cols
         if idx > orig.high():
             inc(col)
             idx = col
@@ -605,6 +616,7 @@ proc toRowMajor*[T](self: Matrix[T]): Matrix[T] =
 proc toColumnMajor*[T](self: Matrix[T]): Matrix[T] =
     ## Converts a row-major matrix to a
     ## column-major one
+    new(result)
     if self.order == ColumnMajor:
         return
     self.order = ColumnMajor
@@ -612,9 +624,9 @@ proc toColumnMajor*[T](self: Matrix[T]): Matrix[T] =
     self.data[] = @[]
     var idx = 0
     var col = 0
-    while col < result.shape.cols:
-        result.data[].add(orig[idx])
-        idx += result.shape.cols
+    while col < self.shape.cols:
+        self.data[].add(orig[idx])
+        idx += self.shape.cols
         if idx > orig.high():
             inc(col)
             idx = col
@@ -626,6 +638,8 @@ proc toColumnMajor*[T](self: Matrix[T]): Matrix[T] =
 iterator items*[T](self: Matrix[T]): MatrixView[T] =
     for row in 0..<self.shape.rows:
         yield self[row]
+    if self.shape.rows == 0:
+        yield self[0]
 
 
 iterator items*[T](self: MatrixView[T]): T =
@@ -637,28 +651,94 @@ iterator pairs*[T](self: Matrix[T]): tuple[i: int, val: MatrixView[T]] =
     var i = 0
     for row in self:
         yield (i, row)
+        inc(i)
 
 
 iterator pairs*[T](self: MatrixView[T]): tuple[i: int, val: T] =
     var i = 0
     for col in self:
         yield (i, col)
+        inc(i)
+
+
+proc `$`*[T](self: MatrixView[T]): string =
+    ## Stringifies the matrix view
+    result = "["
+    for j, e in self:
+        result &= $e
+        if j < self.shape.cols - 1:
+            result &= ", "
+    result &= "]"
+
+
+proc `$`*[T](self: Matrix[T]): string =
+    ## Stringifies the matrix
+    if self.shape.rows == 0:
+        return $self[0]
+    result &= "["
+    for i, row in self:
+        result &= "["
+        for j, e in row:
+            result &= $e
+            if j < self.shape.cols - 1:
+                result &= ", "
+        if i < self.shape.rows - 1:
+            result &= "], \n"
+            result &= " "
+        else:
+            result &= "]"
+    result &= "]"
 
 
 proc dot*[T](self, other: Matrix[T]): Matrix[T] =
     ## Computes the dot product of the two
     ## input matrices
-    
+    when not defined(release):
+        if a.shape.cols != b.shape.rows:
+            raise newException(ValueError, &"incompatible argument shapes for dot product")
+    # TODO
+
+
+proc where*[T](cond: Matrix[bool], x, y: Matrix[T]): Matrix[T] =
+    ## Behaves like numpy.where()
+    when not defined(release):
+        if not (x.shape == y.shape and y.shape == cond.shape):
+            raise newException(ValueError, &"all inputs must be of equal shape for where()")
+    result = x.copy()
+    var 
+        row = 0
+        col = 0
+    if cond.shape.rows == 0:
+        while col < cond.shape.cols:
+            if not cond[0, col]:
+                result[0, col] = y[0, col]
+            inc(col)
+    while row < cond.shape.rows:
+        while col < cond.shape.cols:
+            if not cond[row, col]:
+                result[row, col] = y[row, col]
+            inc(col)
+        inc(row)
+        col = 0
 
 when isMainModule:
-    # var m = newMatrix[int](@[@[1, 2, 3], @[4, 5, 6]])
-    # var k = m.transpose()
-    # assert all(m.transpose() == k), "transpose mismatch"
-    # assert k.sum() == m.sum(), "element sum mismatch"
-    # assert k.sub() == m.sub(), "element sub mismatch"
-    # assert all(k.sum(axis=1) == m.sum(axis=0)), "sum over axis mismatch"
-    # assert all(k.sum(axis=0) == m.sum(axis=1)), "sum over axis mismatch" 
-    # assert all(k.sub(axis=1) == m.sub(axis=0)), "sub over axis mismatch"
-    # assert all(k.sub(axis=0) == m.sub(axis=1)), "sub over axis mismatch"
+    import math
+
+    proc pow(a, b: int): int =
+        return a ^ b
+
+    var m = newMatrix[int](@[@[1, 2, 3], @[4, 5, 6]])
+    var k = m.transpose()
+    assert all(m.transpose() == k), "transpose mismatch"
+    assert k.sum() == m.sum(), "element sum mismatch"
+    assert all(k.sum(axis=1) == m.sum(axis=0)), "sum over axis mismatch"
+    assert all(k.sum(axis=0) == m.sum(axis=1)), "sum over axis mismatch" 
     var y = newMatrix[int](@[1, 2, 3, 4])
-    echo y.sum(axis=0)
\ No newline at end of file
+    assert y.sum() == 10, "element sum mismatch"
+    assert (y + y).sum() == 20, "matrix sum mismatch"
+    assert all(m + m == m * 2), "m + m != m * 2"
+    var z = newMatrix[int](@[1, 2, 3])
+    assert (m * z).sum() == 46, "matrix multiplication mismatch"
+    assert all(z * z == z.apply(pow, 2, axis = -1, copy=true)), "matrix multiplication mismatch"
+    var x = newMatrix[int](@[0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
+    assert (x < 5).where(x, x * 10).sum() == 360, "where mismatch"
\ No newline at end of file