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
This commit is contained in:
Walter Oggioni
329
matrix.nim
329
matrix.nim
@@ -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
|
||||
|
||||
183
matrix_test.nim
183
matrix_test.nim
@@ -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)
|
||||
13
mmath.nimble
Normal file
13
mmath.nimble
Normal file
@@ -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"
|
||||
398
src/mmath/hmatrix.nim
Normal file
398
src/mmath/hmatrix.nim
Normal file
@@ -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
|
||||
|
||||
48
src/mmath/hvector.nim
Normal file
48
src/mmath/hvector.nim
Normal file
@@ -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)
|
||||
|
||||
30
src/mmath/pivot.nim
Normal file
30
src/mmath/pivot.nim
Normal file
@@ -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
|
||||
370
src/mmath/smatrix.nim
Normal file
370
src/mmath/smatrix.nim
Normal file
@@ -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
|
||||
|
||||
44
src/mmath/svector.nim
Normal file
44
src/mmath/svector.nim
Normal file
@@ -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..<elems.len:
|
||||
result[i] = elems[i]
|
||||
|
||||
proc `+`*[SIZE, T](v1 : SVector[SIZE, T], v2: SVector[SIZE, T]) : SVector[SIZE, T] =
|
||||
for i in 0...len(v1):
|
||||
result[i] = v1[i] + v2[i]
|
||||
|
||||
proc `-`*[SIZE, T](v1 : SVector[SIZE, T], v2: SVector[SIZE, T]) : SVector[SIZE, T] =
|
||||
for i in 0...len(v1):
|
||||
result[i] = v1[i] - v2[i]
|
||||
|
||||
proc `*`*[SIZE, T](v1 : SVector[SIZE, T], v2: SVector[SIZE, T]) : T =
|
||||
result = 0
|
||||
for i in 0...len(v1):
|
||||
result += v1[i] * v2[i]
|
||||
|
||||
proc `+=`*[SIZE, T](v1 : var SVector[SIZE, T], v2: SVector[SIZE, T]) =
|
||||
for i in 0...len(v1):
|
||||
v1[i] += v2[i]
|
||||
|
||||
proc `-=`*[SIZE, T](v1 : var SVector[SIZE, T], v2: SVector[SIZE, T]) =
|
||||
for i in 0...len(v1):
|
||||
v1[i] -= v2[i]
|
||||
|
||||
proc `+=`*[SIZE, T](v : var SVector[SIZE, T], value : T) = v.add(value)
|
||||
|
||||
proc norm*[SIZE, T](v : SVector[SIZE, T]) : T =
|
||||
for value in v:
|
||||
result += v * v
|
||||
|
||||
proc abs*[SIZE, T](v : SVector[SIZE, T]) : T =
|
||||
return math.sqrt(v.norm)
|
||||
|
||||
1
tests/config.nims
Normal file
1
tests/config.nims
Normal file
@@ -0,0 +1 @@
|
||||
switch("path", "$projectDir/../src")
|
||||
114
tests/test_hmatrix.nim
Normal file
114
tests/test_hmatrix.nim
Normal file
@@ -0,0 +1,114 @@
|
||||
from random import initRand, rand
|
||||
from nwo/utils import `...`
|
||||
from mmath/hmatrix import det, lu, from_pivot, newHMatrix, HMatrix, lu_det, invert, identity,
|
||||
clone, tril, triu, lu_solve, `*`, `-`, `+`, `+=`, `-=`, `==`, norm2, transpose, lup
|
||||
from mmath/hvector import buildHVector, newHVector, `-`, abs, norm
|
||||
import unittest
|
||||
|
||||
suite "Nim linear algebra library":
|
||||
let test_matrix = newHMatrix[float32](3, 3, [1.0f32, 2f32, 3f32, 4f32, 5f32, 6f32, 8f32, 7f32, 9f32])
|
||||
|
||||
test "+":
|
||||
let mtx1 = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
|
||||
let mtx2 = newHMatrix[int8](3, 3, [-1i8,-2,-3,-4,-5,-6,-7,-8,-9])
|
||||
check(mtx1 + mtx2 == newHMatrix[int8](3,3))
|
||||
|
||||
test "-":
|
||||
let mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
|
||||
check(mtx - mtx == newHMatrix[int8](3, 3))
|
||||
|
||||
test "+=":
|
||||
var mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
|
||||
mtx += mtx
|
||||
check(mtx == newHMatrix[int8](3, 3, [2i8,4,6,8,10,12,14,16,18]))
|
||||
|
||||
test "-=":
|
||||
var mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
|
||||
mtx -= mtx
|
||||
check(mtx == newHMatrix[int8](3, 3))
|
||||
|
||||
test "*":
|
||||
block:
|
||||
let mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,8,7,9])
|
||||
let res = mtx * buildHVector[int8](1i8, 2, 3)
|
||||
check(res == buildHVector[int8](14i8, 32, 49))
|
||||
block:
|
||||
let mtx = newHMatrix[int](2, 5, [1,4,8,2,5,7,3,6,9,0])
|
||||
let res = mtx * mtx.transpose()
|
||||
check(res == newHMatrix[int](2, 2, [110,85,85,175]))
|
||||
block:
|
||||
let mtx = newHMatrix[int](2, 5, [1,4,8,2,5,7,3,6,9,0])
|
||||
let res = mtx.transpose() * mtx
|
||||
check(res == newHMatrix[int](5, 5, [
|
||||
50, 25, 50, 65, 5,
|
||||
25, 25, 50, 35, 20,
|
||||
50, 50, 100, 70, 40,
|
||||
65, 35, 70, 85, 10,
|
||||
5, 20, 40, 10, 25]))
|
||||
|
||||
test "Determinant":
|
||||
check(test_matrix.det() == -9f32)
|
||||
|
||||
test "LU Determinant":
|
||||
check(test_matrix.lu_det() == -9f32)
|
||||
|
||||
test "Inverse":
|
||||
let err = test_matrix * test_matrix.invert() - identity[float32](3)
|
||||
check(err.norm2() < 1e-5)
|
||||
|
||||
test "LU decomposition":
|
||||
var rng = initRand(101325)
|
||||
var arr : array[0..(25 - 1), float32]
|
||||
for i in 0..<arr.len:
|
||||
arr[i] = rng.rand(-4f32..4f32)
|
||||
let mtx = newHMatrix[float32](5, 5, arr)
|
||||
var lu = mtx.clone()
|
||||
let pivot = lu.lup()
|
||||
let l = lu.tril(1.0)
|
||||
let u = lu.triu()
|
||||
let err = pivot * mtx - (l * u)
|
||||
check(err.norm2() < 1e-5)
|
||||
|
||||
test "Linear system solve":
|
||||
var rng = initRand(101325)
|
||||
var arr : array[0..(100 * 100 - 1), float64]
|
||||
for i in 0..<arr.len:
|
||||
arr[i] = rng.rand(-100f32..100f32)
|
||||
var mtx = newHMatrix[float64](100, 100, arr)
|
||||
var lu = mtx.clone()
|
||||
let pivot = lu.lup()
|
||||
var b = newHVector[float64](100)
|
||||
for i in 0...len(b):
|
||||
b[i] = float64(rng.rand(b.len))
|
||||
|
||||
let x = lu.lu_solve(b, pivot)
|
||||
let error = (mtx * x) - b
|
||||
check(error.norm() < 1e-5)
|
||||
|
||||
test "triu":
|
||||
let mtx = newHMatrix[int8](3, 3, [1i8,2,3,4,5,6,7,8,9])
|
||||
block:
|
||||
let upper = mtx.triu()
|
||||
check(upper == newHMatrix[int8](3, 3, [1i8,2,3,0,5,6,0,0,9]))
|
||||
block:
|
||||
let upper = mtx.triu(1)
|
||||
check(upper == newHMatrix[int8](3, 3, [1i8,2,3,0,1,6,0,0,1]))
|
||||
|
||||
test "tril":
|
||||
let mtx = newHMatrix[uint](3, 3, [1u,2,3,4,5,6,7,8,9])
|
||||
block:
|
||||
let lower = mtx.tril()
|
||||
check(lower == newHMatrix[uint](3, 3, [1u,0,0,4,5,0,7,8,9]))
|
||||
block:
|
||||
let lower = mtx.tril(1)
|
||||
check(lower == newHMatrix[uint](3, 3, [1u,0,0,4,1,0,7,8,1]))
|
||||
|
||||
test "transpose":
|
||||
block:
|
||||
let mtx = newHMatrix[int](3, 3, [1, 2, 3, 4, 5, 6, 8, 7, 9])
|
||||
let xpose = newHMatrix[int](3, 3, [1,4,8,2,5,7,3,6,9])
|
||||
check(mtx.transpose() == xpose)
|
||||
block:
|
||||
let mtx = newHMatrix[int](2, 5, [1,4,8,2,5,7,3,6,9,0])
|
||||
let xpose = newHMatrix[int](5, 2, [1,7,4,3,8,6,2,9,5,0])
|
||||
check(mtx.transpose() == xpose)
|
||||
117
tests/test_smatrix.nim
Normal file
117
tests/test_smatrix.nim
Normal file
@@ -0,0 +1,117 @@
|
||||
from random import initRand, rand
|
||||
from nwo/utils import `...`
|
||||
from mmath/smatrix import det, lu, lup, lu_det, lup, det, from_pivot, invert, identity,
|
||||
newSMatrix, SMatrix, clone, tril, triu, lu_solve, `$`, `*`, `-`, `+`, `-=`, `+=`, `==`, `*`, norm2, newSMatrixFromArray,
|
||||
gauss_jordan_high, gauss_jordan_low, transpose
|
||||
from mmath/svector import buildSVector, `-`, abs, norm, Svector
|
||||
from mmath/pivot import SingularMatrixError
|
||||
import unittest
|
||||
|
||||
suite "Nim linear algebra library":
|
||||
let test_matrix = newSMatrixFromArray[3, 3, float32]([1.0f32, 2f32, 3f32, 4f32, 5f32, 6f32, 8f32, 7f32, 9f32])
|
||||
|
||||
test "+":
|
||||
let mtx1 = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
|
||||
let mtx2 = newSMatrixFromArray[3, 3, int8]([-1i8,-2,-3,-4,-5,-6,-7,-8,-9])
|
||||
check(mtx1 + mtx2 == SMatrix[3, 3, int8]())
|
||||
|
||||
test "-":
|
||||
let mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
|
||||
check(mtx - mtx == SMatrix[3, 3, int8]())
|
||||
|
||||
test "+=":
|
||||
var mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
|
||||
mtx += mtx
|
||||
check(mtx == newSMatrixFromArray[3, 3, int8]([2i8,4,6,8,10,12,14,16,18]))
|
||||
|
||||
test "-=":
|
||||
var mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
|
||||
mtx -= mtx
|
||||
check(mtx == SMatrix[3, 3, int8]())
|
||||
|
||||
test "*":
|
||||
block:
|
||||
let mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,8,7,9])
|
||||
let res = mtx * buildSVector[3, int8](1i8, 2, 3)
|
||||
check(res == buildSVector[3, int8](14i8, 32, 49))
|
||||
block:
|
||||
let mtx = newSMatrixFromArray[2, 5, int]([1,4,8,2,5,7,3,6,9,0])
|
||||
let res = mtx * mtx.transpose()
|
||||
check(res == newSMatrixFromArray[2, 2, int]([110,85,85,175]))
|
||||
block:
|
||||
let mtx = newSMatrixFromArray[2, 5, int]([1,4,8,2,5,7,3,6,9,0])
|
||||
let res = mtx.transpose() * mtx
|
||||
check(res == newSMatrixFromArray[5, 5, int]([
|
||||
50, 25, 50, 65, 5,
|
||||
25, 25, 50, 35, 20,
|
||||
50, 50, 100, 70, 40,
|
||||
65, 35, 70, 85, 10,
|
||||
5, 20, 40, 10, 25]))
|
||||
|
||||
|
||||
test "Determinant":
|
||||
check(test_matrix.det() == -9f32)
|
||||
|
||||
test "LU Determinant":
|
||||
check(test_matrix.lu_det() == -9f32)
|
||||
|
||||
test "Inverse":
|
||||
let err = test_matrix * test_matrix.invert() - identity[3, float32]()
|
||||
check(err.norm2() < 1e-5)
|
||||
|
||||
test "LU decomposition":
|
||||
var rng = initRand(101325)
|
||||
var arr : array[0..(25 - 1), float32]
|
||||
for i in 0..<arr.len:
|
||||
arr[i] = rng.rand(-4f32..4f32)
|
||||
let mtx = newSMatrixFromArray[5, 5, float32](arr)
|
||||
var lu = mtx.clone()
|
||||
let pivot = lu.lup()
|
||||
let l = lu.tril(1.0)
|
||||
let u = lu.triu()
|
||||
let err = pivot * mtx - (l * u)
|
||||
check(err.norm2() < 1e-5)
|
||||
|
||||
test "Linear system solve":
|
||||
var rng = initRand(101325)
|
||||
var arr : array[0..(100 * 100 - 1), float64]
|
||||
for i in 0..<arr.len:
|
||||
arr[i] = rng.rand(-100f32..100f32)
|
||||
var mtx = newSMatrixFromArray[100, 100, float64](arr)
|
||||
var lu = mtx.clone()
|
||||
let pivot = lu.lup()
|
||||
var b : SVector[100, float64]
|
||||
for i in 0...len(b):
|
||||
b[i] = float64(rng.rand(b.len))
|
||||
|
||||
let x = lu.lu_solve(b, pivot)
|
||||
let error = (mtx * x) - b
|
||||
check(error.norm() < 1e-5)
|
||||
|
||||
test "triu":
|
||||
let mtx = newSMatrixFromArray[3, 3, int8]([1i8,2,3,4,5,6,7,8,9])
|
||||
block:
|
||||
let upper = mtx.triu()
|
||||
check(upper == newSMatrixFromArray[3, 3, int8]([1i8,2,3,0,5,6,0,0,9]))
|
||||
block:
|
||||
let upper = mtx.triu(1)
|
||||
check(upper == newSMatrixFromArray[3, 3, int8]([1i8,2,3,0,1,6,0,0,1]))
|
||||
|
||||
test "tril":
|
||||
let mtx = newSMatrixFromArray[3, 3, uint]([1u,2,3,4,5,6,7,8,9])
|
||||
block:
|
||||
let lower = mtx.tril()
|
||||
check(lower == newSMatrixFromArray[3, 3, uint]([1u,0,0,4,5,0,7,8,9]))
|
||||
block:
|
||||
let lower = mtx.tril(1)
|
||||
check(lower == newSMatrixFromArray[3, 3, uint]([1u,0,0,4,1,0,7,8,1]))
|
||||
|
||||
test "transpose":
|
||||
block:
|
||||
let mtx = newSMatrixFromArray[3, 3, int]([1, 2, 3, 4, 5, 6, 8, 7, 9])
|
||||
let xpose = newSMatrixFromArray[3, 3, int]([1,4,8,2,5,7,3,6,9])
|
||||
check(mtx.transpose() == xpose)
|
||||
block:
|
||||
let mtx = newSMatrixFromArray[2, 5, int]([1,4,8,2,5,7,3,6,9,0])
|
||||
let xpose = newSMatrixFromArray[5, 2, int]([1,7,4,3,8,6,2,9,5,0])
|
||||
check(mtx.transpose() == xpose)
|
||||
72
vector.nim
72
vector.nim
@@ -1,72 +0,0 @@
|
||||
from oomacro import class
|
||||
from utils import `...`
|
||||
from sequtils import newSeqWith
|
||||
from math import sqrt
|
||||
|
||||
type Vector*[T] = seq[T]
|
||||
|
||||
#SVector*[S : static[int], T] = array[1..S,T]
|
||||
|
||||
proc newVector*[T](size : int, init: T=0) : Vector[T] =
|
||||
result = newSeq[T](size)
|
||||
for i in 0...len(result):
|
||||
result[i] = init
|
||||
|
||||
proc `+`*[T](v1 : Vector[T], v2:Vector[T]) : Vector[T] =
|
||||
result = newVector[T](len(v1))
|
||||
for i in 0...len(v1):
|
||||
result[i] = v1[i] + v2[i]
|
||||
|
||||
proc `-`*[T](v1 : Vector[T], v2:Vector[T]) : Vector[T] =
|
||||
result = newVector[T](len(v1))
|
||||
for i in 0...len(v1):
|
||||
result[i] = v1[i] - v2[i]
|
||||
|
||||
proc `*`*[T](v1 : Vector[T], v2:Vector[T]) : T =
|
||||
result = 0
|
||||
for i in 0...len(v1):
|
||||
result += v1[i] * v2[i]
|
||||
|
||||
proc `+=`*[T](v1 : var Vector[T], v2 : Vector[T]) =
|
||||
for i in 0...len(v1):
|
||||
v1[i] += v2[i]
|
||||
|
||||
proc `-=`*[T](v1 : var Vector[T], v2 : Vector[T]) =
|
||||
for i in 0...len(v1):
|
||||
v1[i] -= v2[i]
|
||||
|
||||
proc `+=`*[T](v : var Vector[T], value : T) = v.add(value)
|
||||
|
||||
proc createVector*[T](elems : varargs[T]) : Vector[T] =
|
||||
result = newSeq[T]()
|
||||
for elem in items(elems):
|
||||
result += elem
|
||||
|
||||
proc norm*[T](v : Vector[T]) : T =
|
||||
for value in v:
|
||||
result += v * v
|
||||
|
||||
proc abs*[T](v : Vector[T]) : T =
|
||||
return math.sqrt(v.norm)
|
||||
|
||||
# let vec = createVector(10,1)
|
||||
|
||||
# echo $vec
|
||||
# var vec2 = newVector[int](10, 1)
|
||||
|
||||
# echo vec2
|
||||
|
||||
# type Foo[T, S] = object of RootObj
|
||||
# data : S
|
||||
|
||||
# proc `+`[T,S](v1 : Foo[T,S], v2 : Foo[T,S]) : Foo[T,S] =
|
||||
# result = Foo[T,S]()
|
||||
# for i in range(len(v1.data)):
|
||||
# result[i] = v1.data[i] + v2.data[i]
|
||||
|
||||
|
||||
# type Bar[T] = Foo[T, array[3,T]]
|
||||
|
||||
# var bar = Bar[int]()
|
||||
# bar.data[1] = 4
|
||||
# bar.data[1] = 4
|
||||
Reference in New Issue
Block a user