SciEx - scientific programming for Elixir (based on bindings to rust's ndarray)

Added a lot of functionality today. SciEx now support arrays of floats of dimension up to 6 (SciEx.Float64.Array1, …, SciEx.Float64.Array6). Most of the rust code (and some of the elixir code) in the bindings is now programmatically generated from the Elixir code as part of a separate build step.

I have added support for the most common floating point functions, including some special functions. With the new programmatic code generation tools, adding support for more functions is trivial. The functions are currently only vectorized in the first argument, but I plan on adding vectorization on the other arguments.

I have added support for drawing random variables according to the given distribution. For example, one can draw 15 values from a Cauchy distribution:

iex> SciEx.Random.draw_from_cauchy(0.0, 1.0, 15)
#SciEx.F64.Array1<[0.5075428433765116, 6.5273500168873, -0.09238828295381869,
-0.5017462511724131, 0.31446532306974234, -0.12084509091920821, -0.35170244950655,
-7.796726646392933, 0.4049382192392038, -1.5996980781784342, -1.0434302700931213, 
0.8270864602500656, -1.7333229127273764, 0.541653889771595, 129.4432667174581],
shape=[15], strides=[1], layout=CFcf (0xf), const ndim=1>

The PDF and CDF for the supported distributions is also implemented. However, they are only vectorized on the first argument.

I’m not sure on whether I’m taking the right approach regarding the rust code, though. I am writing (and compiling) separate rust functions for each supported type combination because that makes the interop with Elixir easier, although I could probably get away with generating less rust code…

I have also added support for arithmetic operators on arrays and scalars. Array types and dimensions must match, otherwise SciEx will raise a mistake. Everything is very explicit (you can use SciEx.Operators).

iex(4)> a = SciEx.Random.draw_from_normal(0.0, 1.0, 10)
#SciEx.F64.Array1<[0.980546308591242, -0.4124865578379766, -0.2678229458353892, 0.13895119159920113, -0.05895938590440134, 0.5463383504051798, -0.47139019152800626, 1.1432875805108502, -1.053279376911445, -0.4890897953025381], shape=[10], strides=[1], layout=CFcf (0xf), const ndim=1>
iex(5)> b = SciEx.Random.draw_from_normal(0.0, 1.0, 10)
#SciEx.F64.Array1<[-0.6112979257481231, 0.04187616378275234, 0.024658043411414837, -0.7423086166282216, 2.5287898587121993, -1.5968438074887625, 0.3028680311098233, 1.1657682958696634, -1.3245816353940532, 1.1693983808728805], shape=[10], strides=[1], layout=CFcf (0xf), const ndim=1>
iex(6)> 0.5 * a + 0.7 * b # This will fail beause it will try to use the Kernel operators
** (ArithmeticError) bad argument in arithmetic expression: 0.5 * #SciEx.F64.Array1<[0.980546308591242, -0.4124865578379766, -0.2678229458353892, 0.13895119159920113, -0.05895938590440134, 0.5463383504051798, -0.47139019152800626, 1.1432875805108502, -1.053279376911445, -0.4890897953025381], shape=[10], strides=[1], layout=CFcf (0xf), const ndim=1>
    :erlang.*(0.5, #SciEx.F64.Array1<[0.980546308591242, -0.4124865578379766, -0.2678229458353892, 0.13895119159920113, -0.05895938590440134, 0.5463383504051798, -0.47139019152800626, 1.1432875805108502, -1.053279376911445, -0.4890897953025381], shape=[10], strides=[1], layout=CFcf (0xf), const ndim=1>)
iex(6)> use SciEx.Operators # this imports the operators in SciEx and hides the ones in the Kernel
SciEx
iex(7)> 0.5 * a + 0.7 * b  
#SciEx.F64.Array1<[0.06236460627193485, -0.17692996427106167, -0.1166508425297042, -0.45014043584015456, 1.7406732081463385, -0.8446214900395438, -0.02368747398712684, 1.3876815973641894, -1.4538468332315597, 0.5740339689597472], shape=[10], strides=[1], layout=CFcf (0xf), const ndim=1>

Because the Rust array library I’m using (GitHub - rust-ndarray/ndarray: ndarray: an N-dimensional array with array views, multidimensional slicing, and efficient operations · GitHub) makes it so easy to parallelize computations across the array, I have added support for parallel computation. It is now automatically activated for each vectorized function (e.g. SciEx.sin(x), SciEx.exp(x), SciEx.erf(x), etc.) if the array size exceeds a given threshold, which was determined experimentally (one doesn’t want to activate parallelism for smal arrays because it brings some overhead). One can either enable choose a different cutoff for each function on a call by call case (or for the current process) or enable or disable paralellization for all sizes (again, on a call-by-call basis or for the current process).

The following benchmark shows the speed gains (this was made on a “busy” laptop with VScode and a web browser open, I plan on doing more rigurous benchmarks with everything but the benchmark turned off). The speed gain is larger for larger arrays, as is expected.

All parallelism is handled by the Rust side (using GitHub - rayon-rs/rayon: Rayon: A data parallelism library for Rust · GitHub, which is integrated as part of ndarray as part of a compile-time option). There are probabbly ways of controlling parallelism in a more fine grained way (i.e. try to “reserve” some CPU cores specifically for Rust or something like that), and I might think about it when I have a bit more time).

