An introduction to Rust programming

Introduction

The purpose of this course is to teach you how to create a Rust program from scratch. You will learn about the basics of Rust programming the Rust compiler and the Rust class system. This will be hard if your only programming experience comes from R and Python ;-)

First of all there is no interactive mode in Rust like there is in R and Python. Rust is no scripting language - you need to compile your program before you can use it.

Practice makes perfect

Programming is like any skill, the more you practice the better you get. It’s really important that you keep using what you have learned after the course is completed otherwise there is a good chance you will forget everything and you’ll be back to square one.

Why Use Rust?

• Rust is a compiled language and very fast.
• Rust is very strict with access rights to variables and this adds to security
• The Rust compiler is EXTREMELY helpful and therefore coding is easier as in e.g. C++

This said coding in Rust is still more complicated than in R or Python. Therefore you should probably only use Rust if a process you frequently use takes a lot of time. Rust in general is easily 10 times faster than R or Python. You can compare Rust with C and C++, but programming in Rust is WAY easier than programming in C or C++. If you want to read more about that I recommend this source: https://combine-lab.github.io/blog/.

A short boil down: Rust is cool - mainly because the Rust compiler/package manager (cargo) is a work of art, the programs are really fast and the memory management makes up to be a very secure language.

Which projects use Rust?

A very short list:

• CellRanger
• alevin-fry

Even this extremely short list points to the main benefit of Rust - speed.

How will this course work?

Rust is OS independent - meaning you all will install it on your computer. Please follow these steps to install it on your system: https://www.rust-lang.org/tools/install.

This document contains code sniplets that will help you getting started, but you will do the coding in a programming IDE. I am working with Sublime Text, but you are welcome to use whichever is to your liking. We will use the terminal to compile and run our programs. I recommend you to use git to store your code. Please create a github or gitlab account if you not already have one.

Introduction

Programming in Rust is totally different from coding in C or C++ - and of cause also very different from coding in R or Python. I will talk about typed variables later on, but want to first highlight the security features Rust implements.

In Rust a variable is ‘owned’ by a function - possibly comparable to R. But instead of R which simply copies the data into every function in Rust the variables can not be changed inside of the function. It is especially prohibited to e.g. give a function a file object to write to. Only the function that did create the file object can write to it.

This is also true for e.g a HashMap that stores objects (some data) in combination with e.g. a string key. If you want to modify the stored object you need the class that ‘owns’ this object make the changes. I have not hidden any of these problems inside this tutorial. That is something you need to learn the hard way.

Variables

Rust uses strong typed variables. This does mean that you need to define your variables as a specific type. As an example you can store integer values as u8, u16, u32, u64, u128, i8, i16, i32, i64, i128, usize or integer. All of these types can not be compared against each other. So if you want to create some data structures it does help to choose the same type for all data (eg. u64). The difference between the integer types are the range of values they can contain: u8 and i8 are 8 bit long whereas u128 and i128 are 128 bit long; the u types are unsigned and can not contain negative values. Usize is the default integer type and is dependent on your computer type. On 64 bit computers it is the same as a u64 - but you still can not compare it directly to u64.

Target Audience

This tutorial does not explain to you how to program. I assume you are absolutely able to do that. In addition there a re tons of other tutorials that explain that to you.

What I found most important for my work is to try to understand code I have not written and hunt bugs I am not responsible for. In addition I assume you all benefit the most of a larger set of code snippets you can go back to and copy code from. Code you can double check how something could work.

Examples

Here I will show you the minimal Rust programming:

• define a variable
• do a small calculation

Integers

Using RStudio at the moment I can show you how you define a variable in Rust:

let x = 5;
let y = 7;
let z = x*y;

rprintln!("{x} * {y} = {z}");

## 5 * 7 = 35

The print statement here is not Rust but specific for the R extension rextendr. The Rust equivalent to the rprintln!() is println!(). Please ignore my usage of rprintln! - it is an artifact from me using Rstudio to compile this tutorial.

The interesting part here is how you define variables: “let <varname>:<type> = <value>;”

Why did we not define the type of integer we want to store?

Here Rust silently assumes you want to use the default integer type usize. If you want to specify the integer type you want to use here you could write the following:

let x:u8 = 5;
let y:u8 = 7;
let z = x*y;

rprintln!("{x} * {y} = {z}");

## 5 * 7 = 35

OK - this works in RStudio - what about breaking this code?

Create a new Rust program

Rust is a compiled language and therefore interacting with it using a jupyter notebook or Rmd file does limit the language a lot. Hence it makes more sense to compile Rust programs and use them on the command line.

As an example of why - I deliberately broke the integer example:

let x:u8 = 5;
let y:u16 = 7;
let z = x*y;

rprintln!("{x} * {y} = {z}");

Compiling this in RStudio breaks the html generation and I can not show the error here using this code. But if I create a stand alone ‘program’ of it it works:

First - create an empty class structure:

cargo new workshop --bin

Open the created file “workshop/srv/main.rs” and replace the “println!(”Hello, world!“);” with the broken code we had before.

Compile the new program (normally you state cargo build -r I had to add here as the jupyter notebook could not handle the terminal color commands):

cd workshop
cargo build -r 2>&1 | sed -r "s/\x1B\[[0-9;]*[mGKH]//g"

