planetis-m / manu Goto Github PK

Nim MAtrix NUmeric package

License: MIT License

Nim 100.00%

linear-algebra scientific nim decomposition matrix-decompositions matrix-functions nim-library determinants simultaneous-linear-equations matrices

manu's Introduction

Manu — Nim Matrix Numeric library

Manu is a pure Nim library, no external dependencies to BLAS frameworks. Supports constructing and manipulating only real, dense matrices. It started as a port of JAMA library, and is adapted to Nim programming paradigm and specific performace considerations.

What is supported:

Compute solutions of simultaneous linear equations, determinants, inverses and other matrix functions.
Generics allow matrices of SomeFloat only.
Arithmetic operators are overloaded to support matrices.
- Broadcast scalars, column and row vectors to work with matrices.
Destructors, with sink annotations, copies can be avoided in some cases.

API documentation

Examples

In the examples directory you will find the following:

showcasing what can already be done.

example2.nim

  import manu

  # Solve a linear system A x = b and compute the residual norm, ||b - A x||.
  let vals = @[@[1.0, 2, 3], @[4.0, 5, 6], @[7.0, 8, 10]]
  let A = matrix(vals)
  let b = randMatrix64(3, 1)
  let x = A.solve(b)
  let r = A * x - b
  let rnorm = r.normInf()
  echo("x =\n", x)
  echo("residual norm = ", rnorm)

Output:

x =
⎡-918.9217543597e-3⎤
⎢      2.1952979104⎥
⎣     -1.0796593055⎦
residual norm = 1.554312234475219e-15

Matrix decompositions

Five matrix decompositions are used to compute solutions of simultaneous linear equations, determinants, inverses and other matrix functions. Theses are:

Cholesky Decomposition of symmetric, positive definite matrices
LU Decomposition (Gaussian elimination) of rectangular matrices
QR Decomposition of rectangular matrices
Eigenvalue Decomposition of both symmetric and nonsymmetric square matrices
Singular Value Decomposition of rectangular matrices

Broadcasting

It is implemented with the help of two distinct types RowVector[T] and ColVector[T]. You can cast any compatible matrix to these and when performing matrix operations, it will be broadcasted to the correct dimensions:

var a = matrix(1, 5, 2.0)
let b = ones64(2, 1)
echo ColVector64(b) + RowVector64(a)
echo 2.0 + a # matrix-scalar ops are implicit

Results in:

⎡3  3  3  3  3⎤
⎣3  3  3  3  3⎦
⎡4  4  4  4  4⎤

If the matrices are not broadcastable an AssertionDefect will be thrown at runtime.

The correct paradigm of usage is to first initialize a matrix, i.e let a = ones64(1, 5) and cast it to RowVector64 where broadcasting is needed: RowVector64(a) + zeros64(5, 5). This system is designed to be more explicit, and since it is type-checked, work well with sink optimizations.

Feature improvements / contributions

Add more tests
Incorporate usefull additions from Apache Commons Math

License

This library is distributed under the MIT license.

manu's People

Contributors

Stargazers

Watchers

Forkers

refaqtor dseeni

manu's Issues

Do inline procs prevent sink optimization?

outofbounds read on matrix seq[seq[T]] constructor

Non-Commutative operators seem wrong, arguments are reversed

Make it compile time

Currently arc/orc is broken with the combination of top-level variables and closures. I would like to contribute a version that uses macrocache to generate a simple case statement in compile-time. I am not sure how receptive the maintainers are to this idea? And even if it would worth it in the end.

matrix-scalar multiplication fails

Examples 1 and 2 fail to compile


# Solve a linear system A x = b and compute the residual norm, ||b - A x||.
let vals = @[@[1.0, 2.0, 3.0], @[4.0, 5.0, 6.0], @[7.0, 8.0, 10.0]]
let A = matrix(vals)
let b = randMatrix(3, 1)
let x = A.solve(b)
let r = A * x - b
let rnorm = r.normInf()
echo("x =\n", x)
echo("residual norm = ", rnorm)

#[
manu_ex2.nim(6, 19) template/generic instantiation of `randMatrix` from here
.nimble/pkgs/manu-2.2.0/manu/matrix.nim(157, 18) Error: cannot instantiate: 'T'
]#

and


template ff(f: float, prec: int = 3): string =
   formatFloat(f, ffDecimal, prec)