Operating System: Linux
CPU Information: AMD Ryzen 7 5800H with Radeon Graphics
Number of Available Cores: 16
Available memory: 6.71 GB
Elixir 1.16.0
Erlang 25.3
JIT enabled: true

Benchmark suite executing with the following configuration:
warmup: 1 s
time: 1 s
memory time: 0 ns
reduction time: 0 ns
parallel: 1
inputs: 10^6, 10^7
Estimated total run time: 8 s

Benchmarking Parallel - sin(x) with input 10^6 ...
Benchmarking Parallel - sin(x) with input 10^7 ...
Benchmarking Sequential - sin(x) with input 10^6 ...
Benchmarking Sequential - sin(x) with input 10^7 ...
Calculating statistics...
Formatting results...

##### With input 10^6 #####
Name                          ips        average  deviation         median         99th %
Parallel - sin(x)          176.49        5.67 ms    ±31.22%        5.04 ms       11.54 ms
Sequential - sin(x)         70.36       14.21 ms     ±0.93%       14.20 ms       14.58 ms

Comparison:
Parallel - sin(x)          176.49
Sequential - sin(x)         70.36 - 2.51x slower +8.55 ms

##### With input 10^7 #####
Name                          ips        average  deviation         median         99th %
Parallel - sin(x)           35.19       28.41 ms     ±8.11%       28.23 ms       35.24 ms
Sequential - sin(x)          7.11      140.62 ms     ±0.36%      140.62 ms      141.15 ms

Comparison:
Parallel - sin(x)           35.19
Sequential - sin(x)          7.11 - 4.95x slower +112.20

