From 09ec3cd46d47eb40c193c1aab3e9da64ce566371 Mon Sep 17 00:00:00 2001
From: Noam <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 16 May 2023 16:57:20 -0700
Subject: [PATCH 1/2] added the direction to the fdr calculation, and perform
 relevant fold calculation with epsilon offset so that fdr will be either the
 min or max with the epsilon. also allows no nonzero results in the case where
 the direction is none. included tests as well

---
 src/fdr.rs | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 66 insertions(+), 3 deletions(-)

diff --git a/src/fdr.rs b/src/fdr.rs
index 9c0440b..9d85dd8 100644
--- a/src/fdr.rs
+++ b/src/fdr.rs
@@ -1,6 +1,8 @@
 use crate::utils::{argsort, argsort_vec, diagonal_product};
 use ndarray::{Array1, Axis};
 
+const EPSILON: f64 = 1e-10;
+
 #[derive(Clone, Copy, Debug)]
 pub enum Direction {
     Less,
@@ -82,7 +84,7 @@ impl<'a> Fdr<'a> {
         let sorted_values = values.select(Axis(0), &order);
         let sorted_ntc = is_ntc.select(Axis(0), &order);
         let sorted_fdr = Self::empirical_fdr(&sorted_ntc);
-        let threshold = Self::threshold(&sorted_values, &sorted_fdr, self.alpha);
+        let threshold = Self::threshold(&sorted_values, &sorted_fdr, self.alpha, self.use_product);
         let unsorted_fdr = sorted_fdr.select(Axis(0), &reorder);
         FdrResult::new(unsorted_fdr, threshold)
     }
@@ -112,7 +114,7 @@ impl<'a> Fdr<'a> {
     }
 
     /// Calculate the p-value threshold
-    fn threshold(values: &Array1<f64>, fdr: &Array1<f64>, alpha: f64) -> f64 {
+    fn threshold(values: &Array1<f64>, fdr: &Array1<f64>, alpha: f64, use_product: Option<Direction>) -> f64 {
         let fdr_pval = fdr
             .iter()
             .zip(values.iter())
@@ -121,13 +123,19 @@ impl<'a> Fdr<'a> {
         if let Some(fp) = fdr_pval {
             *fp.1
         } else {
-            0.0
+            match use_product {
+                Some(Direction::Greater) => values.fold(0.0, |acc: f64, x| acc.max(*x)) + EPSILON,
+                Some(Direction::Less) => values.fold(1.0, |acc: f64, x| acc.min(*x)) - EPSILON,
+                _ => (values.fold(1.0, |acc: f64, x| acc.min(*x)) - EPSILON).max(0.0),
+            }
         }
     }
 }
 
 #[cfg(test)]
 mod testing {
+    use crate::fdr::{Direction, EPSILON};
+
     use super::Fdr;
     use ndarray::{array, Array1};
 
@@ -178,4 +186,59 @@ mod testing {
         let fdr = Fdr::new(&pvalues, &logfc, &ntc_indices, alpha, None).fit();
         assert_eq!(fdr.threshold(), 1. / 3.);
     }
+
+    #[test]
+    fn test_fdr_threshold_direction_lt_saturated() {
+        let m = 10;
+        let pvalues = Array1::linspace(0.1, 1.0, m);
+        let logfc = Array1::linspace(-1.0, -0.01, m);
+        let ntc_indices = vec![0, 3, 5];
+        let alpha = 0.1;
+        let fdr = Fdr::new(&pvalues, &logfc, &ntc_indices, alpha, Some(Direction::Less)).fit();
+        assert_eq!(fdr.threshold(), -(0.1f64).log10() * -1.0 - EPSILON);
+    }
+
+    #[test]
+    fn test_fdr_threshold_direction_gt_saturated() {
+        let m = 10;
+        let pvalues = Array1::linspace(0.1, 1.0, m);
+        let logfc = Array1::linspace(0.1, 1.0, m).iter().rev().cloned().collect::<Array1<f64>>();
+        let ntc_indices = vec![0, 3, 5];
+        let alpha = 0.1;
+        let fdr = Fdr::new(&pvalues, &logfc, &ntc_indices, alpha, Some(Direction::Greater)).fit();
+        assert_eq!(fdr.threshold(), -(0.1f64).log10() * 1.0 + EPSILON);
+    }
+
+    #[test]
+    fn test_fdr_threshold_saturated() {
+        let m = 10;
+        let pvalues = Array1::linspace(0.1, 1.0, m);
+        let logfc = Array1::linspace(-1.0, -0.01, m);
+        let ntc_indices = vec![0, 3, 5];
+        let alpha = 0.1;
+        let fdr = Fdr::new(&pvalues, &logfc, &ntc_indices, alpha, None).fit();
+        assert_eq!(fdr.threshold(), 0.1f64 - EPSILON);
+    }
+
+    #[test]
+    fn test_fdr_threshold_saturated_epsilon() {
+        let m = 10;
+        let pvalues = Array1::linspace(EPSILON, 1.0, m);
+        let logfc = Array1::linspace(-1.0, -0.01, m);
+        let ntc_indices = vec![0, 3, 5];
+        let alpha = 0.1;
+        let fdr = Fdr::new(&pvalues, &logfc, &ntc_indices, alpha, None).fit();
+        assert_eq!(fdr.threshold(), 0.0);
+    }
+
+    #[test]
+    fn test_fdr_threshold_saturated_nonzero() {
+        let m = 10;
+        let pvalues = Array1::linspace(f64::EPSILON, 1.0, m);
+        let logfc = Array1::linspace(-1.0, -0.01, m);
+        let ntc_indices = vec![0, 3, 5];
+        let alpha = 0.1;
+        let fdr = Fdr::new(&pvalues, &logfc, &ntc_indices, alpha, None).fit();
+        assert_eq!(fdr.threshold(), 0.0);
+    }
 }

From 578bcecb42ca0157d48d0184c7fb03c27b04d9ef Mon Sep 17 00:00:00 2001
From: Noam <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 16 May 2023 16:57:44 -0700
Subject: [PATCH 2/2] bump version

---
 Cargo.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Cargo.toml b/Cargo.toml
index a0dbad6..2408132 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "intc"
-version = "0.2.0"
+version = "0.2.1"
 edition = "2021"
 description = "An implementation of the *-INC method to calculate an empirical FDR for non-targeting controls in CRISPR screens "
 license = "MIT"