From cd6ca6f3ec86a144d348b52e64695fad650da582 Mon Sep 17 00:00:00 2001 From: Walter Oggioni Date: Mon, 31 Dec 2018 01:36:45 +0100 Subject: [PATCH] Several improvements: - added some meaningful unit tests - splitted classes into heap-allocated (starting with 'H') and stack-allocated (starting with 'S'), unfortunately there is a lot of duplicated code but I am still unable to find an elegant solution to use the smae code to deal with both the stack-allocated class and the heap-allocated one --- matrix.nim | 329 ---------------------------------- matrix_test.nim | 183 ------------------- mmath.nimble | 13 ++ src/mmath/hmatrix.nim | 398 +++++++++++++++++++++++++++++++++++++++++ src/mmath/hvector.nim | 48 +++++ src/mmath/pivot.nim | 30 ++++ src/mmath/smatrix.nim | 370 ++++++++++++++++++++++++++++++++++++++ src/mmath/svector.nim | 44 +++++ tests/config.nims | 1 + tests/test_hmatrix.nim | 114 ++++++++++++ tests/test_smatrix.nim | 117 ++++++++++++ vector.nim | 72 -------- 12 files changed, 1135 insertions(+), 584 deletions(-) delete mode 100644 matrix.nim delete mode 100644 matrix_test.nim create mode 100644 mmath.nimble create mode 100644 src/mmath/hmatrix.nim create mode 100644 src/mmath/hvector.nim create mode 100644 src/mmath/pivot.nim create mode 100644 src/mmath/smatrix.nim create mode 100644 src/mmath/svector.nim create mode 100644 tests/config.nims create mode 100644 tests/test_hmatrix.nim create mode 100644 tests/test_smatrix.nim delete mode 100644 vector.nim diff --git a/matrix.nim b/matrix.nim deleted file mode 100644 index d793830..0000000 --- a/matrix.nim +++ /dev/null @@ -1,329 +0,0 @@ -from sequtils import newSeqWith, map -from utils import `...`, `-->` -from vector import Vector, newVector, createVector -from options import Option, none -import future -from random import randomize, random - -type - Matrix*[T] = object - rows, columns : int - data : seq[T] - SMatrix[W, H: static[int], T] = - array[1..W, array[1..H, T]] - Pivot* = ref object - data* : Vector[int] - permutations* : int - SingularMatrixError* = object of ValueError - SizeError* = object of ValueError - MatrixRef[T] = ref Matrix[T] - -proc `[]`(p : Pivot, index : int) : int = p.data[index] -proc `[]=`(p : var Pivot, index : int, value : int) = p.data[index] = value -proc len(p : Pivot) : int = p.data.len - -proc newPivot[T](size : int) : Pivot = - result = new (Pivot) - result = Pivot(data: newVector[int](size), permutations:0) - for i in 0...size: - result[i] = i - -type AbtractMAtrix = Matrix or SMAtrix - -iterator iter_walter[T](m: Matrix[T]) : (int,int, T) {.closure.} = - for i in 0...m.rows: - for j in 0...m.columns: - yield (i, j, m[i,j]) - -iterator items*[T](m: Matrix[T]): (int,int, T) = - for i in 0...m.rows: - for j in 0...m.columns: - yield (i, j, m[i,j]) - -proc size*[T](m : Matrix[T]) : (int, int) = (m.rows, m.columns) - -proc newMatrix*[T](rows, columns : int, init : T = 0) : Matrix[T] = - result = Matrix[T](rows: rows, columns:columns, data: newSeq[T](rows * columns)) - for i,j,_ in items(result): - result[i,j] = init - -proc newMatrix*[T](rows, columns : int, values : openarray[T]) : Matrix[T] = - result = Matrix[T](rows: rows, columns:columns, data: newSeq[T](rows * columns)) - for i,j,_ in items(result): - result[i,j] = values[i * columns + j] - -proc identity*[T](sz : int) : Matrix[T] = - result = newMatrix[T](sz,sz,0) - for i in 0...sz: - result[i,i] = 1 - - -proc `[]`*[T](m : Matrix[T], r,c :int) : T = - m.data[r*m.columns + c] - -proc `[]`*[T](m : var Matrix[T], r,c :int) : var T = - m.data[r*m.columns + c] - -proc `[]=`*[T](m : var Matrix[T], r,c :int, newValue : T) = - m.data[r*m.columns + c] = newValue - -# proc `[]=`*[T](m : var Matrix[T], r,c :int, newValue : T) = -# m.data[r*m.columns + c] = newValue - -proc `$`*[T](m : Matrix[T]) : string = - result = "[" - for i,j,v in items(m): - if j == 0: - if i > 0: result &= " " - result &= "[" - result &= $v - if j == (m.columns - 1): - result &= "]" - if i != (m.rows - 1): - result &= ",\n" - else: - result &= "]\n" - else: result &= ", " - -proc `*`*[T](m1 : Matrix[T], m2 : Matrix[T]) : Matrix[T] = - result = newMatrix[T](m1.rows, m2.columns, 0) - for i in 0...result.rows: - for j in 0...result.columns: - for k in 0...m1.columns: - result[i, j] = result[i, j] + m1[i, k] * m2[k, j] - -proc `*`*[T](m1 : Matrix[T], v2 : Vector[T]) : Vector[T] = - result = newVector[T](m1.rows, 0) - for i in 0...m1.rows: - for j in 0...m1.columns: - result[i] = result[i] + m1[i, j] * v2[j] - -proc `+`*[T](m1 : Matrix[T], m2 : Matrix[T]) : Matrix[T] = - result = newMatrix[T](m1.rows, m1.columns) - for i in 0...m1.rows: - for j in 0...m1.columns: - result[i,j] = m1[i,j] + m2[i,j] - -proc `+`*[T](m1 : Matrix[T], v : T) : Matrix[T] = - result = newMatrix[T](m1.rows, m1.columns) - for i in 0...m1.rows: - for j in 0...m1.columns: - result[i,j] = m1[i,j] + v - - -proc `-`*[T](m1 : Matrix[T], m2 : Matrix[T]) : Matrix[T] = - result = newMatrix[T](m1.rows, m2.columns) - for i in 0...m1.rows: - for j in 0...m1.columns: - result[i,j] = m1[i,j] - m2[i,j] - -proc `-`*[T](m1 : Matrix[T], v : T) : Matrix[T] = - result = newMatrix[T](m1.rows, m1.columns) - for i in 0...m1.rows: - for j in 0...m1.columns: - result[i,j] = m1[i,j] - v - -proc `+=`*[T](m1 : var Matrix[T], m2 : Matrix[T]) = - for i in 0...m1.rows: - for j in 0...m1.columns: - m1[i,j] += m2[i,j] - -proc `+=`*[T](m1 : var Matrix[T], v : T) : Matrix[T] = - for i in 0...m1.rows: - for j in 0...m1.columns: - m1[i,j] += v - -proc `-=`*[T](m1 : var Matrix[T], v : T) : Matrix[T] = - for i in 0...m1.rows: - for j in 0...m1.columns: - m1[i,j] -= v - -proc `-=`*[T](m1 : var Matrix[T], m2 : Matrix[T]) = - for i in 0...m1.rows: - for j in 0...m1.columns: - m1[i,j] -= m2[i,j] - -proc `-`*[T](m : Matrix[T]) : Matrix[T] = - result = newMatrix[T](m.rows, m.columns) - for i,j,v in m: - result[i,j] = -v - -proc `==`*[T](m1 : Matrix[T], m2 : Matrix[T]) : bool = - if m1.size() != m2.size(): - return false - for i in 0...m1.rows: - for j in 0...m1.columns: - if m1[i,j] != m2[i,j]: - return false - return true - -proc clone*[T](m : Matrix[T]) : Matrix[T] = newMatrix[T](m.rows, m.columns, m.data) - -proc transpose*[T](m : Matrix[T]) : Matrix[T] = - result = newMatrix[T](m.rows, m.columns) - for i, j, v in m: - result[j, i] = v - -proc det*[T](m : Matrix[T]) : T = - var clone = m.clone() - clone.gauss_jordan_low() - result = 1 - for i in 0...clone.rows: - result *= clone[i, i] - -proc swap_rows[T](m : var Matrix[T], id1 : int, id2 : int, pivot : Pivot=nil, other : MatrixRef[T]=nil) = - for i in 0...m.columns: - let tmp = m[id1, i] - m[id1, i] = m[id2, i] - m[id2, i] = tmp - if other != nil: - other[].swap_rows(id1, id2) - if pivot != nil: - var pv = pivot - let tmp = pv[id1] - pv[id1] = pv[id2] - pv[id2] = tmp - pv.permutations = pv.permutations + 1 - -proc add_row[T](m : var Matrix[T], sourceIndex : int, destIndex : int, factor : T, other : MatrixRef[T]=nil) = - for i in 0...m.columns: - m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor - if other != nil: - other[].add_row(sourceIndex, destIndex, factor) - -proc gauss_jordan_low*[T](m : var Matrix[T], other : MatrixRef[T]=nil) = - var pivot = newPivot[T](m.rows) - for i in 0...m.rows: - if m[i, i] == 0: - for j in (i + 1)...m.columns: - if m[j, i] != 0: - m.swap_rows(i, j, pivot, other) - break - for j in (i + 1)...m.rows: - if m[i, i] != 0: - let factor = -m[j, i] / m[i, i] - m.add_row(i, j, factor, other) - -proc gauss_jordan_high*[T](m : var Matrix[T], other : MatrixRef[T]=nil) = - var pivot = newPivot[T](m.rows) - for i in m.rows-->0: - if m[i, i] == 0: - for j in i-->0: - if m[j, i] != 0: - m.swap_rows(i, j, pivot, other) - break - for j in i-->0: - if m[i, i] != 0: - let factor = -m[j, i] / m[i, i] - m.add_row(i, j, factor, other) - -proc invert*[T](m : Matrix[T]) : Matrix[T] = - var tmp = m.clone() - var res = new(MatrixRef[T]) - res[] = identity[T](tmp.rows) - tmp.gauss_jordan_low(res) - tmp.gauss_jordan_high(res) - for i in 0...res.rows: - let f = tmp[i, i] - for j in 0...res.columns: - res[][i, j] /= f - res[] - -proc triu*[T](m : Matrix[T], diag_replace=nil) : Matrix[T] = - result = Matrix(m.rows, m.columns, 0) - for i in range(m.rows): - for j in range(i, m.columns): - if diag_replace and i == j: - result[i, j] = diag_replace - else: - result[i, j] = m[i, j] - -proc tril*[T](m : Matrix[T], diag_replace=nil) : Matrix[T] = - result = Matrix(m.rows, m.columns, 0) - for i in range(m.rows): - for j in range(i + 1): - if diag_replace and i == j: - result[i, j] = diag_replace - else: - result[i, j] = m[i, j] - -proc lu_row[T](m : var Matrix[T], i : int) = - if m[i, i] == 0: - raise newException(SingularMatrixError, "Matrix is singular") - for j in i...m.columns: - for k in 0...i: - m[i, j] = m[i, j] - m[i, k] * m[k, j] - for j in (i + 1)...m.columns: - for k in 0...i: - m[j, i] = m[j, i] - m[j, k] * m[k, i] - m[j, i] = m[j, i] / m[i, i] - -proc lu_pivot[T](m : var Matrix[T], i : int, pivot : Pivot) = - var max = abs(m[i, i]) - var max_index = i - for j in (i + 1)...m.rows: - if abs(m[j, i]) > max: - max = abs(m[i, j]) - max_index = j - if max_index != i: - m.swap_rows(i, max_index, pivot) - -proc lu*[T](m : var Matrix[T], pivoting=true) : Pivot = - var pivot = newPivot[T](m.rows) - if pivoting: - for i in 0...m.rows: - m.lu_pivot(i, pivot) - m.lu_row(i) - else: - for i in 0...m.rows: - m.lu_row(i) - return pivot - -proc lu_solve*[T](m : Matrix[T], b : Vector[T], p : Pivot = nil) : Vector[T] = - var pivot = p - if pivot == nil: - pivot = new(Pivot) - pivot = newPivot[T](m.rows) - var x = newVector[T](m.rows) - for i in 0...m.rows: - x[i] = b[pivot[i]] - for k in 0...i: - x[i] = x[i] - m[i, k] * x[k] - for i in m.rows-->0: - for k in (i + 1)...m.rows: - x[i] = x[i] - m[i, k] * x[k] - x[i] = x[i] / m[i, i] - return x - -proc lu_invert*[T](m : Matrix[T], pivot : Pivot=nil) : Matrix[T] = - if not pivot: - pivot = newPivot[T](m.rows) - result = newMatrix[T](m.rows, m.columns) - for i in 0...m.rows: - for j in 0...m.rows: - if pivot[j] == i: - result[j, i] = 1 - else: - result[j, i] = 0 - for k in range(j): - result[j, i] -= m[j, k] * result[k, i] - for j in m.rows-->0: - for k in range(j + 1, m.rows): - result[j, i] -= m[j, k] * result[k, i] - result[j, i] = result[j, i] / m[j, j] - -proc lu_det*[T](m : Matrix[T], pivot : Pivot=nil) : T = - if not pivot: - pivot = newPivot[T](m.rows) - result = 1 - for i in 0...m.rows: - result *= m[i, i] - if pivot.permutations mod 2 != 0: - result *= -1 - -proc from_pivot[T](pivot : Pivot): Matrix[T] = - result = Matrix(len(pivot), len(pivot)) - for i in 0...pivot.len: - result[pivot[i], i] = 1 - # result[j, i] = 1 - diff --git a/matrix_test.nim b/matrix_test.nim deleted file mode 100644 index f4feff8..0000000 --- a/matrix_test.nim +++ /dev/null @@ -1,183 +0,0 @@ -from random import randomize, random -from utils import `...` -import matrix -from matrix import det -from vector import createVector, newVector, `-`, abs, norm -import unittest - -suite "Nim linear algebra library": - echo "suite setup: run once before the tests" - randomize() - var mtx : Matrix[float64] - - setup: - echo "run before each test" - let DIM = 100 - var numbers = newSeq[float64](DIM * DIM) - for i in 0...len(numbers): - numbers[i] = float64(random(-DIM..DIM)) - mtx = newMatrix[float64](DIM,DIM,numbers) - - teardown: - echo "run after each test" - - test "LU decomposition": - var lu = mtx.clone() - let pivot = lu.lu() - - test "Linear system solve": - # let nums = [2.0,1.0,3.0,2.0,6.0,8.0,6.0,8.0,18.0] - # let nums = [-3.0, -1.0, 0.0, -2.0, 0.0, 2.0, -3.0, -1.0, 0.0] - let DIM = 100 - var numbers = newSeq[float64](DIM * DIM) - for i in 0...len(numbers): - numbers[i] = float64(random(-DIM..DIM)) - - var mtx = newMatrix[float64](DIM,DIM,numbers) - var lu = mtx.clone() - let pivot = lu.lu() - var b = newVector[float64](DIM) - for i in 0...len(b): - b[i] = float64(random(DIM)) - - let x = lu.lu_solve(b, pivot) - let error = (mtx * x) - b - - check(error.norm() < 1e-5) - # var mtx2 = mtx.clone - # for i in 0...10000: - # var add = random(-DIM..DIM).float64() - # discard mtx2 += add - - # echo mtx - # mtx[1,1] += 10000.0 - # echo mtx - # echo mtx.det() - - # echo lu - - # give up and stop if this fails - - echo "suite teardown: run once after the tests" - -# let m1 = newMatrix[float](3,3,[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0]) -# var m2 = m1.clone() -# let m3 = m2.clone() -# m2[2,1] = -25 -# echo m1 + m2 -# echo m2 -# echo m1.det() - -# let nums = [-63, 3, 70, -23, 55, -# -100, -37, 81, -98, 84, -# -36, -45, -70, 98, -18, -# -15, 92, 82, 85, -2, -# 45, 54, -22, 27, 0 -# ] - -# let s = nums.map(proc(n : int) : float = float(n)) -# let s2 = nums.map(n => float(n)) - -# let m4 = newMatrix[float](5,5,s) -# echo m4 -# echo m4.det - -# var m5 = m4.clone() -# var pivot = m5.lu() -# let b = createVector[float](1.0,2.0,3.0,4.0,5.0) -# let x = m5.lu_solve(b, pivot) -# echo x -# proc new(T: typedesc): ref T = -# echo "allocating " -# new(result) - -# var n = new Vector[int] - -# type -# Index = distinct int - -# proc `==` (a, b: Index): bool {.borrow.} - -# var af = (0, 0.Index) -# var b = (0, 0.Index) - -# echo af == b # works! - - -# type Person = ref object of RootObj -# name : string -# age : int - -# type Employee = ref object of Person -# salary: int - - -# type RecordType = tuple or object - -# proc printFields(rec: RecordType) = -# for key, value in fieldPairs(rec): -# echo key, " = ", value - -# proc printFields(rec: ref[RecordType]) = -# for key, value in fieldPairs(rec[]): -# echo key, " = ", value - - -# let p = Person(name : "Walter", age : 28) -# let e = Employee(name : "Walter", age : 28, salary: 45000) - -# let people : seq[Person] = @[p, e] - -# for person in people: -# printFields person - -# let DIM = 5 -# var numbers = newSeq[float](DIM * DIM) -# for i in 0...len(numbers): -# numbers[i] = float(random(DIM)) - -# let m = newMatrix[float](DIM,DIM, numbers) - -# var m2 = new(Matrix[float]) -# m2[] = m.clone() -# echo m2[] -# m2[][0,0] = 300 -# echo m2[] - -# let inverse = m2[].invert() - -# echo m -# echo m.det -# echo m * inverse -# echo inverse * m - -# let nums = [2.0,1.0,3.0,2.0,6.0,8.0,6.0,8.0,18.0] -# let nums = [-3.0, -1.0, 0.0, -2.0, 0.0, 2.0, -3.0, -1.0, 0.0] -# let DIM = 600 -# var numbers = newSeq[float64](DIM * DIM) -# for i in 0...len(numbers): -# numbers[i] = float64(random(-DIM..DIM)) - -# var mtx = newMatrix[float64](DIM,DIM,numbers) -# var mtx2 = mtx.clone -# for i in 0...10000: -# var add = random(-DIM..DIM).float64() -# discard mtx2 += add - -# echo mtx -# mtx[1,1] += 10000.0 -# echo mtx -# echo mtx.det() - -# var lu = mtx.clone() -# let pivot = lu.lu() -# # echo lu - -# var b = newVector[float64](DIM) -# for i in 0...len(b): -# b[i] = float64(random(DIM)) - -# let x = lu.lu_solve(b, pivot) -# let error = (mtx * x) - b - -# echo abs(error) \ No newline at end of file diff --git a/mmath.nimble b/mmath.nimble new file mode 100644 index 0000000..c4f5e2b --- /dev/null +++ b/mmath.nimble @@ -0,0 +1,13 @@ +# Package + +version = "0.1.0" +author = "Walter Oggioni" +description = "Small linear algebra library" +license = "MIT" +srcDir = "src" + + +# Dependencies + +requires "nim >= 0.18" +requires "nwo >= 0.1" diff --git a/src/mmath/hmatrix.nim b/src/mmath/hmatrix.nim new file mode 100644 index 0000000..8f60cd2 --- /dev/null +++ b/src/mmath/hmatrix.nim @@ -0,0 +1,398 @@ +from sequtils import newSeqWith, map +from nwo/utils import `...`, `-->`, box +from hvector import HVector, newHVector, buildHVector +from pivot import HPivot, newHPivot, `[]`, `[]=`, len, SizeError, SingularMatrixError +from options import Option, none, some +from math import sqrt + +type + HMatrix*[T] = object + rows, columns : int + data : seq[T] + +proc size*[T](m : HMatrix[T]) : (int, int) = (m.rows, m.columns) + +proc `[]`*[T](m : HMatrix[T], r, c :int) : T = + m.data[r * m.columns + c] + +proc `[]`*[T](m : var HMatrix[T], r, c :int) : var T = + m.data[r * m.columns + c] + +proc `[]=`*[T](m : var HMatrix[T], r, c :int, newValue : T) = + m.data[r * m.columns + c] = newValue + +iterator items*[T](m: HMatrix[T]): (int, int, T) = + for i in 0...m.rows: + for j in 0...m.columns: + yield (i, j, m[i, j]) + +proc newHMatrix*[T](rows, columns : int) : HMatrix[T] = + HMatrix[T](rows: rows, columns:columns, data: newSeq[T](rows * columns)) + +proc newHMatrix*[T](rows, columns : int, init : T) : HMatrix[T] = + result = newHMatrix[T](rows, columns) + for i,j,_ in items(result): + result[i,j] = init + +proc newHMatrix*[T](rows, columns : int, values : openarray[T]) : HMatrix[T] = + result = newHMatrix[T](rows, columns) + for i, j, _ in items(result): + result[i,j] = values[i * columns + j] + +proc identity*[T](sz : int) : HMatrix[T] = + result = newHMatrix[T](sz,sz,0) + for i in 0...sz: + result[i,i] = 1 + +proc `$`*[T](m : HMatrix[T]) : string = + result = "[" + for i,j,v in items(m): + if j == 0: + if i > 0: result &= " " + result &= "[" + result &= $v + if j == (m.columns - 1): + result &= "]" + if i != (m.rows - 1): + result &= ",\n" + else: + result &= "]\n" + else: result &= ", " + +proc `*`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m1.rows, m2.columns, 0) + for i in 0...result.rows: + for j in 0...result.columns: + for k in 0...m1.columns: + result[i, j] = result[i, j] + m1[i, k] * m2[k, j] + +proc `*`*[T](m1 : HMatrix[T], v2 : HVector[T]) : HVector[T] = + result = newHVector[T](m1.rows, 0) + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i] = result[i] + m1[i, j] * v2[j] + +proc `+`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m1.rows, m1.columns) + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] + m2[i,j] + +proc `+`*[T](m1 : HMatrix[T], v : T) : HMatrix[T] = + result = newHMatrix[T](m1.rows, m1.columns) + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] + v + + +proc `-`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m1.rows, m2.columns) + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] - m2[i,j] + +proc `-`*[T](m1 : HMatrix[T], v : T) : HMatrix[T] = + result = newHMatrix[T](m1.rows, m1.columns) + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] - v + +proc `+=`*[T](m1 : var HMatrix[T], m2 : HMatrix[T]) = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] += m2[i,j] + +proc `+=`*[T](m1 : var HMatrix[T], v : T) : HMatrix[T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] += v + +proc `-=`*[T](m1 : var HMatrix[T], v : T) : HMatrix[T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] -= v + +proc `-=`*[T](m1 : var HMatrix[T], m2 : HMatrix[T]) = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] -= m2[i,j] + +proc `-`*[T](m : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m.rows, m.columns) + for i,j,v in m: + result[i,j] = -v + +proc `==`*[T](m1 : HMatrix[T], m2 : HMatrix[T]) : bool = + if m1.size() != m2.size(): + return false + for i in 0...m1.rows: + for j in 0...m1.columns: + if m1[i,j] != m2[i,j]: + return false + return true + +proc clone*[T](m : HMatrix[T]) : HMatrix[T] = newHMatrix[T](m.rows, m.columns, m.data) + +proc transpose*[T](m : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m.columns, m.rows) + for i, j, v in items(m): + result[j, i] = v + +proc swap_rows[T](m : var HMatrix[T], id1 : int, id2 : int) = + for i in 0...m.columns: + let tmp = m[id1, i] + m[id1, i] = m[id2, i] + m[id2, i] = tmp + +proc swap_rows[T](m : var HMatrix[T], id1 : int, id2 : int, pivot : var HPivot[T]) = + m.swap_rows(id1, id2) + let tmp = pivot[id1] + pivot[id1] = pivot[id2] + pivot[id2] = tmp + pivot.permutations += 1 + +proc swap_rows[T]( + m : var HMatrix[T], id1 : int, + id2 : int, + pivot : var HPivot[T], + other : var HMatrix[T]) = + m.swap_rows(id1, id2, pivot) + other.swap_rows(id1, id2) + +proc add_row[T]( + m : var HMatrix[T], + sourceIndex : int, + destIndex : int, + factor : T) = + for i in 0...m.columns: + m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor + +proc add_row[T]( + m : var HMatrix[T], + sourceIndex : int, + destIndex : int, + factor : T, + other : var HMatrix[T]) = + add_row(m, source_index, dest_index, factor) + other.add_row(sourceIndex, destIndex, factor) + +proc gauss_jordan_low*[T]( + m : var HMatrix[T], + other : var HMatrix[T]) = + var pivot = newHPivot[T](m.rows) + for i in 0...m.rows: + if m[i, i] == 0: + for j in (i + 1)...m.columns: + if m[j, i] != 0: + m.swap_rows(i, j, pivot, other) + break + for j in (i + 1)...m.rows: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor, other) + +proc gauss_jordan_low*[T]( + m : var HMatrix[T]) = + var pivot = newHPivot[T](m.rows) + for i in 0...m.rows: + if m[i, i] == 0: + for j in (i + 1)...m.columns: + if m[j, i] != 0: + m.swap_rows(i, j, pivot) + break + for j in (i + 1)...m.rows: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor) + +proc gauss_jordan_high*[T]( + m : var HMatrix[T]) = + var pivot = newHPivot[T]() + for i in m.rows-->0: + if m[i, i] == 0: + for j in i-->0: + if m[j, i] != 0: + m.swap_rows(i, j, pivot) + break + for j in i-->0: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor) + +proc gauss_jordan_high*[T]( + m : var HMatrix[T], + other : var HMatrix[T]) = + var pivot = newHPivot[T](m.rows) + for i in m.rows-->0: + if m[i, i] == 0: + for j in i-->0: + if m[j, i] != 0: + m.swap_rows(i, j, pivot, other) + break + for j in i-->0: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor, other) + +proc det*[T](m : HMatrix[T]) : T = + if m.rows != m.columns: + raise newException(SizeError, "Matrix must be square in order to compute the determinant") + var clone = m.clone() + clone.gauss_jordan_low() + result = 1 + for i in 0...clone.rows: + result *= clone[i, i] + +proc invert*[T](m : HMatrix[T]) : HMatrix[T] = + if m.rows != m.columns: + raise newException(SizeError, "Matrix must be square in order to compute the determinant") + var tmp = m.clone() + result = identity[T](m.rows) + tmp.gauss_jordan_low(result) + tmp.gauss_jordan_high(result) + for i in 0...result.rows: + let f = tmp[i, i] + for j in 0...result.columns: + result[i, j] /= f + +proc triu*[T](m : HMatrix[T], diag_replace: T) : HMatrix[T] = + result = newHMatrix[T](m.rows, m.columns, 0) + for i in 0...m.rows: + for j in i...m.columns: + if i == j: + result[i, j] = diag_replace + else: + result[i, j] = m[i, j] + +proc triu*[T](m : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m.rows, m.columns, 0) + for i in 0...m.rows: + for j in i...m.columns: + result[i, j] = m[i, j] + + +proc tril*[T](m : HMatrix[T], diag_replacement : T) : HMatrix[T] = + result = newHMatrix[T](m.rows, m.columns, 0) + for i in 0...m.rows: + for j in 0...(i + 1): + if i == j: + result[i, j] = diag_replacement + else: + result[i, j] = m[i, j] + +proc tril*[T](m : HMatrix[T]) : HMatrix[T] = + result = newHMatrix[T](m.rows, m.columns, 0) + for i in 0...m.rows: + for j in 0...(i + 1): + result[i, j] = m[i, j] + +proc lu_row[T](m : var HMatrix[T], i : int) = + if m[i, i] == 0: + raise newException(SingularMatrixError, "Matrix is singular") + for j in i...m.columns: + for k in 0...i: + m[i, j] = m[i, j] - m[i, k] * m[k, j] + for j in (i + 1)...m.columns: + for k in 0...i: + m[j, i] = m[j, i] - m[j, k] * m[k, i] + m[j, i] /= m[i, i] + +proc lu_pivot[T](m : var HMatrix[T], i : int, pivot : var HPivot[T]) = + var max = abs(m[i, i]) + var max_index = i + for j in (i + 1)...m.rows: + if abs(m[j, i]) > max: + max = abs(m[i, j]) + max_index = j + if max_index != i: + m.swap_rows(i, max_index, pivot) + +proc lup*[T](m : var HMatrix[T]) : HPivot[T] = + result = newHPivot[T](m.rows) + for i in 0...m.rows: + m.lu_pivot(i, result) + m.lu_row(i) + +proc lu*[T](m : var HMatrix[T]) = + for i in 0...m.rows: + m.lu_row(i) + +proc lu_solve*[T](m : HMatrix[T], b : HVector[T], pivot : HPivot[T]) : HVector[T] = + var x = newHVector[T](m.rows) + for i in 0...m.rows: + x[i] = b[pivot[i]] + for k in 0...i: + x[i] = x[i] - m[i, k] * x[k] + for i in m.rows --> 0: + for k in (i + 1)...m.rows: + x[i] = x[i] - m[i, k] * x[k] + if m[i,i] != 0: + x[i] /= m[i, i] + else: + raise newException(SingularMatrixError, "Matrix is singular") + return x + +proc lu_solve*[T]( + m : HMatrix[T], + b : HVector[T]) : HVector[T] = + var pivot = newHPivot[T]() + lu_solve(m, b, pivot) + +proc lu_invert*[T](m : HMatrix[T], pivot : HPivot[T]) : HMatrix[T] = + if m.rows != m.columns: + raise newException(SizeError, "Matrix must be square in order to compute the inverse") + result = newHMatrix[T](m.rows, m.columns) + for i in 0...m.rows: + for j in 0...m.rows: + if pivot[j] == i: + result[j, i] = 1 + else: + result[j, i] = 0 + for k in range(j): + result[j, i] -= m[j, k] * result[k, i] + for j in m.rows-->0: + for k in range(j + 1, m.rows): + result[j, i] -= m[j, k] * result[k, i] + result[j, i] = result[j, i] / m[j, j] + +proc lu_invert*[T](m : HMatrix[T]) : HMatrix[T] = lu_invert(m, newHPivot[T]()) + +proc lu_det*[T](m : var HMatrix[T]) : T = + if m.rows != m.columns: + raise newException(SizeError, "Matrix must be square in order to compute the determinant") + let pivot = m.lup() + result = 1 + for i in 0...m.rows: + result *= m[i, i] + if pivot.permutations mod 2 != 0: + result *= -1 + +proc lu_det*[T](m : HMatrix[T]) : T = + if m.rows != m.columns: + raise newException(SizeError, "Matrix must be square in order to compute the determinant") + var clone = m.clone() + lu_det(clone) + +proc squared_norm2*[T](m : HMatrix[T]): T = + result = T(0) + for i, j, v in items(m): + result += v * v + +proc norm2*[T](m : HMatrix[T]): T = + sqrt(m.squared_norm2()) + +proc `*`*[T](pivot : HPivot[T], m : HMatrix[T]) : HMatrix[T] = + result = m.clone() + var pclone = pivot + for i in 0...pclone.len(): + while i != pclone[i]: + result.swap_rows(i, pclone[i]) + let tmp = pclone[i] + pclone[i] = pclone[tmp] + pclone[tmp] = tmp + +proc from_pivot*[T](pivot : HPivot[T]): HMatrix[T] = + result = newHMatrix[T](len(pivot), len(pivot)) + for i in 0...pivot.len: + result[pivot[i], i] = 1 + diff --git a/src/mmath/hvector.nim b/src/mmath/hvector.nim new file mode 100644 index 0000000..a37adf6 --- /dev/null +++ b/src/mmath/hvector.nim @@ -0,0 +1,48 @@ +from nwo/utils import `...` +from sequtils import newSeqWith +from math import sqrt + +type HVector*[T] = seq[T] + +proc newHVector*[T](size : int, init: T=0) : HVector[T] = + result = newSeq[T](size) + for i in 0...len(result): + result[i] = init + +proc `+`*[T](v1 : HVector[T], v2:HVector[T]) : HVector[T] = + result = newHVector[T](len(v1)) + for i in 0...len(v1): + result[i] = v1[i] + v2[i] + +proc `-`*[T](v1 : HVector[T], v2:HVector[T]) : HVector[T] = + result = newHVector[T](len(v1)) + for i in 0...len(v1): + result[i] = v1[i] - v2[i] + +proc `*`*[T](v1 : HVector[T], v2:HVector[T]) : T = + result = 0 + for i in 0...len(v1): + result += v1[i] * v2[i] + +proc `+=`*[T](v1 : var HVector[T], v2 : HVector[T]) = + for i in 0...len(v1): + v1[i] += v2[i] + +proc `-=`*[T](v1 : var HVector[T], v2 : HVector[T]) = + for i in 0...len(v1): + v1[i] -= v2[i] + +proc `+=`*[T](v : var HVector[T], value : T) = v.add(value) + +proc buildHVector*[T](elems : varargs[T]) : HVector[T] = + result = newSeq[T]() + for elem in items(elems): + result += elem + +proc norm*[T](v : HVector[T]) : T = + for value in v: + result += v * v + +proc abs*[T](v : HVector[T]) : T = + return math.sqrt(v.norm) + diff --git a/src/mmath/pivot.nim b/src/mmath/pivot.nim new file mode 100644 index 0000000..bcd142f --- /dev/null +++ b/src/mmath/pivot.nim @@ -0,0 +1,30 @@ +from hvector import HVector, newHVector +from nwo/utils import `...` + +type + HPivot*[T] = object + data : HVector[int] + permutations* : int + + SPivot*[SIZE : static[int], T] = object + data : array[SIZE, int] + permutations* : int + SingularMatrixError* = object of ValueError + SizeError* = object of ValueError + +proc `[]`*[T](p : HPivot[T], index : int) : int = p.data[index] +proc `[]=`*[T](p : var HPivot[T], index : int, value : int) = p.data[index] = value +proc len*[T](p : HPivot[T]) : int = p.data.len +proc `$`*[T](pivot : HPivot[T]) : string = $pivot.data +proc newHPivot*[T](size : int) : HPivot[T] = + result = HPivot[T](data: newHVector[int](size), permutations:0) + for i in 0...size: + result[i] = i + +proc `[]`*[SIZE, T](p : SPivot[SIZE, T], index : int) : int = p.data[index] +proc `[]=`*[SIZE, T](p : var SPivot[SIZE, T], index : int, value : int) = p.data[index] = value +proc len*[SIZE, T](p : SPivot[SIZE, T]) : int = SIZE +proc `$`*[SIZE, T](pivot : SPivot[SIZE, T]) : string = $pivot.data +proc newSPivot*[SIZE, T]() : SPivot[SIZE, T] = + for i in 0...SIZE: + result[i] = i \ No newline at end of file diff --git a/src/mmath/smatrix.nim b/src/mmath/smatrix.nim new file mode 100644 index 0000000..1e90281 --- /dev/null +++ b/src/mmath/smatrix.nim @@ -0,0 +1,370 @@ +from sequtils import newSeqWith, map +from nwo/utils import `...`, `-->`, box +from svector import SVector +from pivot import SPivot, newSPivot, `[]`, `[]=`, SingularMatrixError, len +from math import sqrt +import random + +type + SMatrix*[ROWS, COLUMNS: static[int], T] = object + data : array[0..(ROWS*COLUMNS - 1), T] + SquareSMatrix[SIZE: static[int], T] = SMatrix[SIZE, SIZE, T] + +proc size*[ROWS, COLUMNS : static[int], T](m : SMatrix) : (int, int) = (m.rows, m.columns) + +proc rows*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : int = ROWS +proc columns*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : int = COLUMNS + +proc `[]`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS,COLUMNS,T], r, c :int) : T = + m.data[r * COLUMNS + c] + +proc `[]`*[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS,COLUMNS,T], r, c :int) : var T = + m.data[r * COLUMNS + c] + +proc `[]=`*[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS,COLUMNS,T], r, c :int, newValue : T) = + m.data[r * COLUMNS + c] = newValue + +iterator items*[ROWS, COLUMNS : static[int], T](m: SMatrix[ROWS,COLUMNS,T]): (int, int, T) = + for i in 0...m.rows: + for j in 0...m.columns: + yield (i, j, m[i, j]) + +proc newSMatrix*[ROWS, COLUMNS : static[int], T](init : T) : SMatrix[ROWS, COLUMNS, T] = + for i,j,_ in items(result): + result[i,j] = init + +proc newSMatrixFromArray*[ROWS, COLUMNS : static[int], T](values : array[0..(ROWS * COLUMNS - 1), T]) : auto = + SMatrix[ROWS, COLUMNS, T](data:values) + +proc identity*[SIZE: static[int], T]() : SquareSMatrix[SIZE, T] = + for i in 0...SIZE: + result[i,i] = 1 + +proc `$`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : string = + result = "[" + for i,j,v in items(m): + if j == 0: + if i > 0: result &= " " + result &= "[" + result &= $v + if j == (m.columns - 1): + result &= "]" + if i != (m.rows - 1): + result &= ",\n" + else: + result &= "]\n" + else: result &= ", " + +proc `*`*[ROWS1, COLUMNS2, COMMON : static[int], T]( + m1 : SMatrix[ROWS1, COMMON, T], + m2 : SMatrix[COMMON, COLUMNS2, T]) : SMatrix[ROWS1, COLUMNS2, T] = + for i in 0...result.rows: + for j in 0...result.columns: + for k in 0...m1.columns: + result[i, j] = result[i, j] + m1[i, k] * m2[k, j] + +proc `*`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], v2 : SVector[ROWS, T]) : SVector[ROWS, T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i] += m1[i, j] * v2[j] + +proc `+`*[ROWS, COLUMNS : static[int], T](addend1 : SMatrix[ROWS, COLUMNS, T], addend2 : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] = + for i in 0...addend1.rows: + for j in 0...addend1.columns: + result[i,j] = addend1[i,j] + addend2[i,j] + +proc `+`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] + v + +proc `-`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS,T], m2 : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS,T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] - m2[i,j] + +proc `-`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + result[i,j] = m1[i,j] - v + +proc `+=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], m2 : SMatrix[ROWS, COLUMNS, T]) = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] += m2[i,j] + +proc `+=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] += v + +proc `-=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], v : T) : SMatrix[ROWS, COLUMNS, T] = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] -= v + +proc `-=`*[ROWS, COLUMNS : static[int], T](m1 : var SMatrix[ROWS, COLUMNS, T], m2 : SMatrix[ROWS, COLUMNS, T]) = + for i in 0...m1.rows: + for j in 0...m1.columns: + m1[i,j] -= m2[i,j] + +proc `-`*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] = + for i,j,v in m: + result[i,j] = -v + +proc `==`*[ROWS, COLUMNS : static[int], T](m1 : SMatrix[ROWS, COLUMNS, T], m2 : SMatrix[ROWS, COLUMNS, T]) : bool = + for i in 0...m1.rows: + for j in 0...m1.columns: + if m1[i,j] != m2[i,j]: + return false + return true + +proc squared_norm2*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]): T = + result = T(0) + for i, j, v in items(m): + result += v * v + +proc norm2*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]): T = + sqrt(m.squared_norm2()) + +proc clone*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] = + newSMatrixFromArray[ROWS, COLUMNS, T](m.data) + +proc transpose*[ROWS, COLUMNS : static[int], T](m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[COLUMNS, ROWS, T] = + for i, j, v in items(m): + result[j, i] = v + +proc swap_rows[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS, COLUMNS, T], id1 : int, id2 : int) = + for i in 0...m.columns: + let tmp = m[id1, i] + m[id1, i] = m[id2, i] + m[id2, i] = tmp + +proc swap_rows[ROWS, COLUMNS : static[int], T](m : var SMatrix[ROWS, COLUMNS, T], id1 : int, id2 : int, pivot : var SPivot[ROWS, T]) = + m.swap_rows(id1, id2) + let tmp = pivot[id1] + pivot[id1] = pivot[id2] + pivot[id2] = tmp + pivot.permutations += 1 + +proc swap_rows[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T], id1 : int, + id2 : int, + pivot : var SPivot[ROWS, T], + other : var SMatrix[ROWS, COLUMNS,T]) = + m.swap_rows(id1, id2, pivot) + other.swap_rows(id1, id2) + +proc add_row[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T], + sourceIndex : int, + destIndex : int, + factor : T) = + for i in 0...m.columns: + m[destIndex, i] = m[destIndex, i] + m[sourceIndex, i] * factor + +proc add_row[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T], + sourceIndex : int, + destIndex : int, + factor : T, + other : var SMatrix[ROWS, COLUMNS, T]) = + add_row(m, source_index, dest_index, factor) + other.add_row(sourceIndex, destIndex, factor) + +proc gauss_jordan_low*[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T], + other : var SMatrix[ROWS, COLUMNS, T]) = + var pivot = newSPivot[ROWS, T]() + for i in 0...m.rows: + if m[i, i] == 0: + for j in (i + 1)...m.columns: + if m[j, i] != 0: + m.swap_rows(i, j, pivot, other) + break + for j in (i + 1)...m.rows: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor, other) + +proc gauss_jordan_low*[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T]) = + var pivot = newSPivot[ROWS, T]() + for i in 0...m.rows: + if m[i, i] == 0: + for j in (i + 1)...m.columns: + if m[j, i] != 0: + m.swap_rows(i, j, pivot) + break + for j in (i + 1)...m.rows: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor) + +proc gauss_jordan_high*[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T]) = + var pivot = newSPivot[ROWS, T]() + for i in m.rows-->0: + if m[i, i] == 0: + for j in i-->0: + if m[j, i] != 0: + m.swap_rows(i, j, pivot) + break + for j in i-->0: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor) + +proc gauss_jordan_high*[ROWS, COLUMNS : static[int], T]( + m : var SMatrix[ROWS, COLUMNS, T], + other : var SMatrix[ROWS, COLUMNS, T]) = + var pivot = newSPivot[ROWS, T]() + for i in m.rows-->0: + if m[i, i] == 0: + for j in i-->0: + if m[j, i] != 0: + m.swap_rows(i, j, pivot, other) + break + for j in i-->0: + if m[i, i] != 0: + let factor = -m[j, i] / m[i, i] + m.add_row(i, j, factor, other) + +proc det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T = + var clone = m.clone() + clone.gauss_jordan_low() + result = 1 + for i in 0...clone.rows: + result *= clone[i, i] + +proc invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] = + var tmp = m.clone() + result = identity[SIZE, T]() + tmp.gauss_jordan_low(result) + tmp.gauss_jordan_high(result) + for i in 0...result.rows: + let f = tmp[i, i] + for j in 0...result.columns: + result[i, j] /= f + +proc triu*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replace: T) : SquareSMatrix[SIZE, T] = + for i in 0...m.rows: + for j in i...m.columns: + if i == j: + result[i, j] = diag_replace + else: + result[i, j] = m[i, j] + +proc triu*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] = + for i in 0...m.rows: + for j in i...m.columns: + result[i, j] = m[i, j] + +proc tril*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], diag_replacement : T) : SquareSMatrix[SIZE, T] = + for i in 0...m.rows: + for j in 0...(i + 1): + if i == j: + result[i, j] = diag_replacement + else: + result[i, j] = m[i, j] + +proc tril*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] = + for i in 0...m.rows: + for j in 0...(i + 1): + result[i, j] = m[i, j] + +proc lu_row[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T], i : int) = + if m[i, i] == 0: + raise newException(SingularMatrixError, "Matrix is singular") + for j in i...m.columns: + for k in 0...i: + m[i, j] = m[i, j] - m[i, k] * m[k, j] + for j in (i + 1)...m.columns: + for k in 0...i: + m[j, i] = m[j, i] - m[j, k] * m[k, i] + m[j, i] /= m[i, i] + +proc lu_pivot[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T], i : int, pivot : var SPivot[SIZE, T]) = + var max = abs(m[i, i]) + var max_index = i + for j in (i + 1)...m.rows: + if abs(m[j, i]) > max: + max = abs(m[i, j]) + max_index = j + if max_index != i: + m.swap_rows(i, max_index, pivot) + +proc lu*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) = + for i in 0...m.rows: + m.lu_row(i) + +proc lup*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) : SPivot[SIZE, T] = + result = newSPivot[SIZE,T]() + for i in 0...m.rows: + m.lu_pivot(i, result) + m.lu_row(i) + +proc lu_solve*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], b : SVector[SIZE, T], pivot : SPivot[SIZE, T]) : SVector[SIZE, T] = + var x : SVector[SIZE, T] + for i in 0...m.rows: + x[i] = b[pivot[i]] + for k in 0...i: + x[i] = x[i] - m[i, k] * x[k] + for i in m.rows --> 0: + for k in (i + 1)...m.rows: + x[i] = x[i] - m[i, k] * x[k] + if m[i,i] != 0: + x[i] /= m[i, i] + else: + raise newException(SingularMatrixError, "Matrix is singular") + return x + +proc lu_solve*[SIZE : static[int], T]( + m : SquareSMatrix[SIZE, T], + b : SVector[SIZE, T]) : SVector[SIZE, T] = + var pivot = newSPivot[SIZE, T]() + lu_solve(m, b, pivot) + +proc lu_invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T], pivot : SPivot[SIZE, T]) : SquareSMatrix[SIZE, T] = + for i in 0...m.rows: + for j in 0...m.rows: + if pivot[j] == i: + result[j, i] = 1 + else: + result[j, i] = 0 + for k in range(j): + result[j, i] -= m[j, k] * result[k, i] + for j in m.rows-->0: + for k in range(j + 1, m.rows): + result[j, i] -= m[j, k] * result[k, i] + result[j, i] = result[j, i] / m[j, j] + +proc lu_invert*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : SquareSMatrix[SIZE, T] = lu_invert(m, newSPivot[SIZE, T]()) + +proc lu_det*[SIZE : static[int], T](m : var SquareSMatrix[SIZE, T]) : T = + let pivot = m.lup() + result = 1 + for i in 0...m.rows: + result *= m[i, i] + if pivot.permutations mod 2 != 0: + result *= -1 + +proc lu_det*[SIZE : static[int], T](m : SquareSMatrix[SIZE, T]) : T = + var clone = m.clone() + lu_det(clone) + +proc `*`*[ROWS, COLUMNS : static[int], T](pivot : SPivot[ROWS, T], m : SMatrix[ROWS, COLUMNS, T]) : SMatrix[ROWS, COLUMNS, T] = + result = m.clone() + var pclone = pivot + for i in 0...pclone.len(): + while i != pclone[i]: + result.swap_rows(i, pclone[i]) + let tmp = pclone[i] + pclone[i] = pclone[tmp] + pclone[tmp] = tmp + +proc from_pivot*[SIZE : static[int], T](pivot : SPivot[SIZE, T]): SquareSMatrix[SIZE, T] = + result = newSMatrix[T](len(pivot), len(pivot)) + for i in 0...SIZE: + result[pivot[i], i] = 1 + diff --git a/src/mmath/svector.nim b/src/mmath/svector.nim new file mode 100644 index 0000000..c98abe7 --- /dev/null +++ b/src/mmath/svector.nim @@ -0,0 +1,44 @@ +from nwo/utils import `...` +from sequtils import newSeqWith +from math import sqrt + +type SVector*[S : static[int], T] = array[0..(S-1), T] + +proc newSVector*[SIZE, T](init: T=0) : SVector[SIZE, T] = + for i in 0...len(result): + result[i] = init + +proc buildSVector*[SIZE, T](elems : varargs[T]) : SVector[SIZE, T] = + for i in 0..