Scientific Computing in Rust

Scientific Computing in Rust

Introduction #

My primary languages right now are Matlab for scientific computing and Rust for everything else. Over the years of using Rust, I’ve been slowly watching the Rust scientific computing ecosystem evolve and in this post I play with see how Rust might take Matlab’s place in my daily workflow. Yes, I am aware there are many other scientific computing languages: Octave, Python, R, and Julia. Except for R, I have used them all many times and expect to continue to run into situations that require them. The point here isn’t about finding a replacement for Matlab as much as is about exploring the capabilities of a language that I love to execute tasks that I love.

A Few Points About My Workflow #

Since I want to use Rust in my workflow, a little discussion about what I need for scientific computing is in order.

  1. Script-based: Usually, any one set of computations fits nicely in one file which makes the story readable and self-documenting.
  2. Portable and Reproducible: For all the hate it gets, Matlab is extremely stable and with few exceptions, I can hand anyone my Matlab script and they will be able to run it.
  3. Code Like Math: Matlab code reads like math which is one of it’s strongest features and safety mechanisms.
  4. Visualization: Interactive plots are an integral part of my workflow.

Crustacean Interpretation of My Workflow #

With the important points of my workflow in mind, here’s my current thoughts on adapting the Rust ecosystem to my pedantic needs.

  1. Script-based: rust-script tool that integrates Cargo’s package management with compilation caching to create a script-like user experience. It is excellent and as of this writing, it is being integrated into Cargo itself.
  2. Portable and Reproducible: Cargo makes any rust-script just as portable and reproducible than any Matlab script.
  3. Code Like Math: is one of the weak points of using Rust. While Rust’s construction is quite mathematical since it draws inspiration from functional languages, the syntax isn’t particularly similar to reading math. There are some very nice crates that add good math tools and increase the math-like syntax. Below are a few of my favorites (thank you to any of the contributors to these awesome projects):
    1. ndarray - $n$-dim arrays
    2. ndarray-linalg - linear algebra on ndarrays
    3. ndarray-stats - statistics on ndarrays
    4. ndarray-rand - random numbers on ndarrays
    5. peroxide - a collection of scientific computing tools
    6. pandas - high-performance DataFrame library
    7. bitvec - bit vectors
  4. Visualization: Plotting has been one of the roughest points in the Rust scientific computing area. Most Rust plotting libraries are either bindings around other non-Rust libraries or static image generation with no interactivity. However, more recently egui-plot has matured into a very nice native plotting solution. Egui-plot is also under active development and we can expect some more improvements. One additional benefit to using egui-plot is that it can create WASM objects of the plots that could be embedded in a website like this; it takes a bit of effort and will be the subject of a later post.

An Example #

To see how this works together, let’s integrate and plot the Lorenz 63 system given by

$$ \begin{split} \frac{dx}{dt} &= 10(y-x) \\ \frac{dy}{dt} &= x(28-z) -y\\ \frac{dz}{dt} &= xy -\frac{8}{3}z. \end{split} $$

We use peroxide for the ODE solver and egui-plot for the interactive plotting.

Code #

#!/usr/bin/env rust-script
//! This is a regular crate doc comment, but it also contains a partial
//! Cargo manifest.  Note the use of a *fenced* code block, and the
//! `cargo` "language".
//!
//! ```cargo
//! [package]
//! rust-version = 1.91 
//!
//! [dependencies]
//! peroxide = {version ="0.39.11", features = ["rayon"]}
//! egui_plot = {version = "0.33.0" }
//! eframe = {version = "0.32.1"}
//! env_logger = { version = "0.11.8", default-features = false, features = [
//!     "auto-color",
//!     "humantime",
//! ] }
//! ```


use peroxide::fuga::*;
use eframe::egui;
use egui_plot::{Legend, Line, Plot, PlotPoint, PlotPoints};


fn main() -> Result<(), Box<dyn Error>> {
    env_logger::init(); // Log to stderr (if you run with `RUST_LOG=debug`).
    
    let initial_conditions = vec![10f64, 1f64, 1f64];
    let rkf45 = RKF45::new(1e-4, 0.9, 1e-6, 1e-3, 100);
    let basic_ode_solver = BasicODESolver::new(rkf45);
    let (_, y_vec) = basic_ode_solver.solve(
        &Lorenz,
        (0f64, 100f64),
        1e-3,
        &initial_conditions,
    )?; // Error handling with `?` - can check constraint violation and etc.
    let y_mat = py_matrix(y_vec);
    let y0 = y_mat.col(0);
    let y2 = y_mat.col(2);

    let options = eframe::NativeOptions {
        viewport: egui::ViewportBuilder::default().with_inner_size([750.0, 500.0]),
        ..Default::default()
    };

    let points: Vec<PlotPoint> = y0.iter().zip(y2.iter()).map(|(x,z)| {
        PlotPoint::new(*x,*z)
    }).collect();

    println!("The Number of Points = {:?}", y0.len());

    Ok(eframe::run_native(
        "Lorenz System",
        options,
        Box::new(|_cc| Ok(Box::new(MyApp { points }))),
    )?)
   
}

struct Lorenz;

impl ODEProblem for Lorenz {
    fn rhs(&self, _: f64, y: &[f64], dy: &mut [f64]) -> anyhow::Result<()> {
        dy[0] = 10f64 * (y[1] - y[0]);
        dy[1] = 28f64 * y[0] - y[1] - y[0] * y[2];
        dy[2] = -8f64 / 3f64 * y[2] + y[0] * y[1];
        Ok(())
    }
}

// The plotting egui App (= Matlab Figure)
#[derive(Default)]
struct MyApp {
    points: Vec<PlotPoint>,
}

impl eframe::App for MyApp {
    fn update(&mut self, ctx: &egui::Context, _frame: &mut eframe::Frame) {
        egui::CentralPanel::default().show(ctx, |ui| {
	        // Add a plot to the App/Figure
            Plot::new("Lorenz Plot")
                .legend(Legend::default())
                .x_axis_label("x(t)")
                .y_axis_label("z(t)")
                .show(ui, |plot_ui| {
                    plot_ui
                        .line(Line::new("curve", PlotPoints::Borrowed(&self.points)).name("Lorenz"));
                });
        });
    }
}

Run the Example #

  • Install rust-script instructions.
  • Then run rust-script ./the_file_name.rs. This will take a little be to execute the first time because it will need to download and install
  • Explore the interactive plot.

Conclusion #

Using Rust for scientific computing is quite doable. It particularly isn’t pretty yet (there is still an awful lot of boilerplate), but I think that will come with time and development. There are some amazing and super skilled developers volunteering their time to make these awesome tools. If you have found any tools or thoughts about scientific computing in Rust, let me know! I always looking for improvements.