##    Compiling workshop v0.1.0 (/home/med-sal/git_Projects/Rust_Programming_1/workshop)
## error[E0308]: mismatched types
##  --> src/main.rs:4:13
##   |
## 4 |   let z = x*y;
##   |             ^ expected `u8`, found `u16`
## 
## error[E0277]: cannot multiply `u8` by `u16`
##  --> src/main.rs:4:12
##   |
## 4 |   let z = x*y;
##   |            ^ no implementation for `u8 * u16`
##   |
##   = help: the trait `Mul<u16>` is not implemented for `u8`
##   = help: the following other types implement trait `Mul<Rhs>`:
##             <u8 as Mul>
##             <u8 as Mul<&u8>>
##             <&'a u8 as Mul<u8>>
##             <&u8 as Mul<&u8>>
## 
## Some errors have detailed explanations: E0277, E0308.
## For more information about an error, try `rustc --explain E0277`.
## error: could not compile `workshop` (bin "workshop") due to 2 previous errors

Here you see the output from the Rust compiler. Do you see the beauty of this? The compiler explains the error. The most reasonable fix here would of cause be to not define the y as u16 in the first place - right?

This becomes even more beautiful if you want to know more about the error: you can simple ask Rust to explain more:

rustc --explain E0308

## Expected type did not match the received type.
## 
## Erroneous code examples:
## 
## ```
## fn plus_one(x: i32) -> i32 {
##     x + 1
## }
## 
## plus_one("Not a number");
## //       ^^^^^^^^^^^^^^ expected `i32`, found `&str`
## 
## if "Not a bool" {
## // ^^^^^^^^^^^^ expected `bool`, found `&str`
## }
## 
## let x: f32 = "Not a float";
## //     ---   ^^^^^^^^^^^^^ expected `f32`, found `&str`
## //     |
## //     expected due to this
## ```
## 
## This error occurs when an expression was used in a place where the compiler
## expected an expression of a different type. It can occur in several cases, the
## most common being when calling a function and passing an argument which has a
## different type than the matching type in the function declaration.

More complicated error hunting:

The bugs so far were not very complicated - right? Now we need to dig deeper and I recommend you to use your programming IDE of choice to look into a new broken program. Rstudio could possibly also do.

Lets look in the simple_problem folder in this github repo.

Download (clone) the repo using git git clone https://github.com/stela2502/Rust_Programming_1. Open the folder ‘simple_error’ in your IDE and open a terminal and cd into the cloned folder.

There you run:

cd simple_error
cargo build -r  2>&1 | sed -r "s/\x1B\[[0-9;]*[mGKH]//g"

##    Compiling simple_error v0.1.0 (/home/med-sal/git_Projects/Rust_Programming_1/simple_error)
## error: expected `;`, found `println`
##  --> src/main.rs:6:31
##   |
## 6 |     let res = sum( data, ids )
##   |                               ^ help: add `;` here
## 7 |
## 8 |     println!("Summing up these ids ({:?}) of the data ({:?}) = {res}", data, ids );
##   |     ------- unexpected token
## 
## error[E0425]: cannot find value `data` in this scope
##   --> src/main.rs:15:16
##    |
## 15 |         sum += data[i];
##    |                ^^^^ not found in this scope
## 
## For more information about this error, try `rustc --explain E0425`.
## error: could not compile `simple_error` (bin "simple_error") due to 2 previous errors

It is almost no fun hiding programming errors as the compiler will already tell you what to do - right? I think this nevertheless will be fun…

Problem 1

Each line needs to end with a ‘;’. All but the last line in a function if you want to return a value.

Promblem 2

We do not have a ‘data’ variable in the sum function. Here the data is called ‘v’ - fix that.

Now the code compiles:

cd simple_fixed
cargo build -r 2>&1 | sed -r "s/\x1B\[([0-9]{1,2}(;[0-9]{1,2})?)?[m|K]//g"
target/release/simple_error

##     Finished release [optimized] target(s) in 0.00s
## Summing up these ids ([1, 2, 3, 4, 5, 6, 7, 8, 9]) of the data ([5, 6, 7, 8]) = 10

Does the output make sense?

Problem 3

No it is wrong: We need to switch the variables in the println! call.

Problem 4

Wait - we tried to sum up the values 6, 7, 8 and 9 - that is not 10! How do you fix that?

This one is rather complicated and I recommend you to add a test case before you try to fix it.

Tests in Rust

Problem 4 is a good example for a really annoying bug. To get them and also make sure you will not re-introduce this bug later it makes a lot of sense to write a test that detects the error you have in your code before you try to fix it:

Rust has a nice way to do so: Add this to the script

#[cfg(test)]
mod tests {
    use crate::sum;
    #[test]
    fn check_sum() {
      
      let data:Vec<usize> = vec![1,2,3,4,5,6,7,8,9];
      let ids:Vec<usize> = vec![5,6,7,8];
      
        assert_eq!( sum( &data, &ids ), 30 );
    }
}

And now you can test your script by

cd simple_fixed
cargo test -r 2>&1 | sed -r "s/\x1B\[([0-9]{1,2}(;[0-9]{1,2})?)?[m|K]//g"

##     Finished release [optimized] target(s) in 0.00s
##      Running unittests src/main.rs (target/release/deps/simple_error-b85c7731fc253ab5)
## 
## running 1 test
## test tests::check_sum ... FAILED
## 
## failures:
## 
## ---- tests::check_sum stdout ----
## thread 'tests::check_sum' panicked at src/main.rs:31:9:
## assertion `left == right` failed
##   left: 10
##  right: 30
## note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
## 
## 
## failures:
##     tests::check_sum
## 
## test result: FAILED. 0 passed; 1 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
## 
## error: test failed, to rerun pass `--bin simple_error`

