numerics in Lisp

[grapesmoker]

It seems like achieving the aims of this project will require some kind of reasonably efficient numerical tools. My day job is working in numpy/pandas/dask, so I'm pretty familiar with that universe of tools, but I don't know what people use these days in Lisp, if anything. I took a brief look around a while back and I came up with:

• GSLL - bindings to the GNU scientific library
• MAGICL - linear algebra package that wraps BLAS

I'm sure there are others I don't know about, but it's hard to find out what they might be. I saw on the CL subreddit some months back that Marco Heisig had done some work on autovectorizing Lisp code that was very promising from a performance standpoint but I'm not sure what the maturity of that project is (and I think it was only for SBCL).

Maybe we don't need to reinvent the numpy wheel (haha) to make this specific project work, but I think it would be cool if there was something that was as ergonomic in CL as numpy is in Python and offered comparable performance. It would certainly help in the implementation of various geometrical algorithms. Would love to hear what others think about this.

[awolven]

I've been working on an interface to Eigen, if that helps.

…

On Wed, Aug 24, 2022 at 7:40 PM Jerry Vinokurov ***@***.***> wrote: It seems like achieving the aims of this project will require some kind of reasonably efficient numerical tools. My day job is working in numpy/pandas/dask, so I'm pretty familiar with that universe of tools, but I don't know what people use these days in Lisp, if anything. I took a brief look around a while back and I came up with: - GSLL - bindings to the GNU scientific library - MAGICL - linear algebra package that wraps BLAS I'm sure there are others I don't know about, but it's hard to find out what they might be. I saw on the CL subreddit some months back that Marco Heisig had done some work on autovectorizing Lisp code that was very promising from a performance standpoint but I'm not sure what the maturity of that project is (and I think it was only for SBCL). Maybe we don't need to reinvent the numpy wheel (haha) to make this specific project work, but I think it would be cool if there was something that was as ergonomic in CL as numpy is in Python and offered comparable performance. It would certainly help in the implementation of various geometrical algorithms. Would love to hear what others think about this. — Reply to this email directly, view it on GitHub <#41>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABGMMJFAZTHA2FIRDQYOG3V226JDANCNFSM57RHDUMA> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[kaveh808]

@awolven Eigen came up in a discussion. How far along are you on your interface?

[kaveh808]

I see three areas in which numerics can be helpful.

1. SIMD for point and matrix operations. I understand SBCL does SIMD now, but I don't know how to use it.
2. Threads for parallel work on large 1D and 2D point arrays across multiple cores.
3. GPU acceleration of compute operations.

[foretspaisibles]

The text Fast Floating-Point Pro cessing in Common Lisp by Fateman, Broughan, Willcock and Rettig discusses how to efficiently program numerics in Lisp and illustrates declarations as well as macros that can make the declarations less tedious.

[kaveh808]

This was also brought to my attention if we want to go the FFI route:

https://github.com/digikar99/bmas
https://github.com/digikar99/cl-bmas

[JMC-design]

Also things like this https://github.com/takagi/avm , i think there's an cl-opencl out there. And for simple advantage of simd in sbcl https://github.com/marcoheisig/Loopus

[JMC-design]

https://github.com/melisgl/mgl-mat

[kaveh808]

I have been advised that for our needs SIMD is more appropriate than CUDA.

https://www.reddit.com/r/GraphicsProgramming/comments/x1rwoo/how_widespread_is_cuda_availability/?

My preference is still a purely CL option, but we could go with an FFI lib if we have to.

[ghost]

From the Reddit article:

Eigen uses vectorisation under the hood, is easy to install (header only) and blazing fast

It also uses C++ template metaprogramming under the hood--this is not a compliment--you don't want to look under its hood except maybe to port it to CL. Compilation of Eigen is "blazing", in that it requires your CPU fan to spin at full RPM.

My preference is still a purely CL option

I agree. Using Eigen would be expedient and maybe a reasonable choice, but Common Lisp is a flexible language with well-specified compile-time metaprogramming facilities. Compiling an Eigen-enabled program has crashed my laptop (admittedly low-end) more than once--CL never.

[awolven]

The eigen interface which I am using basically does what I want it to do, but could be extended fairly easily. It currently needs some work as the original author did not provide finalizations on matrices.

…

