Rust, For Science!

@adamisntdead

Rustfest Rome

\sum_{n=1}^{\infty} n= -\frac{1}{12}
n=1n=112\sum_{n=1}^{\infty} n= -\frac{1}{12}
\frac{\partial L}{\partial f} - \frac{d}{dt} \frac{\partial L}{\partial f'} = 0
LfddtLf=0\frac{\partial L}{\partial f} - \frac{d}{dt} \frac{\partial L}{\partial f'} = 0
\psi (x) = -i \hbar \frac{\partial}{\partial x}(-i \hbar \frac{\partial \psi}{\partial x})
ψ(x)=ix(iψx)\psi (x) = -i \hbar \frac{\partial}{\partial x}(-i \hbar \frac{\partial \psi}{\partial x})
P(A \lvert B ) = \frac{P(B \lvert A) P(A)}{P(B)}
P(AB)=P(BA)P(A)P(B)P(A \lvert B ) = \frac{P(B \lvert A) P(A)}{P(B)}
\sigma = \sqrt{\frac{\sum (x - \bar{x})^2}{n - 1}}
σ=(xxˉ)2n1\sigma = \sqrt{\frac{\sum (x - \bar{x})^2}{n - 1}}
\nabla \times G = J + \frac{\partial D}{\partial t}
×G=J+Dt\nabla \times G = J + \frac{\partial D}{\partial t}
\rho (\frac{\partial V}{\partial t} + V\cdot \nabla V) = \nabla P + \rho g + \mu \nabla^2 V
ρ(Vt+VV)=P+ρg+μ2V\rho (\frac{\partial V}{\partial t} + V\cdot \nabla V) = \nabla P + \rho g + \mu \nabla^2 V
\hat{f}(\zeta) = \int_{-\infty}^{\infty} f(x) e^{-e \pi i x \zeta} dx
f^(ζ)=f(x)eeπixζdx\hat{f}(\zeta) = \int_{-\infty}^{\infty} f(x) e^{-e \pi i x \zeta} dx
F_g = \frac{G m_1 m_2}{r^2}
Fg=Gm1m2r2F_g = \frac{G m_1 m_2}{r^2}
f(x) = \frac{1}{1 + e^{-x}}
f(x)=11+exf(x) = \frac{1}{1 + e^{-x}}
J(\theta_j) \leftarrow 0_j \pm \eta \cdot \frac{\sum (\hat{y} - y)^2 \cdot x}{2n}
J(θj)0j±η(y^y)2x2nJ(\theta_j) \leftarrow 0_j \pm \eta \cdot \frac{\sum (\hat{y} - y)^2 \cdot x}{2n}
\frac{d}{{dx}}f\left( x \right) = \mathop {\lim }\limits_{\Delta \to 0} \frac{{f\left( {x + \Delta } \right) - f\left( x \right)}}{\Delta }
ddxf(x)=limΔ0f(x+Δ)f(x)Δ\frac{d}{{dx}}f\left( x \right) = \mathop {\lim }\limits_{\Delta \to 0} \frac{{f\left( {x + \Delta } \right) - f\left( x \right)}}{\Delta }

@adamisntdead

github.com/adamisntdead

Adam Kelly

What is Scientific Computing?

High-Performance Computing

Data Science

Numerical Analysis & Simulation

Scientific Software Development Process

The Problems With Rewriting

  • Harder to debug
  • Hard to learn
  • Easy to write bad code

How Does Rust Help?

Easy Unit Tests

fn factorial(n: u64) -> u64 {
    match n {
        0 => 1,
        _ => n * factorial(n - 1)
    }
}
#[test]
fn test_factorial() {
    assert_eq!(factorial(2), 2);
    assert_eq!(factorial(3), 6);
    assert_eq!(factorial(4), 24);
}

Writing Better Code

let numbers = [1,2,5,3];

for i in 0..3 {
    println!("{} => {}", i, numbers[i]);
}
help: consider using an iterator
  |
4 |     for (i, <item>) in numbers.iter().enumerate() {
  |         ^^^^^^^^^^^    ^^^^^^^^^^^^^^^^^^^^^^^^^^

Rust's Features

  • Memory and thread safety
  • Built in benchmarking, dependency managment, documentation
  • Libraries such as Rayon and Tokio
  • Easy Integration with existing code
  • WebAssembly Integration

Easier Multithreading

Example: Numerical Integration

fn integrate<F>(f: F, a: f64, b: f64, n: u32) -> f64
where
    F: Fn(f64) -> f64,
{
    let h = (b - a) / n as f64;
    let mut result = 0.5 * f(a) + 0.5 * f(b);

    for i in 1..n {
        result += f(a + i as f64 * h);
    }
    result * h
}

Example: Numerical Integration

extern crate rayon;

use rayon::prelude::*;

Example: Numerical Integration

fn integrate_threaded<F>(f: F, a: f64, b: f64, n: u32) -> f64
where
    F: Fn(f64) -> f64 + Sync,
{
    let h = (b - a) / n as f64;
    let result: f64 = (1..n)
        .into_par_iter()
        .map(|i| f(a + i as f64 * h))
        .sum();

    h * result
}

Other Libraries

  • Built in threading libraries
  • Integration with C libraries (OpenCL, Arrayfire, BLAS)
  • Crossbeam

Integration With Other Tools

Example: MPI

#include <mpi.h>
#include <stdio.h>

int main(int argc, char** argv) {
  // Initialize the MPI environment
  MPI_Init(NULL, NULL);

  // Get the number of processes
  int world_size;
  MPI_Comm_size(MPI_COMM_WORLD, &world_size);

  // Get the processes rank
  int world_rank;
  MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

  // Get the name of the processor
  char processor_name[MPI_MAX_PROCESSOR_NAME];
  int name_len;
  MPI_Get_processor_name(processor_name, &name_len);

  // Print off a hello world message
  printf("Processor %s, rank %d\n", processor_name, world_rank);

  // Finalize the MPI environment
  MPI_Finalize();  
}

Example: MPI

[dependencies]
mpi = "0.5"

Example: MPI

extern crate mpi;

use mpi::traits::*;

fn main() {
    // Initialize the MPI environment
    let universe = mpi::initialize().unwrap();
    let world = universe.world();

    // Get the number of processes
    let size = world.size();

    // Get the processes rank
    let rank = world.rank();

    // Get the processor name
    let processor_name = mpi::environment::processor_name().unwrap();

    // Print off a hello world message
    println!("Processor {}, rank {}", processor_name, rank);
}

Tools

  • Rust Bindgen
  • Integration with other languages (PyO3, neon, ...)
  • Many prexisting bindings crates
  • WebAssembly

Interesting Developments

What's Going On Here?

  • Simulating a Double Pendulum
    • Just plain Rust
    • Uses a struct for state
  • Displaying The Simulation
    • Using GGEZ
    • Takes values from the struct

Separate

Enter WebAssembly

pub struct DoublePendulumLagrangian {
    /// Gravitational Constant.
    pub g: f64,
    /// Mass of the first bob.
    pub m1: f64,
    /// Mass of the second bob.
    pub m2: f64,
    /// Initial angle of the first bob.
    pub t1: f64,
    /// Initial angle of the second bob.
    pub t2: f64,
    /// Angular velocity of the first bob.
    pub dt1: f64,
    /// Angular velocity of the second bob.
    pub dt2: f64,
    /// Length of the rod for the first bob.
    pub l1: f64,
    /// Length of the rod for the second bob.
    pub l2: f64,
}

impl DoublePendulumLagrangian {
    // ...

Enter WebAssembly

const { DoublePendulumLagrangian } = wasm_bindgen;

wasm_bindgen('../out/main_bg.wasm')
    .then(runApp)
    .catch(console.error)

function runApp() {
  const elem = document.getElementById('container')
  const two = new Two({ fullscreen: true })
    .appendTo(elem)

  const g = 9.8
  const m1 = 2.0
  const m2 = 2.0
  const t1 = 2.0
  const t2 = 1.5
  const l1 = 100
  const l2 = 100

  const pendulum = DoublePendulumLagrangian.new(
    g, m1, m2, 
    t1, t2, l1, l2
  )
  drawBobs(two, pendulum)
  // ...

Jupyter Notebooks

What about in Rust?

EXCVR

EXCVR

EXCVR

  • Evaluation Context
  • Jupyter Kernel
  • REPL

So What's The Problem?

Thank you!

@adamisntdead

github.com/adamisntdead

Adam Kelly

Copy of Rust, For Science!

By Adam Kelly

Copy of Rust, For Science!

  • 51
Loading comments...

More from Adam Kelly