Many fixes to matrix library and minor changes

This commit is contained in:
Mattia Giambirtone 2022-12-22 14:55:26 +01:00 committed by Nocturn9x
parent dc44b65e94
commit d6e5e148aa
4 changed files with 309 additions and 232 deletions

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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)
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"