On Tue, Aug 30, 2022 at 6:48 PM Kaveh Kardan ***@***.***> wrote: @awolven <https://github.com/awolven> Eigen came up in a discussion. How far along are you on your interface? — Reply to this email directly, view it on GitHub <#41 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AABGMMOIDDHTRACRMVCHCXDV32MT3ANCNFSM57RHDUMA> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[kaveh808]

Nice example:

https://gist.github.com/digikar99/13ed3868993ca8c71c91b50c8022a64d

[JMC-design]

Their survey is quite informative as well.

[JMC-design]

https://www.codeproject.com/Articles/874396/Crunching-Numbers-with-AVX-and-AVX

[digikar99]

I think the choice of tools depends on what the use cases are and what you want to prioritize. If the team wants to learn SIMD, then one could certainly go the sb-simd way, possibly extending it to other platforms (arm64, M1, etc). But if getting things to work ASAP and reasnably-portably is a priority, I think an FFI would be better, largely because C and fortran libraries seem to have been battling the issue of cross-compiler and cross-platform portability-with-performance for decades, and so not using them seems like we are discarding their efforts; and then if it does turn out that it actually hampers extensibility, then perhaps one could do some slight refactoring to take sb-simd into account wherever appropriate.

For the FFI route, if you are working with at least moderately sized vectors and matrices (more than 100 or 1000 elements), magicl seems like a fine choice in terms of performance as well as stability. I am not exactly sure what its status is on other implementations, but I also don't think they are going terribly beyond the ANSI standard. BLAS is also an industry standard that shouldn't be a wrong choice. In fact, for smaller arrays, one could try putting static-dispatch to recompile magicl appropriately. And wherever magicl happens to be incomplete, one could add more functionality to it; their teams seems fairly comfortable with PRs and requests.

