From bab91b3183b5c83a3b329b9d20033525b556bf00 Mon Sep 17 00:00:00 2001 From: Stephen Waits Date: Mon, 10 Jan 2022 13:54:47 -0700 Subject: [PATCH] initial commit --- .gitignore | 1 + Cargo.lock | 214 ++++++++++++++++++++++++++++++++++++++++++++++++++++ Cargo.toml | 11 +++ LICENSE.txt | 21 ++++++ README.md | 15 ++++ src/main.rs | 166 ++++++++++++++++++++++++++++++++++++++++ tmp | 3 + 7 files changed, 431 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 LICENSE.txt create mode 100644 README.md create mode 100644 src/main.rs create mode 100644 tmp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ea8c4bf --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +/target diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..efc4ab2 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,214 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "anyhow" +version = "1.0.52" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "84450d0b4a8bd1ba4144ce8ce718fbc5d071358b1e5384bace6536b3d1f2d5b3" + +[[package]] +name = "autocfg" +version = "1.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cdb031dd78e28731d87d56cc8ffef4a8f36ca26c38fe2de700543e627f8a464a" + +[[package]] +name = "cfg-if" +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" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6455c0ca19f0d2fbf751b908d5c55c1f5cbc65e03c4225427254b46890bdde1e" +dependencies = [ + "cfg-if", + "crossbeam-epoch", + "crossbeam-utils", +] + +[[package]] +name = "crossbeam-epoch" +version = "0.9.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "97242a70df9b89a65d0b6df3c4bf5b9ce03c5b7309019777fbde37e7537f8762" +dependencies = [ + "cfg-if", + "crossbeam-utils", + "lazy_static", + "memoffset", + "scopeguard", +] + +[[package]] +name = "crossbeam-utils" +version = "0.8.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "cfcae03edb34f947e64acdb1c33ec169824e20657e9ecb61cef6c8c74dcb8120" +dependencies = [ + "cfg-if", + "lazy_static", +] + +[[package]] +name = "either" +version = "1.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457" + +[[package]] +name = "getrandom" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7fcd999463524c52659517fe2cea98493cfe485d10565e7b0fb07dbba7ad2753" +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" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1b03d17f364a3a042d5e5d46b053bbbf82c92c9430c592dd4c064dc6ee997125" + +[[package]] +name = "megapermute" +version = "0.1.0" +dependencies = [ + "anyhow", + "rand", + "rayon", +] + +[[package]] +name = "memoffset" +version = "0.6.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5aa361d4faea93603064a027415f07bd8e1d5c88c9fbf68bf56a285428fd79ce" +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" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb9f9e6e233e5c4a35559a617bf40a4ec447db2e84c20b55a6f83167b7e57872" + +[[package]] +name = "rand" +version = "0.8.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2e7573632e6454cf6b99d7aac4ccca54be06da05aca2ef7423d22d27d4d4bcd8" +dependencies = [ + "libc", + "rand_chacha", + "rand_core", + "rand_hc", +] + +[[package]] +name = "rand_chacha" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" +dependencies = [ + "ppv-lite86", + "rand_core", +] + +[[package]] +name = "rand_core" +version = "0.6.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d34f1408f55294453790c48b2f1ebbb1c5b4b7563eb1f418bcfcfdbb06ebb4e7" +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" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c06aca804d41dbc8ba42dfd964f0d01334eceb64314b9ecf7c5fad5188a06d90" +dependencies = [ + "autocfg", + "crossbeam-deque", + "either", + "rayon-core", +] + +[[package]] +name = "rayon-core" +version = "1.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d78120e2c850279833f1dd3582f730c4ab53ed95aeaaaa862a2a5c71b1656d8e" +dependencies = [ + "crossbeam-channel", + "crossbeam-deque", + "crossbeam-utils", + "lazy_static", + "num_cpus", +] + +[[package]] +name = "scopeguard" +version = "1.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" + +[[package]] +name = "wasi" +version = "0.10.2+wasi-snapshot-preview1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "fd6fbd9a79829dd1ad0cc20627bf1ed606756a7f77edff7b66b7064f9cb327c6" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..6c3d49b --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,11 @@ +[package] +name = "megapermute" +version = "0.1.0" +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" diff --git a/LICENSE.txt b/LICENSE.txt new file mode 100644 index 0000000..180e548 --- /dev/null +++ b/LICENSE.txt @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2022 Stephen Waits + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..5f2dc87 --- /dev/null +++ b/README.md @@ -0,0 +1,15 @@ +# megapermute + +A million permutation hypothesis test utility. + +## 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`. + +- `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). diff --git a/src/main.rs b/src/main.rs new file mode 100644 index 0000000..dfb0ce5 --- /dev/null +++ b/src/main.rs @@ -0,0 +1,166 @@ +use std::{ + fs::File, + io::{BufRead, BufReader}, + iter, +}; + +use anyhow::Result; +use rand::{prelude::SliceRandom, thread_rng}; +use rayon::prelude::*; + +// the number of threads to pass to rayon (which smartly splits it into actual number of threads) +const N_THREADS: usize = 1_000; + +// the number of permutation tests to run per "thread" (above) +const N_PERMUTATIONS_PER_THREAD: usize = 1_000; + +fn load_f64s(filename: &str) -> Result> { + let f = BufReader::new(File::open(filename)?); + Ok(f.lines() + .map(|l| { + l.expect("input error: unable to read line") + .parse::() + .expect("parse error: unable to parse f64") + }) + .collect()) +} + +// accepts an iterator of f64's and computes mean using Welford's online algorithm +fn mean<'a>(iter: impl Iterator) -> f64 { + iter.enumerate() + .fold(0.0, |mu, (i, x)| mu + ((x - mu) / (i + 1) as f64)) +} + +#[derive(Clone, PartialEq, Eq)] +enum Group { + Control, + Treatment, +} + +fn permute(control: &[f64], treatment: &[f64], mu_diff: f64) -> f64 { + let count: f64 = (0..N_THREADS) + .into_par_iter() + .map(|_| { + // create an rng so we can shuffle later + let mut rng = thread_rng(); + + // create an index array of [Control, Treatment] which we'll shuffle repeatedly to make our + // selections during each permutation + let mut index: Vec = iter::repeat(Group::Control) + .take(control.iter().len()) + .chain(iter::repeat(Group::Treatment).take(treatment.iter().len())) + .collect(); + + // do the actual permutations + (0..N_PERMUTATIONS_PER_THREAD) + .filter(|_| { + // 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 + } + }), + ); + + // 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 + } + }), + ); + + // 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 + }) + .count() as f64 + }) + .sum(); + + // p-value is the ratio of permutations where delta(mean) exceeded empircal delta(mean) to + // total permutations + let p_value = count / (N_THREADS * N_PERMUTATIONS_PER_THREAD) as f64; + + // adjust for left or right tail + if mu_diff < 0.0 { + 1.0 - p_value + } else { + p_value + } +} + +fn main() -> Result<()> { + // read data + let control = load_f64s("control.dat")?; + let treatment = load_f64s("treatment.dat")?; + + // compute empircal means + let mean_control = mean(control.iter()); + let mean_treatment = mean(treatment.iter()); + println!(" mu_control = {}", mean_control); + println!(" mu_treatment = {}", mean_treatment); + 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); + println!(" p-value = {}", pvalue); + + Ok(()) +} + +#[cfg(test)] +mod tests { + // import names from main scope + use super::*; + + const MEAN_EPSILON: f64 = 0.000001; + const PVALUE_EPSILON: f64 = 0.001; + + // for float comparison + fn approx_eq(a: f64, b: f64, epsilon: f64) -> bool { + (a - b).abs() < epsilon + } + + #[test] + fn test_mouse_data() { + // test data from Table 2.1 in "An Introduction to the Bootstrap" (book) + let control = vec![52.0, 104.0, 146.0, 10.0, 51.0, 30.0, 40.0, 27.0, 46.0]; + let treatment = vec![94.0, 197.0, 16.0, 38.0, 99.0, 141.0, 23.0]; + + // compute empircal means + let mean_control = mean(control.iter()); + assert!(approx_eq(mean_control, 56.22222222222222, MEAN_EPSILON)); + let mean_treatment = mean(treatment.iter()); + assert!(approx_eq(mean_treatment, 86.85714285714286, MEAN_EPSILON)); + assert!(approx_eq( + mean_treatment - mean_control, + 30.63492063492064, + MEAN_EPSILON + )); + + // run permutation test to compute p-value + let pvalue = permute(&control, &treatment, mean_treatment - mean_control); + assert!(approx_eq(pvalue, 0.13896357, PVALUE_EPSILON)); + } +} diff --git a/tmp b/tmp new file mode 100644 index 0000000..a5ba572 --- /dev/null +++ b/tmp @@ -0,0 +1,3 @@ +17.11 +17.1 +16.2