The main limitation of ndarray is that it seems like it currently can’t work on arrays larger than the available RAM memory. There used to be some talk about making ndarray run on top of an arrow backend (which is what polars dos (and by extension, GitHub - elixir-explorer/explorer: Series (one-dimensional) and dataframes (two-dimensional) for fast and elegant data exploration in Elixir · GitHub), but it seems like the arrow backend doesn’t play well with multidimensional arrays, which is the whole point of ndarray. While it’s not obvious what the use of 4D+ arrays is, 2D and 3D arrays have a natural real-world usage when encoding gridded data in the real word (and 2D arrays can always be implemented as matrices, which have even more real world applications).

ndarray has linear algebra extensions, which means one could write matrix-manipulating functions as part of SciEx while reusing the same backend, but since linear algebra usually operates exclusively on 1D and 2D arrays (matrices) it’s probably better to implement specific matrix type backed by GitHub - sarah-quinones/faer-rs: Linear algebra foundation for the Rust programming language · GitHub (it seems easy to convert from 2D arrays to faer’s matrix). I’ll leave that for later, as most of the things I plan on using SciEx for don’t require matrices for now. I do understand that if I want to gain traction (i.e. people to help develop the library) I’ll have to implement a bit more than what I intend to use myself.

I wonder if in the future it could be possible to integrate SciEx with Explorer dataframes without using Elixir as a bridge, as Explorer is really useful to read columnar data, which is an important format for 1D arrays (and in some cases for 2D arrays).

I have added some new improvements (already available in the GitHub repo).

Discrete cosine transforms (the classical demonstration is image compression or gaussian blur, but I haven’t yet implemented support for importing images as an array of pixels; there’s already an ndarray package that does it, I just havdn’t had time to bind it):

iex> use SciEx.Operators # activate the SciEx arithmetic operators, which are vectorized
SciEx
iex> alias SciEx.Float64.Array1 # We'll use 64-bit floats, but we could have used 32-bit floats instead
SciEx.Float64.Array1
iex> x = Array1.linspace(0.0, 10.0, 11) # create equally spaced values as an example 
#SciEx.Float64.Array1<[0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], shape=[11], strides=[1], layout=CFcf (0xf), const ndim=1>
iex> y = SciEx.FFT.dct2(x) # Apply type-II discrete cosine transform
#SciEx.Float64.Array1<[110.0, -48.87159242255254, 5.218048215738236e-15, -5.271101174621571, 4.163336342344337e-15, -1.7623006857965466, 5.218048215738236e-15, -0.7639325745301744, 2.220446049250313e-15, -0.3060225567178963, 1.7763568394002505e-15], shape=[11], strides=[1], layout=CFcf (0xf), const ndim=1>
iex> x0 = SciEx.FFT.dct3(y) / (2 * 11) # invert the DCT and correct for normalization             
#SciEx.Float64.Array1<[2.5837917664003644e-15, 0.9999999999999997, 1.9999999999999998, 3.0000000000000004, 4.000000000000001, 4.999999999999999, 5.999999999999999, 6.999999999999999, 7.999999999999999, 9.0, 10.0], shape=[11], strides=[1], layout=CFcf (0xf), const ndim=1>

I have also added support for complex numbers. With complex numbers I can bind the general complex-to-complex FFT (and not only the rea-to-real DCT):

iex> alias SciEx.Complex64, as: C # use complex numbers where the real and imaginary part are 64-bit floats
SciEx.Complex64
iex> C.new(0.0, 0.0)
#SciEx.Complex64<0 + 0i>
iex> C.new(1.0, 2.0)
#SciEx.Complex64<1 + 2i>
iex> C.new(0.0, 0.0) + C.new(1.0, 2.0) # the SciEx operators work on complex numbers
#SciEx.Complex64<1 + 2i>
iex> C.new(1.0, 2.0) / C.new(0.5, 1.0)
#SciEx.Complex64<2 + 0i>

These complex number are not actual normal Elixir structs with a real and imaginary field. They are implemented as opaque Rustler resources (even the string representation shown here is generated by Rust!), so that I can offload all operations on complex numbers to rust instead of reimplementing them badly in Elixir. These rust complex numbers play well with arrays (that is, one can have arrays of complex numbers) and implement the FFT with that.

I’m generating most of the rustler bindings programatically from Elixir using EEx templates and some custom utilities for more advanced stuff like vectorizing functions on some but not all arguments. Vectorization is interesting because it creates a combinatorial explosion on the function arguments. For example, if you have a function of 3 arguments (let’s say, the PDF for a normal distribution - normal(x, mu, sigma), you have 2^3=8 possibilities for the combinations of array and scalar arguments. This most straightforward implementation uses requires 8 different rustler bindings (i.e. rust functions), which I generate programmatically from Elixir based on the argument types. I could probabbly avoid this combinatorial explosion by being a bit smarter with the encoding adn decoding between Rust and Elixir and by doing a bit more work on the rust side, but I’m actually quite happy on how this turns out - having 8 functions instead of 1 is not such a big difference. And in any case, the Rust crate compiles in < 2min, and the extra number of functions shouldn’t cause performance problems.

If anyone wants to contribute, I believe the best way would be to make it possible to use the package with RustlerPrecompiled, because making it so that the user can use the package without a rust compiler really lowers the barrier to adoption. I’m not at all familiar with rustler precompiled.

Vix implements write_to_binary/2 which returns a binary of arbitrary image data. I’d be happy to collaborate on an efficient way for that data (plus the required metadata) to be passed to SciEx.

Due to the work by @akash-akya, we already have bi-drectional integration with evision and Nx so I’m confident this could be done at the NIF level without copying.

Thanks, but I don’t think that will be necessary. You can already do that from an extension of ndarray: ndarray_image - Rust

Supposedly (according to the docs), you feed it an image and it gives you the appropriate ndarray for free. I just haven’t had the time to bind it.

Not necessary, agreed. However:

• ndarray_image works on image paths. Sometimes working on image binary is useful - for example if you’ve already done a bunch of pre-processing on an image before numerical processing. Thats quite common.
• The number of image types is limited (to my eyes)
• Users of vix/image may well want to make use of your library and without a bridge thats not going to be practical (I consider having to write an image to an intermediate file so that ndarray_image can ready it back to be impractical - perhaps thats not fair).

Yes, you’re right, that makes sense. Feel free to contribute in any way you might find useful, thanks!

Happy to take on providing a binary and the associated metadata (dimensions, data type, band order) but I’ve not got the competence to do the work on the Rust side to convert the binary to an ndarray. I’ll work (hopefully with the support of @akash-akya) on the libvips side.

I have added support for loading 2D-arrays from grayscale PNG images and support for contouring “heightmaps” (represented as 2D arrays) using the marching-squares algorithm. @kip I’m sorry, but I haven’t been able to work on integrating SciEx with your image processing tools yet.

I plan on wrapping my own kernel density estimation (KDE) functions and expose them to Elixir as part of SciEx. Once this is done, I plan on using SciEx for the most computationally heavy parts of my Quartz plotting library.