Add my linear algebra library

6 files changed, 427 insertions, 0 deletions

diff --git a/lalge/README.md b/lalge/README.md
new file mode 100644
index 0000000..640654d
--- /dev/null
+++ b/lalge/README.md
@@ -0,0 +1,46 @@
+# lalge (name tbd)
+
+A functional, dead simple linear algebra library written in Nim. Mostly for learning purposes.
+
+This is also serving as my finals review, haha.
+
+If you're looking for an actually robust linear algebra library for Nim, check out [Neo](https://github.com/andreaferretti/neo).
+
+## What's this language?
+
+Nim is a compiled whitespace-sensitive systems language with a strong type system and interesting memory management capabilities.
+
+It doesn't quite beat Python's usability or Rust's features, but I like it quite a lot, and think it's a good complement to other languages I use.
+
+Nim's [UFCS](https://en.wikipedia.org/wiki/Uniform_Function_Call_Syntax) (TL;DR: `a.add(b) == add(a, b)`) and operator overloading support is also a major draw for this type of project.
+
+## Architecture
+
+lalge uses Nim's `seqs` - dynamic arrays, equivalent to Rust's `Vec`, Go's `slice`, or Java's `ArrayList` - to implement matrices. This may have some performance overhead. On the other hand, it has the benefit of making the code much simpler.
+
+This also has the benefit of being able to declare matrices as such:
+```nim
+example: Matrix = @[
+  @[1, 0, 0],
+  @[0, 1, 0],
+  @[0, 0, 1],
+]
+```
+
+## Goals
+- [ ] Cover the basic operations associated with the MATH 221 / typical introductory linear algebra curriculum
+  - [x] Addition / Subtraction
+  - [x] (Scalar / Vector) Multiplication
+  - [x] Inverse
+  - [x] Transpose
+  - [ ] Eigenvalues / Eigenvectors
+  - [ ] everything that i am forgetting
+- [ ] Helpful characteristic functions
+  - [ ] Dimension
+  - [ ] Rank / Nullity
+  - [ ] Spans
+  - [ ] Determinant
+- [x] Support arbitrary-sized matrices with internal seq representation
+- [x] Support specialized "vectors", as opposed to Nx1 matrices
+
+After I get to a reasonable point of completeness, I'd love to go through and compare my library to Eigen (generally praised as a Very Good library)
diff --git a/lalge/src/matrices.nim b/lalge/src/matrices.nim
new file mode 100644
index 0000000..8e3e680
--- /dev/null
+++ b/lalge/src/matrices.nim
@@ -0,0 +1,171 @@
+import std/math, std/sugar
+import types
+
+# A few notes about Nim:
+
+# "result" lets us operate on the return type
+# without (necessarily) constructing it.
+# this is very useful, and as such i use it a lot.
+
+# nim also supplies many functional tools which
+# i also made liberal use of, occasionally to the
+# detriment of readability. sorry.
+
+# the "*" operator on procs is the public marker.
+# oh, and functions with "side effects" are called "procs".
+
+# TODO: implement a better error handling system
+const not_square = "The matrix is not square"
+const mismatch = "Mismatched rows and columns"
+
+## Overrides the echo() proc to provide more human-readable output
+proc echo*(A: Matrix) =
+  echo "@["
+  for row in A:
+    echo "  ", row
+  echo "]"
+
+## Generate a new, filled XxX Matrix
+func gen*(rows, columns: Natural, element: float = 0.0): Matrix =
+  result.setLen(rows)
+  for i in 0 ..< rows:
+    result[i].setLen(columns)
+    for j in 0 ..< columns:
+      result[i][j] = element
+
+## Returns the number of rows in a matrix
+func rows*(A: Matrix): int =
+  return A.len()
+
+## Returns the number of columns in a matrix
+func cols*(A: Matrix): int =
+  if A.rows() == 0:
+    return 0
+  else:
+    return A[0].len()
+
+## Simple rows iterator for uniformity
+iterator rows*(A: Matrix): RowVector =
+  for row in A:
+    yield row
+
+## Simple column iterator for matrices
+iterator cols*(A: Matrix): Vector =
+  for i in 0 ..< A.cols():
+    var row: Vector
+    for j in 0 ..< A.rows():
+      row.add(A[i][j])
+    yield row
+
+## Apply an arbitrary function to every element of a matrix
+func map*(A: Matrix, oper: proc(i, j: int): float): Matrix =
+  result = gen(A.rows(), A.cols())
+  for i in 0 ..< A.rows():
+    for j in 0 ..< A.cols():
+      # Note that the position is _actively_ provided
+      result[i][j] = oper(i, j)
+
+# TODO: check if these meet assert requirements
+
+## Matrix addition
+func `+`*(A, B: Matrix): Matrix =
+  return map(A, (i, j) => A[i][j] + B[i][j])
+
+## Matrix subtraction
+func `-`*(A, B: Matrix): Matrix =
+  return map(A, (i, j) => A[i][j] - B[i][j])
+
+## Scalar-Matrix multiplication
+func `*`*(a: float, B: Matrix): Matrix =
+  return map(B, (i, j) => a * B[i][j])
+
+## Scalar-Matrix multiplication
+# func `*`*(A: Matrix, b: int): Matrix =
+#   return b * A
+
+## Matrix multiplication
+func `*`*(A, B: Matrix): Matrix =
+  assert A.rows() == B.cols(), mismatch
+  result = gen(A.rows(), B.cols())
+  for i in 0 ..< A.rows():
+    for j in 0 ..< B.cols():
+      for k in 0 ..< B.rows():
+        result[i][j] += A[i][k] * B[k][j]
+
+## Absolute value of a matrix
+func abs*(A: Matrix): Matrix =
+  return map(A, (i, j) => abs(A[i][j]))
+
+## Simple row function for uniformity
+func row*(A: Matrix, r: int): RowVector =
+  return A[r]
+
+## Returns a "column vector" of a matrix as a Xx1 Matrix
+# func col*(A: Matrix, j: int): Vector =
+#   result = gen(A.rows(), 1)
+#   for i in 0 ..< A.rows():
+#     result[i] = @[A[i][j]]
+
+## Alternate implementation of col, returns a Vector
+func col*(A: Matrix, c: int): Vector =
+  for row in A:
+    result.add(@[row[c]])
+
+## Produce a vector of the diagonal entries of a matrix
+func diag*(A: Matrix): Vector =
+  assert A.rows() == A.cols(), not_square
+  for i in 0 ..< A.rows():
+    result.add(A[i][i])
+
+## Transpose the rows and columns of a matrix
+func transpose*(A: Matrix): Matrix =
+  result = gen(A.cols(), A.rows())
+  return map(result, (i, j) => (A[j][i]))
+
+# func transpose*(a: Vector): Matrix =
+#   result = gen(a.len(), 1)
+#   return map(result, (i, j) => (a[j]))
+
+## Generate an arbitary sized identity matrix
+func identity*(size: Natural): Matrix =
+  return map(gen(size, size), (i, j) => (if i==j: 1.0 else: 0.0))
+
+## Oft-used identity matrices
+let
+  I1*: Matrix = identity(1)
+  I2*: Matrix = identity(2)
+  I3*: Matrix = identity(3)
+  I4*: Matrix = identity(4)
+  I5*: Matrix = identity(5)
+
+## Calculates the determinant of a matrix through Laplace expansion
+func det*(A: Matrix): float =
+  assert A.rows() == A.cols(), not_square
+  # Shortcut by formula
+  if A.rows() == 2:
+    return (A[0][0]*A[1][1]) - (A[0][1]*A[1][0])
+  # Actual factual formula
+  for i in 0 ..< A.rows():
+    var sub: Matrix
+    for a in 1 ..< A.rows():
+      var row: RowVector
+      for b in 0 ..< A.cols():
+        if b == i: continue
+        row.add(A[a][b])
+      sub.add(row)
+    result += (-1.0)^i * A[0][i] * det(sub)
+
+# func spans(A, B: Matrix): bool =
+#   return false
+# func subspace(a, b: Matrix): bool =
+#   return false
+
+# Matrix Properties
+
+## Calculate the dimension of a matrix by Gaussian Elimination
+# func dim(a: Matrix): int =
+#   return 1
+# func rank(a: Matrix): int =
+#   return 1
+# func nullity(a: Matrix): int =
+#   return 1
diff --git a/lalge/src/types.nim b/lalge/src/types.nim
new file mode 100644
index 0000000..321042a
--- /dev/null
+++ b/lalge/src/types.nim
@@ -0,0 +1,8 @@
+type
+  # Any generic Vector is assumed to be a column vector.
+  Vector* = seq[float]
+  RowVector* = Vector # can this be distinct?
+  Matrix* = seq[RowVector]
+
+type
+  Radian* = float
diff --git a/lalge/src/vectors.nim b/lalge/src/vectors.nim
new file mode 100644
index 0000000..a9880d9
--- /dev/null
+++ b/lalge/src/vectors.nim
@@ -0,0 +1,57 @@
+import std/math, std/sequtils, std/sugar
+import types
+
+const dim_mismatch = "Dimensional mismatch - check your vectors"
+
+## Vector addition
+func `+`*(a, b: Vector): Vector =
+  assert a.len() == b.len(), dim_mismatch
+  result.setLen(a.len)
+  for i in 0 ..< a.len():
+    result[i] = (a[i] + b[i])
+
+## Vector subtraction
+func `-`*(a, b: Vector): Vector =
+  assert a.len() == b.len(), dim_mismatch
+  result.setLen(a.len)
+  for i in 0 ..< a.len():
+    result[i] = (a[i] - b[i])
+
+## Vector dot product
+func `*`*(a, b: Vector): float =
+  assert a.len() == b.len(), dim_mismatch
+  for i in 0 ..< a.len():
+    result += a[i] * b[i]
+
+func dot*(a, b: Vector): float =
+  return a * b
+
+## Scalar-Vector multiplication
+func `*`*(a: float, b: Vector): Vector =
+  return map(b, (x) => (a*x))
+
+## Produce the length (in space) of a vector
+func length*(a: Vector): float =
+  return sqrt(a * a)
+
+## Produce the number of elements in a vector
+func size*(a: Vector): int =
+  return len(a)
+
+## Returns whether two vectors are orthagonal
+func ortho*(a, b: Vector): bool =
+  return a * b == 0
+
+## Produce the angle between two vectors, in radians
+func angle*(a, b: Vector): Radian =
+  return arccos((a * b) / (a.length * b.length))
+
+## Produce the cross product between two 3D vectors
+func cross*(a, b: Vector): Vector =
+  assert a.len() == 3, "The first vector is not three-dimensional"
+  assert b.len() == 3, "The second vector is not three-dimensional"
+  return @[
+    (a[1]*b[2]) - (b[1]*a[2]),
+    (a[2]*b[0]) - (b[2]*a[0]),
+    (a[0]*b[1]) - (b[0]*a[1])
+  ]
diff --git a/lalge/tests/nim.cfg b/lalge/tests/nim.cfg
new file mode 100644
index 0000000..85bf6c4
--- /dev/null
+++ b/lalge/tests/nim.cfg
@@ -0,0 +1 @@
+--path:"../src/"
diff --git a/lalge/tests/tests.nim b/lalge/tests/tests.nim
new file mode 100644
index 0000000..e798c39
--- /dev/null
+++ b/lalge/tests/tests.nim
@@ -0,0 +1,144 @@
+import std/math
+import types, vectors, matrices
+
+let A: Matrix = @[
+  @[10.0, 2.0, 3.0, 5.0, 6.0],
+  @[8.0, 7.0, 6.0, 4.0, 3.0],
+  @[4.0, 6.0, 0.0, 9.0, 0.0],
+  @[6.0, 7.0, 9.0, 3.0, 9.0],
+  @[3.0, 0.0, 7.0, 9.0, 9.0],
+]
+let B: Matrix = @[
+  @[1.0, 2.0, 3.0, 4.0, 5.0],
+  @[5.0, 3.0, 1.0, 22.0, 3.0],
+  @[5.0, 21.0, 4.0, 6.0, 3.0],
+  @[12.0, 1.0, 5.0, 0.0, 9.0],
+  @[6.0, 7.0, 1.0, 3.0, 5.0],
+]
+
+let L: Matrix = @[
+  @[1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
+  @[-1.0, 1.0, 0.0, 0.0, 0.0, 0.0],
+  @[0.0, -1.0, 1.0, 0.0, 0.0, 0.0],
+  @[0.0, 0.0, -1.0, 1.0, 0.0, 0.0],
+  @[0.0, 0.0, 0.0, -1.0, 1.0, 0.0],
+  @[0.0, 0.0, 0.0, 0.0, -1.0, 1.0],
+]
+let U: Matrix = @[
+  @[1.0, -1.0, 0.0, 0.0, 0.0, 0.0],
+  @[0.0, 1.0, -1.0, 0.0, 0.0, 0.0],
+  @[0.0, 0.0, 1.0, -1.0, 0.0, 0.0],
+  @[0.0, 0.0, 0.0, 1.0, -1.0, 0.0],
+  @[0.0, 0.0, 0.0, 0.0, 1.0, -1.0],
+  @[0.0, 0.0, 0.0, 0.0, 0.0, 1.0],
+]
+
+let
+  P13: Matrix = @[
+    @[0.0, 0.0, 1.0],
+    @[0.0, 1.0, 0.0],
+    @[1.0, 0.0, 0.0],
+  ]
+  P23: Matrix = @[
+    @[1.0, 0.0, 0.0],
+    @[0.0, 0.0, 1.0],
+    @[0.0, 1.0, 0.0],
+  ]
+
+let foo = @[
+  @[1.0, 1.0, 1.0],
+  @[3.0, 4.0, 5.0],
+  @[40.0, 51.0, 12.0]
+]
+
+let eliminationMatrix: Matrix = @[
+  @[1.0, 0.0, 0.0],
+  @[-2.0, 1.0, 0.0],
+  @[0.0, 0.0, 1.0]
+]
+
+let
+  ones: Matrix = @[
+    @[1.0, 2.0, 3.0],
+    @[4.0, 5.0, 6.0]
+  ]
+  twos: Matrix = @[
+    @[7.0, 8.0],
+    @[9.0, 10.0],
+    @[11.0, 12.0]
+  ]
+
+# Addition and subtraction
+assert I2 + I2 == 2 * I2
+assert I3-I3 == gen(3, 3, 0.0)
+
+# Num. of rows / columns
+assert rows(I2) == 2
+assert I2.rows() == 2
+assert I2.rows == 2
+assert cols(I3) == 3
+assert I3.cols() == 3
+assert I3.cols == 3
+
+# Identity Matrices
+assert I1 == @[
+  @[1.0],
+]
+assert I2 == @[
+  @[1.0, 0.0],
+  @[0.0, 1.0],
+]
+assert I3 == @[
+  @[1.0, 0.0, 0.0],
+  @[0.0, 1.0, 0.0],
+  @[0.0, 0.0, 1.0],
+]
+assert I4 == @[
+  @[1.0, 0.0, 0.0, 0.0],
+  @[0.0, 1.0, 0.0, 0.0],
+  @[0.0, 0.0, 1.0, 0.0],
+  @[0.0, 0.0, 0.0, 1.0],
+]
+assert I5 == @[
+  @[1.0, 0.0, 0.0, 0.0, 0.0],
+  @[0.0, 1.0, 0.0, 0.0, 0.0],
+  @[0.0, 0.0, 1.0, 0.0, 0.0],
+  @[0.0, 0.0, 0.0, 1.0, 0.0],
+  @[0.0, 0.0, 0.0, 0.0, 1.0],
+]
+
+assert det(@[
+  @[2.0, -1.0, 0.0, -1.0],
+  @[-1.0, 2.0, -1.0, 0.0],
+  @[0.0, -1.0, 2.0, -1.0],
+  @[-1.0, 0.0, -1.0, 2.0],
+]) == 0
+
+assert det(A * B) == (det(A) * det(B))
+assert det(I3) == 1
+assert det(2*I4) == 16
+
+# Vectors
+
+let
+  alice: Vector = @[1.0, 1.0, 1.0]
+  bob = @[3.0, 4.0, 5.0]
+  carol = @[40.0, 51.0, 12.0]
+  dan: Vector = @[10.0, 20.0, 30.0, 40.0]
+  eve = @[2.0, 3.0, 4.0]
+
+assert alice + bob == bob + alice
+assert alice + carol == carol + alice
+assert carol - eve != eve - carol
+assert alice.length() == sqrt(3.0)
+
+# echo transpose(ones)
+# echo P13 * 8
+# echo P13 + P23
+# echo ones * twos
+# assert col(eliminationMatrix, 0) == col(eliminationMatrix, 0)
+# echo identity(0)
+# echo identity(2)
+# echo (1/0)
+# echo L*U
+# echo U*L