I hope you can find and fix the bug.

A short break

How do you get data into Rust?

Getting data into a compiled langue has always been a problem for me. Hence I’ll tell you how I finally managed to get data into Rust - just copy paste the solution - I do the same.

The git repo contains a folder ‘read_data’ which is a simple Rust package. Looking into it you see a Cargo.toml. This file describes your package and is comparable to the R DESCRIPTION file. Here you define the package and the dependencies. Just look into the example. It should be self explanatory.

The src folder contains the code files with main.rs defining the executable and the lib.rs defining the library part of the package. The additional data.rs file contains the code for the library in this package and needs to be mentioned in the lib.rs file to be compiled.

Let’s dig through the code.

Start with the main.rs file.

This file will be compiled to be the executable. The program starts with the main function.

• lines 0..7 import the used classes
• lines 13..22 define the command line options
• lines 24..55 define the main function

In the main function most of the code is a timer functionality that I do think is extremely useful.

let now = SystemTime::now();

// do the stuff

match now.elapsed() {
    Ok(elapsed) => {
        let mut milli = elapsed.as_millis();
        
        let mil = milli % 1000;
        milli= (milli - mil) /1000;

        let sec = milli % 60;
        milli= (milli -sec) /60;

        let min = milli % 60;
        milli= (milli -min) /60;

        println!("finished in {milli} h {min} min {sec} sec {mil} milli sec");
    },
    Err(e) => {println!("Error: {e:?}");}
}

Using ‘let <var_name>: <type> = <value>’ defines variables in Rust. The type does not need to be included in the definition if the compiler can identify the type itself.

The ‘match’ part is the main Rust error handling functionality. It is comparable to the C switch statement as it allows to match a value against multiple patterns and execute code based on the first matching pattern.

The function ‘now.elapsed()’ can fail and therefore returns a “Result<Duration, SystemTimeError>” instead of only a Duration. This allows for the match to check if an error has occurred using the two functions Ok and Err that both take a Result as input:

match <statement> {
    Ok(<varname>) => { <do something with the variable> },
    Err(e) => { <do something with the error> },
};

I am sure you will like this ;-)

The reading of the command line options:

let opts: Opts = Opts::parse();
let mut sep = '\t';
if &opts.sep != "\\t"{
    println!("I set sep to {}", opts.sep );
    sep = opts.sep.chars().next().unwrap(); 
}

The first line parses the command line. The rest is my way to handle the problem that you can not give Rust a ‘\t’ on the command line as the used char class can only contain one literal. So I have used the tab as default instead of dealing with this problem in any other way.

The reading and printing of the data table

let data = Data::read_file( &opts.data , sep );
data.print();

Quite a lot of overhead to just read a tab separated file into a script - or? But you have not even seen most of it as what you see here first defines a data variable using the Data::read_file() function and then uses the Data::print() function to print the contents of the data object. Therefore the real code loading the data is hidden in the library.

The lib.rs file

pub mod data;

Oh that is boring. It only imports the code from data.rs …

The data.rs file

This file now defines the Data class.

• lines 0..5 define the dependencies
• lines 8..14 define the data part of the class
• lines 17..110 code the functions
• lines 112..126 implement the tests for this class

A class in Rust is defined by creating the data structure for the class (struct Data {}) and then implementing functions for the class (impl Data {}). Functions can be defined with

<accessability> fn <function name> (<variables>) -> <return value> { 
    <code> 
} 

Lets look into the new function:

pub fn new( rows:usize, cols:usize, data: Vec::<f64>, rownames: Vec::<String> ) -> Self {

    let ret = Array::from_iter(&mut data.into_iter());
    let data = ret.into_shape([rows, cols]).unwrap();

    Self {
        rows, 
        cols,
        rownames,
        data,
    }
}

First this function can be called from outside the class - it is public. And this is what it does: it takes ownership of the data using the into_iter() function. This means that we can no longer access data in the main function. This is part of Rusts more secure memory management: Only one function can have access to the data. We then create an Array object from this data (from_iter) and finally we re-format the data into a two dimensional table (into_shape). The ending unwrap() call handles any possible error for the reshape process.

The Self{} statement at the end creates the object from it’s parts. Here we have the rows, cols, the rownames (from the function call) and the data array. This line is not ending with a ‘;’. Here the compiler instead returns the Self object.

This is quite simple - or?

The read_file function

The read_file function shows a general problem of compiled languages. Normally you should know how much data you want to store in the objects as the memory needs to be requested per variable. Adding more ‘space’ to a variable is quite costly. Hence we iterate over the data twice. Lines 32..52 read the data and ‘only’ collect the column and row counts.

The ‘line’ object contains one line of the data file and the reader.lines() function returns a ‘Result’ - meaning it is normal that this function might not work. And - again - this is the common way to treat a Result in Rust:

let mut header =true;
for line in reader.lines() {
    if header{
        // just drop the header
        header = false;
        continue;
    }
    match line {
        Ok(line) => {
            header =true;
            for mut val in line.split( sep ).collect::<Vec<&str>>(){ ## here is the mut!
                if header{
                    names.push( val.to_string() );
                    header = false;
                }else {
                    val = val.trim();
                    let v = match val.parse::<f64>() {
                        Ok( v ) => v,
                        Err(_err) => {
                            match val.parse::<usize>(){
                                Ok(v) =>  { 
                                    v as f64
                                },
                                Err(err) => {
                                    eprintln!("I could not parse '{val}' to usize or f64 {err:?}");
                                    0.0
                                },
                            }
                        },
                    };
                    arr.push( v );
                }
            }
        },
        Err(err) => {
            panic!("Unexpected error reading the csv file: {err:?}");
        }
    }
}