proc magic(n: int): Matrix =
   var m = matrix(n, n)
   # Odd order
   if n mod 2 == 1:
      let a = (n + 1) div 2
      let b = n + 1
      for j in 0 ..< n:
         for i in 0 ..< n:
            m[i, j] = float(n * ((i + j + a) mod n) + ((i + 2 * j + b) mod n) + 1)
   # Doubly Even Order
   elif n div 4 == 0:
      for j in 0 ..< n:
         for i in 0  ..< n:
            if (i + 1) div 2 mod 2 == (j + 1) div 2 mod 2:
               m[i, j] = float(n * n - n * i - j)
            else:
               m[i, j] = float(n * i + j + 1)
   # Singly Even Order
   else:
      let p = n div 2
      let k = (n - 2) div 4
      let a = magic(p)
      for j in 0 ..< p:
         for i in 0 ..< p:
            let aij = a[i, j]
            m[i, j] = aij
            m[i, j + p] = aij + float(2 * p * p)
            m[i + p, j] = aij + float(3 * p * p)
            m[i + p, j + p] = aij + float(p * p)
      for i in 0 ..< p:
         for j in 0 ..< k:
            swap(m[i, j], m[i + p, j])
         for j in n - k + 1 ..< n:
            swap(m[i, j], m[i + p, j])
      swap(m[k, 0], m[k + p, 0])
      swap(m[k, k], m[k + p, k])
   return m

proc main() =
   # Tests LU, QR, SVD and symmetric Eig decompositions.
   # 
   #   n       = order of magic square.
   #   trace   = diagonal sum, should be the magic sum, (n^3 + n)/2.
   #   max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
   #   rank    = linear algebraic rank,
   #             should equal n if n is odd, be less than n if n is even.
   #   cond    = L_2 condition number, ratio of singular values.
   #   lu_res  = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
   #   qr_res  = test of QR factorization, norm1(Q*R-A)/(n*eps).

   echo("    Test of Matrix object, using magic squares.")
   echo("    See MagicSquareExample.main() for an explanation.")
   echo("      n     trace       max_eig   rank        cond      lu_res      qr_res\n")

   let start_time = epochTime()
   let eps = pow(2.0, -52.0)
   var buf = ""
   for n in 3 .. 32:
      buf.setLen(0)
      buf.add(align($n, 7))

      let M = magic(n)
      let t = int(M.trace())
      buf.add(align($t, 10))

      let E =
         eig((M + M.transpose()) * 0.5)
      let d = E.getRealEigenvalues()
      buf.add(align(ff(d[n-1]), 14))

      let r = M.rank()
      buf.add(align($r, 7))

      let c = M.cond()
      if c < 1.0 / eps:
         buf.add(align(ff(c), 12))
      else:
         buf.add("         Inf")

      let LU = lu(M)
      let L = LU.getL()
      let U = LU.getU()
      let p = LU.getPivot()
      var R = L * U - M[p, 0 .. n-1]
      var res = R.norm1() / (float(n) * eps)
      buf.add(align(ff(res), 12))

      let QR = qr(M)
      let Q = QR.getQ()
      R = QR.getR()
      R = Q * R - M
      res = R.norm1() / (float(n) * eps)
      buf.add(align(ff(res), 12))

      echo buf

   let stop_time = epochTime()
   let etime = (stop_time - start_time) / 1000.0
   echo("\nElapsed Time = ", align(ff(etime), 12), " seconds")
   echo("Adios")

main()

#[
manu_ex1.nim(67, 20) template/generic instantiation of `magic` from here
manu_ex1.nim(6, 21) Error: cannot instantiate 'Matrix[T: SomeFloat]' inside of type definition: 'magic'; Maybe generic arguments are missing?
]#

nim --version :
Nim Compiler Version 1.5.1 [Linux: amd64]
Compiled at 2021-01-07
Copyright (c) 2006-2020 by Andreas Rumpf

git hash: 89a21e4ec71e705833d2aacd069e291cf41a19c6
active boot switches: -d:release

 PASS: [11/35] kdtree c                                                     ( 6.04 sec)
  FAIL: [12/35] manu c                                                       ( 7.09 sec)
  Test "manu" in category "nimble-packages"
  Failure: reBuildFailed
  nimble test
    Executing task test in /home/runner/work/Nim/Nim/pkgstemp/manu/manu.nimble
  sh: 1: nake: not found
  stack trace: (most recent call last)
  /tmp/nimblecache-946033505/nimscriptapi_2731018556.nim(187, 16)
  /home/runner/work/Nim/Nim/pkgstemp/manu/manu.nimble(12, 8) testTask
  /home/runner/work/Nim/Nim/lib/system/nimscript.nim(273, 7) exec
  /home/runner/work/Nim/Nim/lib/system/nimscript.nim(273, 7) Error: unhandled exception: FAILED: nake test [OSError]
       Error: Exception raised during nimble script execution
  
  PASS: [13/35] msgpack4nim c                                                ( 4.06 sec)