[WIP] Float types as parametrics

[mratsim]

This is work in progress to solve: #5

This feature is key to implement machine learning algorithms without code duplication for float32, float64 and float (pointer-sized, 32 bits on 32-bit, 64 bits on 64-bit machine)

A lesser alternative is exporting AnyMatrix = Matrix32 or Matrix64 or DMatrix32 or DMatrix64 but the proc will lose dimension information.

Furthermore should Nim or a library implement half-precision float16 in the future, supporting it would be trivial. (float16 allows 2x TFLOPS on some hardware).

[andreaferretti]

Thank you for taking time to do this refactoring. There are, unfortunately, two things that I do not like about it.

For one thing, you removed all static[int] constraints. I am not sure why this is needed, and I think it makes things less clear.

Second, this makes DVector andDMatrix into generic types. Unfortunately, their type parameter can only be specialized to float32 or float64 (because BLAS only supports those), while the type may suggest that other combinations - say DVector[int] are supported. Not only that, but now this breaks compatibility with actual generic implementations (because the compiler will not be able to decide which one is more specific). There is one in Emmy, but that is just an example. Still, the point stands that it makes sense to have matrices over something other than floating point numbers with a generic implementation, and two generic implementations will not coexist easily.

As such, I am not sure I am going to merge this. On the other hand, if this massively simplifies the implementation by reducing the duplication, it may be worth it.

Let's discuss it a bit.

By the way, I have a prototype neural network implementation based on this library, and it has its fair share of duplication. I am going to open it in the next days, so may be it will be a good starting point to see the effect of this change in concrete

[andreaferretti]

By the way, I see that you moved all constraints to the types definitions. I had missed that, so part of what I said above does not make sense.

For some reason, I always found issues in the compiler when putting type constraints there - this is why I was putting them in functions. I will experiment a bit with your new version, and let you know any other comment I have

[mratsim]

Just to make sure, the type definition currently is:

type
  Vector*[N: static[int], T: SomeReal] = ref array[N, T]
  Matrix*[M, N: static[int], T: SomeReal] = object
    order: OrderType
    data: ref array[N * M, T]

Regarding your first point, static[int] should be enforced on all proc due to being enforced at the type level.

For your second point, T is declared with T: SomeReal so it's specialized over SomeReal = float | float32 | float64, so there will be no issues regarding BLAS. SomeReal doc.

Is that OK for both?

Regarding the generalization to DVector[int], I see your point. After refactoring, I will explore using is to help the compiler choose between the generic or specialized proc. is operator doc

Is operator

The is operator checks for type equivalence at compile time. It is therefore very useful for type specialization within generic code

[mratsim]

Bad news, I realized that when launching tests I was importing the working linalg library instead of my modifications.

It seemed to good to be true to pass all the tests whatever I did. So I will probably have to rework my commits from scratch.

[mratsim]

I've reset the branch to rework from scratch while making sure I load the WIP version with the test suite.

I'll keep the static keyword on procs to avoid introducing regression.

With the current very minimal modifications, I fail the following tests:

funcs.nim

funcs.nim(26, 18) Error: type mismatch: got (Matrix64[2, 2], Matrix64[2, 1])
but expected one of:
proc solve[M, N: static[int]](a: Matrix32[M, M]; b: Matrix32[M, N]): Matrix32[M, N]
proc solve[M: static[int]](a: Matrix32[M, M]; b: Vector32[M]): Vector32[M]
proc solve[M, N: static[int]](a: Matrix64[M, M]; b: Matrix64[M, N]): Matrix64[M, N]
proc solve[M: static[int]](a: Matrix64[M, M]; b: Vector64[M]): Vector64[M]

ops and row_major_ops but for Error: ordinal type expected (#13)

[mratsim]

Regarding the solve issue, I've reproduced a minimum test case. It seems like compiler bug nim-lang/Nim#4554:

import linalg

let a = matrix(@[
  @[3.0, 1.0],
  @[1.0, -2.0]
], 2, 2)
let b = matrix(@[@[1.0],
                 @[0.0]],
                2, 1)

#let x = a.solve(b)
proc test*(a: Matrix64, b: Matrix64) =
  echo a
  echo b

test(a,a) #works
test(b,b) #works
test(a,b) #fails

Error: type mismatch: got (Matrix64[2, 2], Matrix64[2, 1])
but expected one of:
proc test[Matrix64](a: Matrix64; b: Matrix64)

[andreaferretti]

Oh, I am sorry to hear that, but in fact I was surprised that the compiler allowed type constraints in the type declaration section. I think you should be good just by running nimble test on your cloned copy. Unless I have messed something up, it should allow you to make modifications and test them out of the box.

Regarding the compiler issues with static - it is unfortunate that a regression was introduced. I hope that after this PR gets merged, stabilizing concepts, @zah will tackle static, which is the last thing that regularly gives me issues! :-)

Let me know if I can be of any help

[mratsim]

For now I'm stuck as all types of workarounds for the solve proc do not work:
Reproducing with a minimal function actual show that I will have more problems down the line

import linalg