I also don't think ANSI CL is sufficient for a good number crunching library that is performant for small arrays as well as large, unless you invent a whole DSL itself (I'm looking at coalton!). In effect, cross-implementation portability gets thrown out of the window, unless the implementation decides to at least support CFFI and CLTL2 (and perhaps CLOSER-MOP as well). Even if you do have CFFI, writing reasonably maintable/modular code requires dispatching on specialized array types. Here, generic functions do not suffice (and I haven't found any good way to use CLOSER-MOP - or even SB-PCL - to achieve the same results on SBCL). Writing a small layer of abstraction for dispatching on types using only lambda-lists with all arguments as required and positions (rather than optional, keyword, and rest) takes less than 300 LOC. But when you put in the code for parsing optional, keyword, and rest, code blows up pretty quickly - and this was the main goal of polymorphic-functions. However that only gets you runtime dispatch. To get compile-time dispatch in order to be performant on small arrays, you need CLTL2, through something like cl-environments and cl-form-types. Without them, I don't think I'd recommend Common Lisp as the "ultimate" language for number crunching. (With them though, I think its potential goes beyond julia, if only we had more people writing lisp :/.)

[kaveh808]

I'm thinking we need our numerics to do the following:

1. fast vector/matrix math for our internal calculations such as procedural modeling (mesh extrusions, heightfields, etc)
2. fast data transfer to our draw routines (OpenGL 4.1, Vulkan, etc.)
3. fast data transfer to the renderer, especially if a real-time renderer such as Filament
4. fast data transfer to external C libs such as physics engines
5. ... anything I'm missing?

Case (1) would benefit from vectorization and concurrency on our part. The others would do their own such things under the hood I imagine.

So it appears to me that most of our needs involve data transfer to the outside (C/C++) world. Which makes me wonder if my preference to keep everything in CL is misguided.

Should we store our vector and math data structures as C arrays or structs with simple CL interfaces?

[ghost]

we should start another discussion about a roadmap

OK, if you want to create a new roadmap discussion, I'll continue to contribute my ideas.

Back on topic, without regard to exactly what parts of that list you hope to accomplish, Common Lisp has some numeric issues which need to be addressed, again for example:

https://www.researchgate.net/publication/339814323_Why_You_Cannot_Yet_Write_an_Interval_Arithmetic_Library_in_Common_Lisp

Which makes me wonder if my preference to keep everything in CL is misguided.
...
Should we store our vector and math data structures as C arrays or structs with simple CL interfaces?

This is a hell of a good point to ponder. The reality is we are running on C machines not Lisp machines. Until we migrate back, we won't be able to "keep everything in CL". Maybe conceding a little to the C model would allow us to keep more in Common Lisp ultimately.

[kaveh808]

Made roadmap discussion.

[digikar99]

Should we store our vector and math data structures as C arrays or structs with simple CL interfaces?

As far as arrays of implementation supported numerics are concerned, CFFI has a cffi-sys:with-pointer-to-vector-data which provides pointers to cl:vector object's data, that I have found useful. Although, I'm not aware about its limitations.

But packed structures is an open issue. I had come across suggestions about structures of arrays rather than arrays of structures, but which of the two suits better would depend on the operations involved!

[digikar99]

https://www.researchgate.net/publication/339814323_Why_You_Cannot_Yet_Write_an_Interval_Arithmetic_Library_in_Common_Lisp

That looks like a particularly useful paper, thanks!

Any idea if the issues will persist even if we let go off portability for this particular aspect and focus on SBCL?

[JMC-design]

@digikar99 thanks for pointing out that cffi bit, I knew how to do it in sbcl but not port ably. Will be replacing my code with that.
On arrays of structs, arrays of arrays etc... I found a lot of concern about it in the past. but not so much in the present.
Just like separate buffers for vertices and indices or one buffer, Carmack says do whatever no big difference.

[mtreis86]

I am not sure if it would fit, but I have been experimenting with and using Coalton and it offers a bunch of things CL is missing in terms of building tools that work with numbers.

[digikar99]

Coalton could be helpful for solving the arrays of structures issue, but I don't think it can by itself resolve the issues about numerics mentioned in the paper above.

But also, could Coalton on SBCL help solve both the issues?

[JMC-design]

perhaps all this is a bit premature. Shouldn't we be looking at a framework wherein these possibilities exist? So coming up with a protocol for what needs to be done. e.g. object->vertex-array, assemble-object-for-display, submit-objet-for-display, translate-object-for-client :gl/:vulkan, etc...
Obviously not those names or breakdown, just trying to give a sense. It all seems a bit premature optimization when we don't know what the bottle necks are. ( i mean, besides legacy gl)

[kaveh808]

I'm in favor of taking small steps. Based on the timings that were done for packing points, looks like the Origin package is a good one to switch to for our vector and matrix code.

Once that is done we can test packing our geometry into vertex arrays by replacing the functions in opengl.lisp.

Then we can do more tests and see how best to vectorize/optimize the code.

[JMC-design]

So I just wanted to be sure you're ok with the not being able to dispatch on them right? well, i guess you can and typecase...
Then I'm guessing we should change all the (p* vector scalar) back to (p:scale vector scalar)? and have p: be a local nickname for origin.vec3?

[kaveh808]

We currently do a ctypecase on p+, p-, p*, p/.

We can leave it as is or switch back to p-scale. Whatever you think best.

[JMC-design]

While I like the genericness of * for everything(use it myself), scale gives us more information when we're skimming the code. and if being pedantic, it is the right term.

[JMC-design]

the p! is so ubiquitous I've just kept it as is.

[JMC-design]

nevermind, I forgot origin only has operations between vectors.
will keep the p+ p- etc..

[JMC-design]

So I've got the basics done over here https://github.com/JMC-design/kons-9/tree/point-origin
I'm running into a lot of The value 1.0000001 is not of type (SINGLE-FLOAT -1.0 1.0) in demo.lisp which is all I've tested so far. Not sure if something should be clamped somewhere or where that expectation of being between -1 and 1 comes from. Any pointers to where that happens would help.
All occurences of the clamp issue are commented with something about clamps

[JMC-design]

nevermind, origin angle doesn't like points being equal

[JMC-design]

all demos run.
major differences p:x p:y p:z accessors, work with setf, p-set! still usuable.
p:length instead of magnitude
p:scale p:* p:+ p:- p:/ where appropriate, p/ p* etc when used with scalars or was just too ambiguous
p:lerp p1 p2 factor at end.
p:random lo hi when two scalars used, all other rand functions still there.
check point-origin.lisp for what stayed, and all changed to use origin where appropriate.

[kaveh808]

Excellent. Send a pull request to main.

[JMC-design]

will do. I'm trying to figure out how to squash the commits because I messed up rebasing against wrong branch so the history might be horrible.
I guess subbing in the matrix stuff should be a different commit.

[kaveh808]

I agree. Let's do the matrix stuff separately.