What does this ‘mut’ mean? Rust differentiates between variables you only store a value in and variables that you can modify. So whenever you want to assign a value to a variable more than once you need to declare the variable as mut.

for mut val in line.split( sep ).collect::<Vec<&str>>()

Translates to: split the line by ‘sep’ and give me ownership of the variables you create as Vec<&str>; iterate over the variables in the vector and I want to be able to modify these values.

Where do we modify the value? val = val.trim(); does remove whitespace from the strings. This is necessary as ” 2” can not be parsed as 2:<usize>. Tried that and it broke ;-) Rust also does not want to parse an int as float and vice versa. So therefore I need to differentiate between the two options here, too. The only other interesting part here is the let v = match val.parse::<f64>() {. This tries to convert the data into a f64 and if that fails it converts it to a usize. And if both ties fail it throws an error, but does not break - so more like a warning. In this case it returns 0.0 (a f64).

And now try the tool:

cd read_data
cargo build -r 2>&1 | sed -r "s/\x1B\[([0-9]{1,2}(;[0-9]{1,2})?)?[m|K]//g"
target/release/read_data -h

##     Finished release [optimized] target(s) in 0.02s
## read_data 1.0.0
## Stefan L. <stefan.lang@med.lu.se>
## Example of how to load numeric data from a tab separated text file into Rust
## 
## USAGE:
##     read_data [OPTIONS]
## 
## OPTIONS:
##     -d, --data <DATA>    the data (text file) [default: testData/Spellman_Yeast_Cell_Cycle.tsv]
##     -h, --help           Print help information
##     -s, --sep <SEP>      the column separator for the file [default: \t]
##     -V, --version        Print version information

And we can even run the test with some data file:

cd read_data
target/release/read_data -d testData/Spellman_Yeast_Cell_Cycle.tsv

## I got 256 rows and 16 cols in this data
## [[-0.35, 0.1, 0.47, 0.39, 0.17, ..., 0.5, 0.05, -0.2, -0.11, 0.3],
##  [-2.06, -1.98, -1.64, -0.21, 1.21, ..., -0.02, 0.88, 1.3, 1.29, 1.29],
##  [-1.61, -2.55, -1.09, -0.55, 0.98, ..., -0.16, 0.77, 1.27, 1.17, 1.3],
##  [-0.8, -1.25, -0.06, 0.27, 0.51, ..., 0.43, 0.25, 0.21, 0.43, -0.04],
##  [-0.39, -0.58, 1.08, 1.21, 0.52, ..., 0.34, 0.39, 0.1, -0.47, -0.05],
##  ...,
##  [-0.41, -0.57, 1.15, 0.69, 0.61, ..., 0.64, 0.42, -0.09, -0.03, -0.5],
##  [-0.79, -0.25, 1.18, 0.84, 0.19, ..., 0.8, 0.08, -0.24, -0.14, -0.27],
##  [-0.68, -0.19, 1.23, 0.84, 0.14, ..., 0.98, 0.07, -0.24, -0.35, -0.22],
##  [-0.48, -0.42, 0.87, 0.92, 0.67, ..., 0.8, 0.45, 0.16, -0.09, -0.35],
##  [-0.8, -1.22, -0.77, -1.11, 0.02, ..., -0.24, -0.11, 0.29, 0.58, 0.81]]
## finished in 0 h 0 min 0 sec 0 milli sec

Now we know how to load a numeric table into Rust. Let’s use that!

Use the Data class to implement a kmeans algorithm

Now you know how to create a package and even know how to load numeric tables into Rust. But how do we interact with this table and how does a kmeans algorithm work?

I could of cause explain that in detail here, but that would not help you a lot - or? Instead I show you how I would find this info nowadays: using ChatGPT. So why not try it here, too.

This was my conversation with ChatGPT:

Hi Chatgpt. Can you tell me how to code kmeans in Rust if the data is stored in 2d ndarrays  and the rows should be clustered.

Sure, I can help you with that! Here's an implementation of the KMeans algorithm in Rust using the ndarray crate to handle the 2D data:

rust

use ndarray::{Array2, Axis};

fn kmeans(data: &Array2<f64>, k: usize) -> Array2<f64> {
    let n = data.nrows();
    let d = data.ncols();
    let mut centroids = data.slice(s![..k, ..]).to_owned();
    let mut labels = Array2::zeros((n, 1));

    loop {
        let mut new_labels = Array2::zeros((n, 1));
         for i in 0..n {
            let mut min_distance = f64::INFINITY;
            for j in 0..k {
                let distance = (data.slice(s![i, ..]) - centroids.slice(s![j, ..])).norm();
                if distance < min_distance {
                    min_distance = distance;
                    new_labels[[i, 0]] = j as f64;
                }
            }
        }
        if new_labels == labels {
            break;
        }
        labels.assign(&new_labels);
        for j in 0..k {
            let mask = labels.iter().enumerate().filter(|(_, &l)| l == j as f64)
                .map(|(i, _)| i).collect::<Vec<_>>();
            let cluster_data = data.select(Axis(0), &mask);
            if !cluster_data.is_empty() {
                centroids.slice_mut(s![j, ..]).assign(&cluster_data.mean_axis(Axis(0)).unwrap());
            }
        }
    }
    labels
}

This implementation takes as input an Array2<f64> representing the data, and an integer k representing the number of clusters to find. It returns an Array2<f64> of the same shape as the input data, containing the cluster labels for each row.

The algorithm works as follows:

    Initialize centroids to k randomly chosen rows from the input data.
    Loop until convergence:
    a. Assign each row to the cluster with the closest centroid.
    b. Update the centroid of each cluster to be the mean of its assigned rows.
    c. If the cluster assignments haven't changed, exit the loop.
    Return the cluster labels for each row.

The implementation should be fairly self-explanatory, but I'll provide some additional details:

    We use the norm() method to compute the Euclidean distance between a row of data and a centroid.
    The slice() method is used to extract subsets of the data and centroids arrays.
    The select() method is used to extract the rows of data belonging to a given cluster.
    The mean_axis() method is used to compute the mean of the rows belonging to a given cluster.

Hope that helps! Let me know if you have any questions.

Implement the kmeans algorithm using the AI output?

Sure - why not. Let’s recapitulate what the kmers algorithm should do and check if we find the respective lines in the code.

The data should already be in memory and it should be a ndarray. Here it is obviously called “data”. But I would not like to have an Array2<f64> back. This function should rather return the cluster identities for each row as a Vec<usize>. Keep that in mind for later.

The next step in the kmeans algorithm is to select a random set of initial ‘centromers’. This step is totally missing in the AI function. To select a random set of starting rows we need a random number generator. In Rust this works like that:

use rand::Rng;
let mut rng = rand::thread_rng();
rng.gen_range(0..100); // gives a random number in this range.

Next we loop until convergence is reached. That is fine. BUT we do not want the program to end in an endless loop. Add another variable called max_it:usize and use this in a ‘for’ loop.

Let’s walk through:

Change the data type to the type we have in our read_data, the return value to a Vec<usize> and add a way to randomly select the initial set of centroids.

use ndarray::ArrayBase;
use rand::Rng;


fn kmeans(data: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>>, k: usize) -> Vec<usize> {
    let n = data.nrows();
    let d = data.ncols();
    let mut rng = rand::thread_rng();

    let mut centroids = Array::zeros((k, data.len_of(Axis(1))));
    for i in  0..k{
        let sample_idx = rng.gen_range(0..n);
        centroids.row_mut(i).assign(&self.data.row(sample_idx));
    }
    let mut labels = Vec<usize>::with_capacity(n);

    .
    .
    .
}

You see the creation of the initial labels changed quite a bit from the AI generated let mut centroids = data.slice(s![..k, ..]).to_owned();. We instead use a for loop that creates k random row ids and asignes the data from the respective row to the i’th centrioid.

Hence we now get exactly k random rows of the data and not just the first k rows as in the AI code. At least the part of how to get data from the data object into the centroids object is more or less a copy from ChatGPT.

Should we have more functions?

The next step is to calculate the closest centroid for each row of the data.

We will need this functionality also in the loop. As we use the same functionality more than once I would implement a function to do that. This function should compare the centroids to every row of the data and I would call it assign_labels:

fn assign_labels(data: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>>, centroids: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>> ) -> Vec<usize> {
    let n = data.nrows();
    let mut labels = Vec<usize>::with_capacity( n );
    
    for id in 0..n{
        let mut min_distance = f64::INFINITY;
        for j in 0..k {
            let distance = (data.slice(s![i, ..]) - centroids.slice(s![j, ..])).norm();
            if distance < min_distance {
                min_distance = distance;
                labels[i] = j;
            }
        }
    }
    labels
}

The assign_labels function is mainly a straight copy from the initial AI result. In short we use each data row once and calculate the distance to all centroids. If the distance we get is smaller than all others before we assign the centroid id to the lable for this data row.

Let’s add this logic to the main function:

fn kmeans(data: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>>, k: usize) -> Vec<usize> {
    let n = data.nrows();
    let d = data.ncols();
    let mut rng = rand::thread_rng();

    let mut centroids = Array::zeros((k, data.len_of(Axis(1))));
    for i in  0..k{
        let sample_idx = rng.gen_range(0..n);
        centroids.row_mut(i).assign(&self.data.row(sample_idx));
    }
    let mut labels = assign_labels( data, centroids );

    .
    .
    .
}

Add a function to calculate the new centroids

The next time we get the centroids will be more complicated. I like functions as one can also create specific tests for each function.

So lets add a “calculate_centroids” function:

fn calculate_centroids( data: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>>, labels:Vec<usize>, 
            k:usize, &mut centroids: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>> ) {

    for j in 0..k {
        let mask = labels.iter().enumerate()
            .filter(|(_, &l)| l == j )
            .map(|(i, _)| i)
            .collect::<Vec<_>>();
        let cluster_data = data.select(Axis(0), &mask);
        if !cluster_data.is_empty() {
            centroids.slice_mut(s![j, ..]).assign(&cluster_data.mean_axis(Axis(0)).unwrap());
        }
    } 
}

Oops - that one is complicated!

let mask = labels.iter().enumerate()
  .filter(|(_, &l)| l == j )
  .map(|(i, _)| i)
  .collect::<Vec<_>>(); 

Let’s break it down:

• labels.iter().enumerate() This creates the data we work on: The ids and variables in the label vector (.iter(): values and .enumerate(): with the ids ).

• filter(|(_, &l)| l == j ) This filter only uses the second variable, the label id and checks it against the cluster ‘j’. Only labels that have this id reach the next step.

• map(|(i, _)| i) Here we only get the id and value of rows that are in cluster ‘j’. So we simply need to ‘cut out’ the row id here and ignore the value.

So the mask variable in the end contains the row_ids that are in cluster j.

We collect this data into a new let cluster_data and if that is not empty we calculate the new centroid as the mean of this cluster_data. Simple?

In the calculate centroids I now also give the centroids object to the function. This way we do not need to create a temporary variable holding the centroids but instead change the values in the variable owned by the kmeans function. That is quite slick and we should probably also have done that in the assign_labels function. Would it even make sense there?

Let’s add this logic to the main function:

fn kmeans(data: &ArrayBase<ndarray::OwnedRepr<f64>, Dim<[usize; 2]>>, k: usize, max_it:usize) -> Vec<usize> {
    let n = data.nrows();
    let d = data.ncols();
    let mut rng = rand::thread_rng();

    let mut centroids = Array::zeros((k, data.len_of(Axis(1))));
    for i in  0..k{
        let sample_idx = rng.gen_range(0..n);
        centroids.row_mut(i).assign(&self.data.row(sample_idx));
    }
    let mut labels = assign_labels( data, centroids );
    let mut new_labels:Vec<usize>;

    for it in 0..max_it{
        calculate_centroids( data, labels, k, centroids );
        new_labels = assign_labels( data, centroids );
        if new_labels == labels{
            print("converged after {it} itherations");
            break;
        }
        labels = new_labels;
    }
    labels
}

And the TESTING!

As we now have so many nice functions that could be tested - lets add some tests:

To be able to test this we now neeed some data in this. :-( Let’s combine the kmeans functionality with the Data object we had as an example before.

Copy the read_data folder and name it kmeans. Now you need to change the Cargo.toml: read_data to kmeans.

You could compile it and check if it is still working.

Add the kmeans functionality to the Data class or to the main.rs script. This will take some time - I will help you, but there will be no help in this document.

We go on with the testing here:

#[cfg(test)]
mod tests {

    use crate::data::Data;
    use ndarray::Axis;
    use ndarray::s;

    #[test]
    fn check_assign_labels() {
        let data = Data::read_file( &"testData/CellexalVR_TestData_tsne.csv".to_string(), ',' );
        let centroids = data.data.slice(s![..10, ..]).to_owned();
        let kmeans = data.assign_labels( &centroids );
        let mut min = usize::MAX;
        let mut max = usize::MIN;
        for k in &kmeans {
            if min > *k{
                min = *k;
            }
            if max < *k{
                max = *k;
            }
        }
        assert_eq!( [min, max], [0,9] );
    }

     #[test]
    fn check_kmeans() {

        let mut data = Data::read_file( &"testData/CellexalVR_TestData_tsne.csv".to_string(), ',' );
        data.scale();
        let kmeans = data.kmeans( 15, 15000 );
        assert_eq!( kmeans.len(), data.data.len_of(Axis(0)) );

        //let exp:Vec<usize> = vec![0; data.data.len_of(Axis(0)) ];

        let mut min = usize::MAX;
        let mut max = usize::MIN;
        for k in kmeans {
            if min > k{
                min = k;
            }
            if max < k{
                max = k;
            }
        }
        assert_eq!( [min , max ], [0, 14] );
    }

}

Add these tests BEFORE you implement the functions. This way you have a simple test to check if it does work ;-)

Please try to integrate the kmeans function into the read_data crate.

You can either add the kmeans logics to the main.rs script or you augment the read_table Data class.

This will be a hands on training session. And not a too simple one. Please do not look into the kmeans_example folder. This would ruin the exercise for you.

Instead copy the read_data folder from the github repo and rename the project to kmeans_ in the Cargo.toml. Think about where you want to add the functions and do it. Compile and dig through the compiler errors. Or you try to clean out the problems first.

Compile and test the kmeans program

cd kmeans_example
cargo build -r 2>&1 | sed -r "s/\x1B\[([0-9]{1,2}(;[0-9]{1,2})?)?[m|K]//g"
target/release/kmeans -h

##     Finished release [optimized] target(s) in 0.02s
## kmeans 1.0.0
## Stefan L. <stefan.lang@med.lu.se>
## Calculate a kmeans clustering based on a numeric matrix
## 
## USAGE:
##     kmeans [OPTIONS] --clusters <CLUSTERS>
## 
## OPTIONS:
##     -c, --clusters <CLUSTERS>    the number of clusters
##     -d, --data <DATA>            the data (text file) [default:
##                                  testData/Spellman_Yeast_Cell_Cycle.tsv]
##     -h, --help                   Print help information
##     -m, --max-it <MAX_IT>        the max iterations [default: 1000000]
##     -o, --outfile <OUTFILE>      the grouping outfile [default: testData/Clustering.txt]
##     -s, --sep <SEP>              the column separator for the file [default: \t]
##     -V, --version                Print version information

cd kmeans_example
target/release/kmeans -d testData/CellexalVR_TestData_tsne.csv -c 15 -s ','

## I got 1654 rows and 3 cols in this data
## finished after 17 iterations
## finished in 0 h 0 min 0 sec 4 milli sec

Ok that did not take long. But how fast would R do that?

Did this make sense? Let’s plot this data using R:

clusters = read.delim('kmeans_example/testData/Clustering.txt', sep=",", row.names=1, header=F)
head(clusters)

##          V2
## HSPC_001  6
## HSPC_002  6
## HSPC_003  9
## HSPC_004  5
## HSPC_006 12
## HSPC_008  1

data = read.delim( 'kmeans_example/testData/CellexalVR_TestData_tsne.csv', sep=",", row.names=1)
head(data)

##                V1          V2         V3
## HSPC_001 24.35282  12.8947641  3.8425078
## HSPC_002 18.38284  12.0568382  1.5107061
## HSPC_003  2.78855 -12.6130571 -6.8210608
## HSPC_004 24.57608   0.5695849 -0.9598198
## HSPC_006 16.85040 -16.1784775  3.2230418
## HSPC_008 26.04562   2.7516960  7.5135636

plot( data[,1:2], col=rainbow(15)[clusters[,1]+1], pch=16)

But what about speed? Can we compare that to R? First check R like we did it in the R programming 1:

# from ChatGPT:

e.dist <- function(x, y) {
  sqrt(sum((x - y)^2))
}


kmeans_r = function( yyc, k=8 ){
    n = nrow(yyc)
    centroids <- yyc[sample(1:n, k, replace=F),] # make a vector of random clusters.

    dists <- matrix(0, nrow = n, ncol = k) # make an empty matrix to fill with distances.

    clusters <- NULL # make an empty variable to catch the clusters in the loop below
    old_clusters <- rep(0, n)
    for(iteration in 1:100){ # 100 iterations

      for(gene in 1:n){
        for(cl in 1:k){
          dists[gene, cl] <- e.dist(yyc[gene,], centroids[cl,]) # for each gene calculate the distance to each centroid (cl).
        }
      }

      clusters <- apply(dists, 1, which.min) # assign a  cluster according to which centroid is nearest
      if ( all.equal(old_clusters, clusters) ==TRUE ){
        return (clusters);
      }
      old_clusters =clusters

      for(cl in 1:k){
        centroids[cl,] <- apply(yyc[which(clusters == cl),], 2, mean) # define new centroids.
      }
    }
    clusters
}

Let’s time that call:

start_time <- Sys.time()
#data = read.delim( 'kmeans_example/testData/CellexalVR_TestData_tsne.csv', sep=",", row.names=1)
#data = read.delim( './Spellman_Yeast_Cell_Cycle.tsv', sep=",", row.names=1)
df <- data.frame(age = c(18, 21, 22, 24, 26, 26, 27, 30, 31, 35, 39, 40, 41, 42, 44, 46, 47, 48, 49, 54),
    spend = c(10, 11, 22, 15, 12, 13, 14, 33, 39, 37, 44, 27, 29, 20, 28, 21, 30, 31, 23, 24)
)
cl = kmeans_r( df, 3)
end_time <- Sys.time()
print(end_time - start_time)

## Time difference of 0.1014931 secs

plot( df[,1:2], col=rainbow(15)[cl], pch=16)

OK even this tiny example (20 rows) took more time than what the Rust program needed for a much bigger matrix (1654 rows). So this simple R code is far slower than the Rust code. But what about the kmeans implementation in R? This should be highly optimized.

Does the R internal kmeans (likely programmed in C) outperform Rust?

start_time <- Sys.time()
data = read.delim( 'kmeans_example/testData/CellexalVR_TestData_tsne.csv', sep=",", row.names=1)
cl = kmeans( data, 15)
end_time <- Sys.time()
print(end_time - start_time)

## Time difference of 0.006499529 secs

plot( data[,1:2], col=rainbow(15)[cl$cluster], pch=16)

Take home message:

Rust totally outclasses simple R code and ends up in the same speed range if optimized (c or c++) code is used in R.

But trust me - you do not want to write optimized c or c++ code!

How would we now implements a DNA mapper like STAR?

Fist we would need to see how nucleotide sequences can be handled in Rust. For this we define that one single DNA nucleotide can either be ‘A’, ‘C’, ‘G’ or ‘T’ - in other words 4 different options. So - we could also code for them using 2 bits ‘A’ = 00, ‘C’ = 01, ‘G’ = 10, ‘T’ = 11.

Each char tahes a u8 to store it internally. That is 8 bits wide and as our very limited DNA code can be coded in only 2 bits this should be 4x as efficient - a single u8 can contain the same info as a 4bp DNA code.

With being an integer this representation of a DNA sequence opens up possibilities that you do simply not have in R.

pub type Base = u8;
pub const A: Base = 0;
pub const C: Base = 1;
pub const G: Base = 2;
pub const T: Base = 3;
pub fn encode_binary(c: char) -> Base {
    match c {
        'A' | 'a' => A,
        'C' | 'c' => C,
        'G' | 'g' => G,
        'T' | 't' => T,
        _ => panic!("cannot decode {c} into 2 bit encoding"),
    }
}

let s = String::from("AGCTACGT");
println!("How Rust can look at the sequence {s}");
let t = s.as_bytes();
println!("as_bytes {t:?}");
let st = std::str::from_utf8(t).unwrap();
println!("as str {st:?}");

let mut w = 0_u8;
for c in st.chars().rev() {
  w <<= 2;
  w |= encode_binary(c);
}
println!("and as single u8 integer {w}");
println!("which looks like this in binary form: {w:b}");

let faster =std::str::from_utf8( b"AGCTAGCT").unwrap();

for c in faster.chars().rev() {
  w <<= 2;
  w |= encode_binary(c);
}
println!("and now in one line from string to binary: {w:b}");

What does this do? In the beginning it defined a Base variable of class type. This variable is later used to define which kind of data you can process using this logic (e.g. 8, 16, 32 or even 64 bp of sequence info). The next four lines define which sequences should be converted into which ‘number’. It helps tremendously if you try to think of 0-3 as binary numbers: #00, #01, #10 and #11.

“Defining something using const instead of let can help make the code more reliable, efficient, and maintainable by enforcing immutability and compile-time evaluation, and by providing global scope and better type checking.” For me this translates to - just do it this way - it does certainly not hurt.

The ‘let part’ of the script defines the string we want to convert, converts the string to a bytes vector ([u8]). The for loop finally converts the str into a Base integer. Can you see what it does - does the initial string really get represented in the u8? Do 8 bp really fit into this u8? And should that not throw an error?

cd sequences
cargo build -r 2>&1 | sed -r "s/\x1B\[([0-9]{1,2}(;[0-9]{1,2})?)?[m|K]//g"
target/release/sequences

##     Finished release [optimized] target(s) in 0.00s
## How Rust can look at the sequence AGCTACGT
## as_bytes [65, 71, 67, 84, 65, 67, 71, 84]
## as str "AGCTACGT"
## and as single u8 integer 216
## which looks like this in binary form: 11011000
## and now in one line from string to binary: 11011000

As said - this is not scope of this course, but I assume it will become interesting for you later on.

So what about we try to implement a mapper using this knowledge? This will likely be a horrible experience for you :-D

Implement the DNA mapper

Each of the coding examples have been one Rust project so far. But in this project we possibly need more binaries than only one. We will also create more library files and therfore we need to separate the binaries from the library. In Rust you do that like this:

cargo new my_mapper 2>&1 | sed -r "s/\x1B\[([0-9]{1,2}(;[0-9]{1,2})?)?[m|K]//g"
cd my_mapper
mkdir src/bin
touch src/lib.rs
mv src/main.rs src/bin/first_executable.rs

find ./

##      Created binary (application) `my_mapper` package
## ./
## ./src
## ./src/main.rs
## ./src/bin
## ./src/bin/first_executable.rs
## ./src/lib.rs
## ./Cargo.toml

Each file in the bin folder will code for a new binary.

Now we need one class that stores the target sequences as 2bit, one that makes them searchable and one or two that implement the glue between the two. Why three or four classes? Trust me we will need it!

The DNA data class

Before we have learned how to convert a DNA string into a u8 integer. But we can only store 4bp in one u8. How do we store longer sequences? A vector of u8s?

Which functions does that class need?

1. We definitely need an encode() function that converts a UTF8 &[u8] (that is what a str internally is in Rust) into our own 2bit version of Vec::<u8>.
2. The matching DNA fragmanets likely do only match a sub sequence of the here store DNA sequence -> we need a slice() function.
3. Reading 2bit encoded DNA sequences is a pain. We likely need a decode and a print function (e.g. for debug).
4. Our other class needs some keys it can use to map to this DNA sequence - we need to get an iterator implemented here that can create these keys.

With this fourth function you now see why we wanted two classes. If we have two classes - one for the matcher and one for the DNA data we can use the DNA data class for both the target sequences as well as for the reads. Nice.

Should we call it DnaData hence forth?

And our mapper class?

That thing needs to store our indexed DNA fragments and have a fast way to find them.

First the data structures we need:

1. The indexed DNA fragments are of class DnaData and we only need to fill that once - so and array/vector would be ok.
2. Each key from the indexed DNA fragments needs to link to the DNA fragment the key came from and the position on the fragment that the key did start. Oh - do we possibly need three classes? An additional IndexElement?

There are multiple ways to define the data structure for the index keys. But if we remember that we (1) have not defined the keys and (2) each of our DNA fragments is internally stored as integer - why not use a vector, too?

Each x bp DNA fragment would be seen as an integer and represent the array position the data is stored in.

Lets see - 4bp (8 bits) only code for 256 different entries. 256 * 8 byte is no memory at all, but there is not enough variability here. 16bp (32 bit) - that codes for 4.294.967.296 differnt DNA strings - quite cool, but does that fit into memory? (2 ** 32 * 8byte /1024 /1024 / 1024 = 32 Gb). Would likely fit, but that seams too excessive. So what about a u16 - 8bp or 16 bit: 2 ** 32 * 8byte /1024 /1024 / 1024 = 0.5 Gb for 32.768 different 8bp DNA strings.

So that seams to be the best option nowadays? This defines the iterator’s return value: u16.

And finally the IndexElement

Why do we need that? The Mapper needs to store information on all the different u16 indices. What it needs is a storage for the id of the DnaData and the startposition on that data. And as all of these 8bp DNA regions can of cause also occur multiple times in the genome we need to store multiple of these (dna_data_id, start) touples. This touple would be our fourth class. We likely need a length function and Rust also wants an is_empty to that.

There will be more when we actually implement this mapper.

Final Remarks

I think implementing a mapper here would likely blow the scope. But you get the idea - right?

I have a GenomeMapper class in my Rustody library (https://github.com/stela2502/Rustody) that more or less implements this logics. You can take a look and use that if you need it.

The Rustody project also used Traits (Rust way of inheritance) and a folder structure for the library itself. I hope you learn something reading it - if needed.