let a = matrix(@[
  @[3.0, 1.0],
  @[1.0, -2.0]
], 2, 2)
let b = matrix(@[@[1.0],
                 @[0.0]],
                2, 1)

#let x = a.solve(b)
proc test*[M,N: static[int]](a: Matrix64[M,N], b: Matrix64[M,N]) =
  echo a
  echo b

test(a,a) #works
test(b,b) #works
#test(a,b) # fails type mismatch
#Error: type mismatch: got (Matrix64[2, 2], Matrix64[2, 1])
#but expected one of:
#proc test[M, N](a: Matrix64[M, M]; b: Matrix64[M, N])
echo "Test #1 DONE"

proc test2*[M,N: static[int]](a: Matrix64[M,M], b: Matrix64[M,N]) =
    echo a
    echo b

test2(a,a) #works
#test2(b,b) #b is [2,1], dim error should be catched by Nim but it's catched by clang
#Error: execution of an external compiler program 'clang -c  -w  -I/usr/local/Cellar/nim/HEAD-875e344/nim/lib -o nimcache/linalg_playground.o nimcache/linalg_playground.c' failed with exit code: 1

#nimcache/linalg_playground.c:495:31: error: passing 'Matrix_zhuhOuhVfbLJLsmoc8sThg' (aka 'struct Matrix_zhuhOuhVfbLJLsmoc8sThg') to parameter of incompatible type 'Matrix_4UORjnPzvIMBohnnojhPgw' (aka 'struct Matrix_4UORjnPzvIMBohnnojhPgw')
#        test2_A4PrUCOWHf2urCkKnNhZzg(b_ONM9cc62rV8Rg9cSV3Ca9c3lA, b_ONM9cc62rV8Rg9cSV3Ca9c3lA);
#                                     ^~~~~~~~~~~~~~~~~~~~~~~~~~~
#nimcache/linalg_playground.c:400:77: note: passing argument to parameter 'a' here
#N_NIMCALL(void, test2_A4PrUCOWHf2urCkKnNhZzg)(Matrix_4UORjnPzvIMBohnnojhPgw a, Matrix_zhuhOuhVfbLJLsmoc8sThg b) {
#                                                                            ^
#1 error generated.
echo "Test #2 DONE"

proc test3*[M,N: static[int], T](a: Matrix[M,N,T], b: Matrix[M,N,T]) =
  echo a
  echo b

test3(a,a) #fails
#Error: ambiguous call; both linalg.$(m: Matrix32[$.M, $.N]) and linalg.$(m: Matrix64[$.M, $.N]) match for: (Matrix[2, 2])
#test3(a,b) #same failure as just above
echo "Test #3 DONE"

For the last part, it seems like $ does not detect the float type at compile time.
I tried changing $ in display.nim to

proc `$`*[M, N: static[int]](m: Matrix32[M,N] or Matrix64[M,N]): string =
  result = "[ "
  for i in 0 .. < M - 1:
    result &= toStringHorizontal(m.row(i)) & "\n  "
  result &= toStringHorizontal(m.row(M - 1)) & " ]"

Alas I get caught in overloading resolution issues with system.$

Error: ambiguous call; both system.$(x: T) and linalg.$(m: Matrix[$.M, $.N, $.T] or Matrix32[$.M, $.N]) match for: (Matrix64[2, 2])

Resolution order

1. Exact match: a and f are of the same type.
2. Literal match: a is an integer literal of value v and f is a signed or unsigned integer type and v is in f's range. Or: a is a floating point literal of value v and f is a floating point type and v is in f's range.
3. Generic match: f is a generic type and a matches, for instance a is int and f is a generic (constrained) parameter type (like in [T] or [T: int|char].

I thought nim-lang/Nim#1991 regarding type class overloading and fixed 16 hours ago would help but unfortunately no.

When trying to completely remove Matrix32 and Matrix64 and use the new underlying Matrix, I get hit with static bug:

proc `$`*[M, N: static[int], T](m: Matrix[M, N, T]): string =
  result = "[ "
  for i in 0 .. < M - 1:
    result &= toStringHorizontal(m.row(i)) & "\n  "
  result &= toStringHorizontal(m.row(M - 1)) & " ]"

Error: cannot instantiate: 'M'

I have yet to try the syntax mentioned in the bug nim-lang/Nim#1991 but it will probably need generic generic parameters https://github.com/nim-lang/Nim/issues/3856

type B = object of A

proc p[X](x: X): string = "generic"
proc p[X: A](x: X): string = "specialized"

[mratsim]

Opened nim-lang/Nim#5641, nim-lang/Nim#5643, nim-lang/Nim#5644

Edit nim-lang/Nim#5645 as well

[andreaferretti]

Thank you!

[mratsim]

Following the resolution of almost all Nim related issues, I have now a minimal working rewrite.
It compiles and passes all the tests, except the ones that uses seq[float] instead of a DVector. Reported here: #15

Now all remaining procs can be converted one by one.

[andreaferretti]

It is great to hear that this is now working! I am not sure what problems seq[float] vs DVector might cause, because they are in fact the same type!

[andreaferretti]

As mentioned in the README, an expanded rewrite is here