add github mirroring

Signed-off-by: Stephen Waits <steve@waits.net>

.build.yml

@@ -0,0 +1,11 @@
+image: alpine/edge
+secrets:
+  - 88cce185-d2ec-4573-9546-e33f36f79bbf
+sources:
+  - git+ssh://git@git.sr.ht/~swaits/megapermute
+tasks:
+  - mirror-to-github: |
+      cd ~/megapermute
+      ssh-keyscan -t rsa github.com >> ~/.ssh/known_hosts
+      git remote add github git@github.com:swaits/megapermute.git
+      git push --mirror github

.gitignore

@@ -1 +1,21 @@
 /target
+#--------------------------------------------------#
+# The following was generated with gitignore.nvim: #
+#--------------------------------------------------#
+# Gitignore for the following technologies: Rust
+
+# Generated by Cargo
+# will have compiled files and executables
+debug/
+target/
+
+# Remove Cargo.lock from gitignore if creating an executable, leave it for libraries
+# More information here https://doc.rust-lang.org/cargo/guide/cargo-toml-vs-cargo-lock.html
+Cargo.lock
+
+# These are backup files generated by rustfmt
+**/*.rs.bk
+
+# MSVC Windows builds of rustc generate these, which store debugging information
+*.pdb
+

Cargo.lock

@@ -4,15 +4,15 @@ version = 3
 
 [[package]]
 name = "anyhow"
-version = "1.0.52"
+version = "1.0.75"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "84450d0b4a8bd1ba4144ce8ce718fbc5d071358b1e5384bace6536b3d1f2d5b3"
+checksum = "a4668cab20f66d8d020e1fbc0ebe47217433c1b6c8f2040faf858554e394ace6"
 
 [[package]]
 name = "autocfg"
-version = "1.0.1"
+version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "cdb031dd78e28731d87d56cc8ffef4a8f36ca26c38fe2de700543e627f8a464a"
+checksum = "d468802bab17cbc0cc575e9b053f41e72aa36bfa6b7f55e3529ffa43161b97fa"
 
 [[package]]
 name = "cfg-if"
@@ -20,21 +20,11 @@ version = "1.0.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
 
-[[package]]
-name = "crossbeam-channel"
-version = "0.5.2"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e54ea8bc3fb1ee042f5aace6e3c6e025d3874866da222930f70ce62aceba0bfa"
-dependencies = [
- "cfg-if",
- "crossbeam-utils",
-]
-
 [[package]]
 name = "crossbeam-deque"
-version = "0.8.1"
+version = "0.8.3"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "6455c0ca19f0d2fbf751b908d5c55c1f5cbc65e03c4225427254b46890bdde1e"
+checksum = "ce6fd6f855243022dcecf8702fef0c297d4338e226845fe067f6341ad9fa0cef"
 dependencies = [
  "cfg-if",
  "crossbeam-epoch",
@@ -43,64 +33,48 @@ dependencies = [
 
 [[package]]
 name = "crossbeam-epoch"
-version = "0.9.6"
+version = "0.9.15"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "97242a70df9b89a65d0b6df3c4bf5b9ce03c5b7309019777fbde37e7537f8762"
+checksum = "ae211234986c545741a7dc064309f67ee1e5ad243d0e48335adc0484d960bcc7"
 dependencies = [
+ "autocfg",
  "cfg-if",
  "crossbeam-utils",
- "lazy_static",
  "memoffset",
  "scopeguard",
 ]
 
 [[package]]
 name = "crossbeam-utils"
-version = "0.8.6"
+version = "0.8.16"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "cfcae03edb34f947e64acdb1c33ec169824e20657e9ecb61cef6c8c74dcb8120"
+checksum = "5a22b2d63d4d1dc0b7f1b6b2747dd0088008a9be28b6ddf0b1e7d335e3037294"
 dependencies = [
  "cfg-if",
- "lazy_static",
 ]
 
 [[package]]
 name = "either"
-version = "1.6.1"
+version = "1.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457"
+checksum = "a26ae43d7bcc3b814de94796a5e736d4029efb0ee900c12e2d54c993ad1a1e07"
 
 [[package]]
 name = "getrandom"
-version = "0.2.3"
+version = "0.2.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "7fcd999463524c52659517fe2cea98493cfe485d10565e7b0fb07dbba7ad2753"
+checksum = "c85e1d9ab2eadba7e5040d4e09cbd6d072b76a557ad64e797c2cb9d4da21d7e4"
 dependencies = [
  "cfg-if",
  "libc",
  "wasi",
 ]
 
-[[package]]
-name = "hermit-abi"
-version = "0.1.19"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33"
-dependencies = [
- "libc",
-]
-
-[[package]]
-name = "lazy_static"
-version = "1.4.0"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646"
-
 [[package]]
 name = "libc"
-version = "0.2.112"
+version = "0.2.144"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "1b03d17f364a3a042d5e5d46b053bbbf82c92c9430c592dd4c064dc6ee997125"
+checksum = "2b00cc1c228a6782d0f076e7b232802e0c5689d41bb5df366f2a6b6621cfdfe1"
 
 [[package]]
 name = "megapermute"
@@ -113,39 +87,28 @@ dependencies = [
 
 [[package]]
 name = "memoffset"
-version = "0.6.5"
+version = "0.9.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "5aa361d4faea93603064a027415f07bd8e1d5c88c9fbf68bf56a285428fd79ce"
+checksum = "5a634b1c61a95585bd15607c6ab0c4e5b226e695ff2800ba0cdccddf208c406c"
 dependencies = [
  "autocfg",
 ]
 
-[[package]]
-name = "num_cpus"
-version = "1.13.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "19e64526ebdee182341572e50e9ad03965aa510cd94427a4549448f285e957a1"
-dependencies = [
- "hermit-abi",
- "libc",
-]
-
 [[package]]
 name = "ppv-lite86"
-version = "0.2.16"
+version = "0.2.17"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "eb9f9e6e233e5c4a35559a617bf40a4ec447db2e84c20b55a6f83167b7e57872"
+checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de"
 
 [[package]]
 name = "rand"
-version = "0.8.4"
+version = "0.8.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "2e7573632e6454cf6b99d7aac4ccca54be06da05aca2ef7423d22d27d4d4bcd8"
+checksum = "34af8d1a0e25924bc5b7c43c079c942339d8f0a8b57c39049bef581b46327404"
 dependencies = [
  "libc",
  "rand_chacha",
  "rand_core",
- "rand_hc",
 ]
 
 [[package]]
@@ -160,55 +123,41 @@ dependencies = [
 
 [[package]]
 name = "rand_core"
-version = "0.6.3"
+version = "0.6.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d34f1408f55294453790c48b2f1ebbb1c5b4b7563eb1f418bcfcfdbb06ebb4e7"
+checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c"
 dependencies = [
  "getrandom",
 ]
 
-[[package]]
-name = "rand_hc"
-version = "0.3.1"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d51e9f596de227fda2ea6c84607f5558e196eeaf43c986b724ba4fb8fdf497e7"
-dependencies = [
- "rand_core",
-]
-
 [[package]]
 name = "rayon"
-version = "1.5.1"
+version = "1.8.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "c06aca804d41dbc8ba42dfd964f0d01334eceb64314b9ecf7c5fad5188a06d90"
+checksum = "9c27db03db7734835b3f53954b534c91069375ce6ccaa2e065441e07d9b6cdb1"
 dependencies = [
- "autocfg",
- "crossbeam-deque",
  "either",
  "rayon-core",
 ]
 
 [[package]]
 name = "rayon-core"
-version = "1.9.1"
+version = "1.12.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d78120e2c850279833f1dd3582f730c4ab53ed95aeaaaa862a2a5c71b1656d8e"
+checksum = "5ce3fb6ad83f861aac485e76e1985cd109d9a3713802152be56c3b1f0e0658ed"
 dependencies = [
- "crossbeam-channel",
  "crossbeam-deque",
  "crossbeam-utils",
- "lazy_static",
- "num_cpus",
 ]
 
 [[package]]
 name = "scopeguard"
-version = "1.1.0"
+version = "1.2.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd"
+checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49"
 
 [[package]]
 name = "wasi"
-version = "0.10.2+wasi-snapshot-preview1"
+version = "0.11.0+wasi-snapshot-preview1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "fd6fbd9a79829dd1ad0cc20627bf1ed606756a7f77edff7b66b7064f9cb327c6"
+checksum = "9c8d87e72b64a3b4db28d11ce29237c246188f4f51057d65a7eab63b7987e423"

Cargo.toml

@@ -6,6 +6,6 @@ edition = "2021"
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
-anyhow = "1.0.52"
-rand = "0.8.4"
-rayon = "1.5.1"
+anyhow = "1.0.75"
+rand = "0.8.5"
+rayon = "1.8.0"

README.md

@@ -2,14 +2,109 @@
 
 A million permutation hypothesis test utility.
 
+## Installation
+
+1. If you don't already have rust installed, [install it](https://www.rust-lang.org/tools/install).
+2. Then run: `cargo install --git https://git.sr.ht/~swaits/megapermute`
+
 ## Usage
 
-`megapermute` requires two input files, `control.dat` and `treatment.dat`. Each is a text
-file containing a list of numbers, treated internally as `f64`.
+`megapermute` is a command-line utility which reads two input files,
+`control.dat` and `treatment.dat`. Each is a text file containing a list of
+numbers, treated internally as 64-bit floating point numbers.
 
-- `control.dat` the control sample
-- `treatment.dat` the treatment sample
+- `control.dat`: the control sample
+- `treatment.dat`: the treatment sample
 
-With these two files in the current directly, running `megapermute` computes
-the empircal difference in means and p-value for the null hypothesis (H_0 ==
-H_a).
+With these two files in the current directory, running `megapermute` computes
+the empircal difference in means and p-value for the null hypothesis, that H_0
+is equivalent to H_a. You can think of p-value like the probability that the
+null hypothesis is true.
+
+### Example
+
+Given the following two input files:
+
+```bash
+~/tmp
+❯ bat control.dat
+───────┬───────────────────────
+       │ File: control.dat
+───────┼───────────────────────
+   1   │ 52
+   2   │ 104
+   3   │ 146
+   4   │ 10
+   5   │ 51
+   6   │ 30
+   7   │ 40
+   8   │ 27
+   9   │ 46
+───────┴───────────────────────
+
+~/tmp
+❯ bat treatment.dat
+───────┬───────────────────────
+       │ File: treatment.dat
+───────┼───────────────────────
+   1   │ 94
+   2   │ 197
+   3   │ 16
+   4   │ 38
+   5   │ 99
+   6   │ 141
+   7   │ 23
+───────┴──────────────────────
+```
+
+Running `megapermute` should yield results like this:
+
+```bash
+~/tmp
+❯ megapermute
+                 mu_control = 56.22222222222222
+               mu_treatment = 86.85714285714286
+(mu_treatment - mu_control) = 30.63492063492064
+                    p-value = 0.139556
+```
+
+## Performance
+
+My goal with this was to take advantage of multiple cores by using
+[`rayon`](https://crates.io/crates/rayon) iterators to compute the permutation
+hypothesis test quickly. As of January 2022, this application can do one
+million permutations on the example dataset listed above in 39.4 ms (N=69,
+SD=2.6 ms), which is on the order of 25M permutations/second.
+
+```bash
+~/tmp
+❯ hyperfine megapermute
+Benchmark 1: megapermute
+  Time (mean ± σ):      39.4 ms ±   2.6 ms    [User: 499.3 ms, System: 10.2 ms]
+  Range (min … max):    35.1 ms …  48.4 ms    69 runs
+```
+
+Benchmark machine info:
+
+```bash
+~/tmp
+❯ neofetch
+                    'c.          stephenwaits@SWaits-Mac-01.local
+                 ,xNMM.          --------------------------------
+               .OMMMMo           OS: macOS 11.6.1 20G224 x86_64
+               OMMM0,            Host: MacBookPro16,1
+     .;loddo:' loolloddol;.      Kernel: 20.6.0
+   cKMMMMMMMMMMNWMMMMMMMMMM0:    Uptime: 31 days, 21 hours, 23 mins
+ .KMMMMMMMMMMMMMMMMMMMMMMMWd.    Packages: 149 (brew)
+ XMMMMMMMMMMMMMMMMMMMMMMMX.      Shell: zsh 5.8
+;MMMMMMMMMMMMMMMMMMMMMMMM:       Resolution: 3840x2160, 1920x1080
+:MMMMMMMMMMMMMMMMMMMMMMMM:       DE: Aqua
+.MMMMMMMMMMMMMMMMMMMMMMMMX.      WM: Amethyst
+ kMMMMMMMMMMMMMMMMMMMMMMMMWd.    Terminal: WezTerm
+ .XMMMMMMMMMMMMMMMMMMMMMMMMMMk   CPU: Intel i9-9880H (16) @ 2.30GHz
+  .XMMMMMMMMMMMMMMMMMMMMMMMMK.   GPU: Intel UHD Graphics 630, AMD Radeon Pro 5500M
+    kMMMMMMMMMMMMMMMMMMMMMMd     Memory: 21514MiB / 32768MiB
+     ;KMMMMMMMWXXWMMMMMMMk.
+       .cooc,.    .,coo:.
+
+```

src/main.rs

@@ -14,30 +14,62 @@ const N_THREADS: usize = 1_000;
 // the number of permutation tests to run per "thread" (above)
 const N_PERMUTATIONS_PER_THREAD: usize = 1_000;
 
+// load a file of numbers as f64 and panic if we run into any problems.
+// This function takes a filename as a string and returns a vector of f64s.
+// If there is an error, it will return a Result with an error message.
 fn load_f64s(filename: &str) -> Result<Vec<f64>> {
+    // Open the file and read it into a buffer.
     let f = BufReader::new(File::open(filename)?);
+    // Map each line of the file to a f64.
     Ok(f.lines()
         .map(|l| {
+            // If there is an error reading the line, panic.
             l.expect("input error: unable to read line")
+                // If there is an error parsing the line, panic.
                 .parse::<f64>()
                 .expect("parse error: unable to parse f64")
         })
+        // Collect the f64s into a vector.
         .collect())
 }
 
-// accepts an iterator of f64's and computes mean using Welford's online algorithm
+// This function accepts an iterator of f64's and computes mean using Welford's online algorithm
 fn mean<'a>(iter: impl Iterator<Item = &'a f64>) -> f64 {
+    // The function uses the enumerate method to get the index and value of each element in the iterator
+    // The fold method is used to iterate through the iterator and add the values to the accumulator
+    // The accumulator is initialized to 0.0
+    // The accumulator is updated by adding the difference between the current value and the accumulator divided by the index plus 1
     iter.enumerate()
         .fold(0.0, |mu, (i, x)| mu + ((x - mu) / (i + 1) as f64))
 }
 
+// this enum is used to denote group memebership during each permutation
 #[derive(Clone, PartialEq, Eq)]
 enum Group {
     Control,
     Treatment,
 }
 
-fn permute(control: &[f64], treatment: &[f64], mu_diff: f64) -> f64 {
+// Run the permutation tests.
+//
+// Accepts the control and treatment arrays along with the difference in empircal means (treatment
+// minus control).
+//
+// Runs N_PERMUTATIONS_PER_THREAD on each of N_THREADS using Rayon's `par_iter()`.
+//
+// The original data (arrays) is never copied. Each of the iterations creates an index array of
+// `enum Group` denoting `Control` or `Treatment` at each index. That `index` array is shuffled to
+// permute group memebership. Finally, each of the group's means are computed and compared, and
+// filtered to only count differences which are larger than the `mu_diff` parameter.
+//
+// The final p-value is the number of permutations where the diff in means exceeded `mu_diff`
+// divided by (N_THREADS * N_PERMUTATIONS_PER_THREAD).
+//
+// Note: the left tail or right tail is chosen automatically based on whether the empircal mean
+// delta is positive or negative.
+fn permutation_test(control: &[f64], treatment: &[f64], mu_diff: f64) -> f64 {
+    // use Rayon to divide this work across N_THREADS, ultimately counting the number of
+    // permutations where delta(means) > mu_diff
     let count: f64 = (0..N_THREADS)
         .into_par_iter()
         .map(|_| {
@@ -57,39 +89,28 @@ fn permute(control: &[f64], treatment: &[f64], mu_diff: f64) -> f64 {
                     // shuffle our group selections
                     index.shuffle(&mut rng);
 
-                    // comput emean of this "control permutation"
-                    let permuted_mean_control = mean(
-                        control
-                            .iter()
-                            .chain(treatment.iter())
-                            .enumerate()
-                            .filter_map(|(i, x)| {
-                                if index[i] == Group::Control {
-                                    Some(x)
-                                } else {
-                                    None
-                                }
-                            }),
-                    );
+                    // variables for mean computation
+                    let (mut mu_control, mut n_control) = (0.0, 0.0);
+                    let (mut mu_treatment, mut n_treatment) = (0.0, 0.0);
 
-                    // comput emean of this "treatment permutation"
-                    let permuted_mean_treatment = mean(
-                        control
-                            .iter()
-                            .chain(treatment.iter())
-                            .enumerate()
-                            .filter_map(|(i, x)| {
-                                if index[i] == Group::Treatment {
-                                    Some(x)
-                                } else {
-                                    None
-                                }
-                            }),
+                    // walk the combined iterator and add each element to the corresponding mean
+                    // using Welford's online algorithm
+                    control.iter().chain(treatment.iter()).enumerate().for_each(
+                        |(i, x)| match index[i] {
+                            Group::Control => {
+                                n_control += 1.0;
+                                mu_control += (x - mu_control) / n_control
+                            }
+                            Group::Treatment => {
+                                n_treatment += 1.0;
+                                mu_treatment += (x - mu_treatment) / n_treatment
+                            }
+                        },
                     );
 
                     // select this permutation if the diff in these permuted means is more than the
                     // empirical diff of means
-                    (permuted_mean_treatment - permuted_mean_control) > mu_diff
+                    (mu_treatment - mu_control) > mu_diff
                 })
                 .count() as f64
         })
@@ -107,6 +128,18 @@ fn permute(control: &[f64], treatment: &[f64], mu_diff: f64) -> f64 {
     }
 }
 
+// convert p-value to conventional language
+fn pvalue_to_string(p: f64) -> String {
+    match p {
+        p if p < 0.01 => "very strong evidence against null hypothesis",
+        p if p < 0.025 => "strong evidence against null hypothesis",
+        p if p < 0.05 => "reasonably strong evidence against null hypothesis",
+        p if p < 0.10 => "borderline evidence against null hypothesis",
+        _ => "no evidence against null hypothesis",
+    }
+    .to_string()
+}
+
 fn main() -> Result<()> {
     // read data
     let control = load_f64s("control.dat")?;
@@ -116,15 +149,18 @@ fn main() -> Result<()> {
     let mean_control = mean(control.iter());
     let mean_treatment = mean(treatment.iter());
     println!("                 mu_control = {}", mean_control);
+    println!("                  N_control = {}", control.len());
     println!("               mu_treatment = {}", mean_treatment);
+    println!("                N_treatment = {}", treatment.len());
     println!(
         "(mu_treatment - mu_control) = {}",
         mean_treatment - mean_control
     );
 
     // run permutation test to compute p-value
-    let pvalue = permute(&control, &treatment, mean_treatment - mean_control);
+    let pvalue = permutation_test(&control, &treatment, mean_treatment - mean_control);
     println!("                    p-value = {}", pvalue);
+    println!("                     result = {}", pvalue_to_string(pvalue));
 
     Ok(())
 }
@@ -142,6 +178,8 @@ mod tests {
         (a - b).abs() < epsilon
     }
 
+    // This is a test file for the permutation test.
+    // It tests the permutation test on the mouse data from Table 2.1 in "An Introduction to the Bootstrap" (book)
     #[test]
     fn test_mouse_data() {
         // test data from Table 2.1 in "An Introduction to the Bootstrap" (book)
@@ -160,7 +198,7 @@ mod tests {
         ));
 
         // run permutation test to compute p-value
-        let pvalue = permute(&control, &treatment, mean_treatment - mean_control);
+        let pvalue = permutation_test(&control, &treatment, mean_treatment - mean_control);
         assert!(approx_eq(pvalue, 0.13896357, PVALUE_EPSILON));
     }
 }