Mojo Scientific Library (MSL) v0.1.0 release!

Community Showcase
1

Mojo Scientific Library (MSL) v0.1.0 :fire:

Hello there! I am back again :). I am sharing the first release of MSL (Mojo Scientific Library), a comprehensive collection of scientific computation routines written in pure Mojo, derived from the GNU Scientific Library (GSL). It’s been a fun project to work on porting GSL routines has been a great way to learn the math behind numerical computing (people have come up with some crazy tricks over the years!). I think it’s now reached the point where it has enough routines to be genuinely useful, so I’m happy to share it here!

MSL is designed as the low-level scalar backend for SciJo (New release coming soon) and is now available on the modular-community channel.

Install

pixi add msl

Or use pixi-backend (checkout README for details).

What’s included in v0.1.0

There’s a lot of stuff packed in here!

  • Special functions - Airy, Bessel (J/Y/I/K orders 0, 1, and integer n), Gamma, Beta, Error functions, Legendre polynomials, Digamma, Incomplete Gamma/Beta
  • Numerical integration - Non-adaptive QNG (10→21→43→87-point), fixed Gauss-Kronrod rules (QK15–QK61), general adaptive QAG, and QAGS with Wynn epsilon extrapolation
  • Numerical differentiation - Central (5-point), forward and backward (4-point) with automatic step refinement and separate truncation/rounding error estimates
  • Interpolation - Piecewise linear, natural cubic spline, and Akima spline with eval/deriv/deriv2/integral
  • Root-finding - Bisection, Brent’s method, Newton-Raphson, and secant
  • Minimization - Brent’s method and golden section search
  • ODE solvers - Classical RK4 (fixed step) and RKF45 (adaptive step with error control)
  • Probability distributions - 14 distributions with samplers and PDFs: Gaussian, Uniform, Exponential, Gamma, Beta, Chi-squared, Poisson, Student-t, Log-normal, Weibull, Binomial, Negative Binomial, Cauchy, Laplace
  • RNG - Mersenne Twister (MT19937) with an extensible RNGAlgorithm trait so you can plug in your own generator
  • Linear algebra - Vector and Matrix with math operations, BLAS Level 1/2/3 backed by mojoBLAS, and zero-copy non-owning views for interop with NuMojo NDArrays
  • Permutations - combinatorial permutation operations

Links

MSL is now part of the MojoMath organisation alongside SciJo, StatMojo, and Matmojo. As always, contributions and feedback are very welcome! If there’s a routine you need, feel free to open an issue or PR! Happy computing!

2

Very cool! That’s a whole ton of math and other useful stuff in the new version.

I’ve also written some numerical optimization routines in Mojo for some work I’m doing, including Golden Section, Brent’s method, Powell’s method, and cyclic coordinate descent (CCD). Happy to compare notes if you’re interested.

We have some SIMD-optimized multi-dimensional interpolation too, but it’s fairly specialized for sampling FFT data so I’m afraid the code itself isn’t very general. I’m not sure if any one else would be interested in something like that, but let me know. Maybe the overall approach could be generalized.

3

Thanks!

I originally intended for this library to be a port of routines from GNU scientific library. But after some back and forth between this library and my other SciJo project (This uses routines from some other libraries), I decided it’s better to include all the routines in one place for better organization and easy of access.

Therefore, I think It’ll be great if we can generalise your algorithms and add them to the library. You can just start by adding the basic algorithms for now; we can generalize them later when the library grows, there’s no rush :slight_smile:

I will also go through your algorithms in the meantime, I never came across the CCD and SIMD optimised interpolation :slight_smile:

4

Ok. The optimization stuff should be fairly general and maybe simple to copy over. Although there may be some dependencies on other math types (like Vec or Dimension) in other parts of the library? I can’t remember off-hand.

The interpolation bits aren’t general at all. But the idea itself should be general: Re-organize memory so each sample uses a single SIMD load instruction, (not strided, by making several copies of each pixel) then apply the per-pixel weights using SIMD arithmetic. It uses quite a bit more memory than the source data (cursed by dimensionality there), but the runtime performance boost is big! It’s tradeoffs all the way down.

In the case of FFT data, we’re sampling complex pixels, so there’s actually two loads I think. One each for the real and imaginary components. I’m not sure it’s possible to generalize over complex-or-real, so if you want a real-valued version of this, that’s probably another copy of the code?