solo: EmptyDrops_CR on CellRanger's actual statistics (SGT ambient profile, libc++ sampler)

CHANGELOG.md

@@ -97,6 +97,16 @@ Sections commonly used: Features, Bug fixes, Other changes.
   glibc's malloc and per-thread heaps that return whole segments to
   the OS when abandoned, so allocator cache size stays bounded.
 
+- **`--soloCellFilter EmptyDrops_CR` now uses CellRanger's actual
+  statistics.** The ambient profile is smoothed with Simple Good-Turing,
+  as CellRanger and STAR do, instead of an approximation that reserved
+  unseen mass from the singleton rate and spread the remainder in
+  proportion to raw counts. The Monte-Carlo null is drawn with libc++'s
+  `std::mt19937` and `std::discrete_distribution`, seeded
+  `19760110 * (isim + 1)` per simulation as STAR seeds it, replacing a
+  SplitMix64 stream that could not agree with STAR's over an arbitrary
+  number of draws. Cell calls move as a result.
+
 ### Bug fixes
 
 - **STARsolo `Gene` assignment now requires exon concordance**, matching

DIVERGENCE.md

@@ -79,6 +79,18 @@ For `--quantMode TranscriptomeSAM`, rustar-aligner builds the per-transcript exo
 
 rustar-aligner uses an in-tree splitmix64 (`src/rng.rs`) rather than the `rand` crate, avoiding the `getrandom`/`zerocopy`/`ppv-lite86` dependency chain. This is the generator underlying §1.1; it is called out separately because it is a dependency/implementation choice independent of the tie-break policy.
 
+### 1.2 `EmptyDrops_CR` Simple-Good-Turing with fewer than five distinct frequencies
+
+**What STAR does.** The ambient profile for `--soloCellFilter EmptyDrops_CR` is smoothed with Simple Good-Turing (Elworthy's `SimpleGoodTuring/sgt.h`). `analyse()` returns early, doing nothing, when the frequency spectrum has fewer than five distinct counts — Elworthy's `MinInput` guard. `PZero`, the mass reserved for genes unseen in the ambient droplets, is neither assigned in that case nor initialised at construction, so a caller that asks for it reads whatever the stack held.
+
+**What rustar-aligner does.** `PZero` is zero from construction.
+
+**Why.** There is nothing to reproduce: the value STAR reads is not a decision it made. With fewer than five distinct frequencies there is no basis for reserving unseen mass, and zero says so. Reproducing STAR would mean writing code whose correct behaviour is to emit an uninitialised value, and a test asserting it.
+
+**Impact.** Degenerate inputs only — any dataset large enough to reach the significance test has far more than five distinct frequencies. On those inputs an uninitialised read can place arbitrary mass on unseen genes, which makes the multinomial log-probabilities meaningless; zero keeps them defined.
+
+**Source.** `src/solo/sgt.rs`, locked by `solo::sgt::tests::too_few_frequencies_leaves_the_unseen_mass_at_zero` (asserting the exact bit pattern, since the point is that nothing was written). STAR: `SoloFeature_emptyDrops_CR.cpp`, `SimpleGoodTuring/sgt.h`.
+
 ---
 
 ## 5. Known residual single-read differences

src/solo/count.rs

@@ -959,20 +959,60 @@ fn emptydrops_called(
         return Ok(called);
     }
 
-    // Ambient probabilities with a Good-Turing P0 unseen-mass correction.
-    let n1 = ambient.iter().filter(|&&x| (x - 1.0).abs() < 0.5).count() as f64;
-    let p0 = (n1 / amb_total).clamp(1e-12, 0.5);
-    let n_zero = ambient.iter().filter(|&&x| x == 0.0).count().max(1) as f64;
-    let amb_p: Vec<f64> = ambient
-        .iter()
-        .map(|&x| {
-            if x > 0.0 {
-                (1.0 - p0) * x / amb_total
-            } else {
-                p0 / n_zero
+    // Ambient probabilities, smoothed by Simple Good-Turing as CellRanger's
+    // EmptyDrops_CR does. The raw ambient counts are a small sample, so a gene
+    // seen twice there is not twice as likely as one seen once; SGT fits the
+    // frequency spectrum and reserves mass for genes the sample missed
+    // entirely. The previous approximation here reserved mass from the
+    // singleton rate alone and spread the rest in proportion to raw counts,
+    // which is the right shape but the wrong numbers.
+    let amb_p: Vec<f64> = {
+        // Frequency of frequencies over the integer ambient counts.
+        let counts: Vec<u32> = ambient.iter().map(|&x| x as u32).collect();
+        let mut fof: std::collections::BTreeMap<u32, u32> = std::collections::BTreeMap::new();
+        for &c in &counts {
+            *fof.entry(c).or_insert(0) += 1;
+        }
+        let n_zero = fof.get(&0).copied().unwrap_or(0);
+
+        let mut sgt = crate::solo::sgt::Sgt::new();
+        for (&freq, &n) in &fof {
+            if freq != 0 {
+                sgt.add(freq, n);
             }
-        })
-        .collect();
+        }
+        let fitted = sgt.analyse();
+
+        // Per-gene probability. A gene absent from the ambient set shares the
+        // reserved unseen mass; one present takes its smoothed estimate. When
+        // the spectrum is too small to fit (D17), there is no reserved mass and
+        // the raw proportions are all that is left.
+        let unseen_each = if n_zero > 0 {
+            sgt.pzero() / f64::from(n_zero)
+        } else {
+            0.0
+        };
+        let raw: Vec<f64> = counts
+            .iter()
+            .map(|&c| {
+                if c == 0 {
+                    unseen_each
+                } else if fitted {
+                    sgt.estimate(c).unwrap_or(f64::from(c) / amb_total)
+                } else {
+                    f64::from(c) / amb_total
+                }
+            })
+            .collect();
+        // Renormalise: the smoothed estimates are per-event probabilities and
+        // only sum to one over the whole spectrum, not over this gene set.
+        let total: f64 = raw.iter().sum();
+        if total > 0.0 {
+            raw.into_iter().map(|p| p / total).collect()
+        } else {
+            raw
+        }
+    };
     let amb_logp: Vec<f64> = amb_p.iter().map(|&p| p.max(1e-300).ln()).collect();
 
     // Observed multinomial log-prob per candidate.
@@ -1004,32 +1044,32 @@ fn emptydrops_called(
     // log-prob at each count; compare each candidate against sim[*][its total].
     let nonzero: Vec<usize> = (0..n_features).filter(|&g| amb_p[g] > 0.0).collect();
     let weights: Vec<f64> = nonzero.iter().map(|&g| amb_p[g]).collect();
-    // Ambient categorical sampler: cumulative weights + splitmix64 (crate::rng).
-    // WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls.
-    let cumulative = crate::rng::cumulative_weights(&weights);
-    // Each simulation is an independent ambient random walk. Seed a dedicated RNG
-    // per simulation (splitmix-derived from the base seed) so the result is
-    // deterministic regardless of how the work is scheduled across threads, then
-    // run the simulations in parallel. Each walk records the running log-prob at
-    // every count level; `walks[s][k]` is the log-prob of simulation `s` after `k`
-    // draws. (This matches STAR's per-thread-RNG approach; the per-sim seeding
-    // gives different draws than a single sequential stream, but the same
-    // distribution — p-values are stable to Monte-Carlo error.)
+    // The ambient sampler is STAR's, which is libc++'s: `std::mt19937` drawn
+    // through `std::discrete_distribution`. Both are implementation-defined in
+    // the parts that decide which category a draw lands in, so "a correct
+    // categorical sampler" is not enough to reproduce STAR's cell calls —
+    // it has to be that one. See `solo::libcxx_rng`.
+    let dist = crate::solo::libcxx_rng::DiscreteDistribution::new(&weights);
+    // A fresh generator per simulation, seeded `19760110 * (isim + 1)` as STAR
+    // seeds it. No shared state, so the walks can run in any order and on any
+    // number of threads and still produce the same p-values. Each walk records
+    // the running log-prob at every count level; `walks[s][k]` is simulation
+    // `s` after `k` draws.
     use rayon::prelude::*;
     const BASE_SEED: u64 = 19_760_110;
     let walks: Vec<Vec<f64>> = (0..sim_n)
         .into_par_iter()
         .map_init(
             || (vec![0u32; n_features], Vec::<usize>::new()),
             |(curr, touched), s| {
-                let seed = BASE_SEED ^ (s as u64).wrapping_mul(0x9E37_79B9_7F4A_7C15);
-                let mut rng = crate::rng::SplitMix64::seed(seed);
+                let seed = BASE_SEED.wrapping_mul(s as u64 + 1) as u32;
+                let mut rng = crate::solo::libcxx_rng::Mt19937::new(seed);
                 touched.clear();
                 let mut walk = Vec::with_capacity(max_count + 1);
                 walk.push(0.0);
                 let mut lp = 0f64;
                 for ic in 1..=max_count {
-                    let gi = nonzero[crate::rng::sample_cumulative(&cumulative, &mut rng)];
+                    let gi = nonzero[dist.sample(&mut rng)];
                     if curr[gi] == 0 {
                         touched.push(gi);
                     }

src/solo/libcxx_rng.rs

@@ -0,0 +1,219 @@
+//! Bit-exact ports of the three libc++ random facilities STARsolo's
+//! `EmptyDrops_CR` depends on.
+//!
+//! STAR's Monte-Carlo rescue draws from `std::mt19937`, converts to doubles
+//! with `std::generate_canonical<double, 53>`, and samples categories with
+//! `std::discrete_distribution`. All three are implementation-defined in the
+//! parts that matter: the standard fixes `mt19937`'s output but not how
+//! `generate_canonical` consumes it, and says nothing about how
+//! `discrete_distribution` maps a uniform draw onto categories. So "port the
+//! algorithm" is not enough — it has to be *libc++'s* algorithm, because that
+//! is what STAR was built against and what its numbers come out of.
+//!
+//! Every value asserted in the tests below was produced by compiling a C++
+//! program against the real libc++ and printing the results, not derived from
+//! reading the source. That is the only way to be sure.
+//!
+//! `solo::count` currently samples with a `SplitMix64` stream, under a comment
+//! calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops
+//! cell calls". That claim cannot hold in general: two unrelated generators
+//! cannot agree on an arbitrary number of draws, so it is true of the cases
+//! that happened to be checked and unknown everywhere else. These types remove
+//! the guesswork. Wiring them into the EmptyDrops path is the next step and is
+//! deliberately separate, since it moves cell calls.
+
+/// libc++'s `std::mt19937`.
+///
+/// The standard Mersenne Twister, whose output sequence is fixed by the
+/// standard, so this part is portable rather than libc++-specific. It is here
+/// because the two facilities that follow are not.
+#[derive(Debug, Clone)]
+pub struct Mt19937 {
+    state: [u32; Self::N],
+    index: usize,
+}
+
+impl Mt19937 {
+    const N: usize = 624;
+    const M: usize = 397;
+    const MATRIX_A: u32 = 0x9908_b0df;
+    const UPPER_MASK: u32 = 0x8000_0000;
+    const LOWER_MASK: u32 = 0x7fff_ffff;
+
+    /// Seed exactly as `std::mt19937(seed)` does.
+    pub fn new(seed: u32) -> Self {
+        let mut state = [0u32; Self::N];
+        state[0] = seed;
+        for i in 1..Self::N {
+            state[i] = 1_812_433_253u32
+                .wrapping_mul(state[i - 1] ^ (state[i - 1] >> 30))
+                .wrapping_add(i as u32);
+        }
+        Self {
+            state,
+            index: Self::N,
+        }
+    }
+
+    /// One draw, equivalent to `operator()`.
+    pub fn next_u32(&mut self) -> u32 {
+        if self.index >= Self::N {
+            self.twist();
+        }
+        let mut y = self.state[self.index];
+        self.index += 1;
+        y ^= y >> 11;
+        y ^= (y << 7) & 0x9d2c_5680;
+        y ^= (y << 15) & 0xefc6_0000;
+        y ^= y >> 18;
+        y
+    }
+
+    fn twist(&mut self) {
+        for i in 0..Self::N {
+            let y = (self.state[i] & Self::UPPER_MASK)
+                | (self.state[(i + 1) % Self::N] & Self::LOWER_MASK);
+            let mut next = self.state[(i + Self::M) % Self::N] ^ (y >> 1);
+            if y & 1 != 0 {
+                next ^= Self::MATRIX_A;
+            }
+            self.state[i] = next;
+        }
+        self.index = 0;
+    }
+
+    /// libc++'s `std::generate_canonical<double, 53>`.
+    ///
+    /// This is where implementations diverge. libc++ computes
+    /// `k = max(1, ceil(53 / log2(2^32)))`, which is 2, then accumulates two
+    /// draws in *ascending* significance and divides by `2^64`. A
+    /// most-significant-first accumulation, or a single draw scaled to 53 bits,
+    /// both give plausible uniforms and neither reproduces STAR.
+    pub fn canonical_f64(&mut self) -> f64 {
+        // r = 2^32, k = 2, so the base is r^k = 2^64.
+        let base = 2f64.powi(64);
+        let mut sum = 0f64;
+        let mut factor = 1f64;
+        for _ in 0..2 {
+            sum += f64::from(self.next_u32()) * factor;
+            factor *= 2f64.powi(32);
+        }
+        sum / base
+    }
+}
+
+/// libc++'s `std::discrete_distribution`.
+///
+/// Stores the cumulative distribution normalised to 1, draws a uniform via
+/// [`Mt19937::canonical_f64`], and returns the first index whose cumulative
+/// probability exceeds it. The final category is the fallback, so a draw of
+/// exactly 1.0 cannot fall off the end.
+#[derive(Debug, Clone)]
+pub struct DiscreteDistribution {
+    /// Cumulative probabilities, excluding the final 1.0.
+    cumulative: Vec<f64>,
+}
+
+impl DiscreteDistribution {
+    /// Build from unnormalised weights, as `discrete_distribution(w)` does.
+    pub fn new(weights: &[f64]) -> Self {
+        let total: f64 = weights.iter().sum();
+        let mut cumulative = Vec::with_capacity(weights.len().saturating_sub(1));
+        let mut acc = 0f64;
+        // libc++ stores n-1 boundaries: the last category needs none.
+        for w in weights.iter().take(weights.len().saturating_sub(1)) {
+            acc += w;
+            cumulative.push(if total > 0.0 { acc / total } else { 0.0 });
+        }
+        Self { cumulative }
+    }
+
+    /// One sample.
+    pub fn sample(&self, rng: &mut Mt19937) -> usize {
+        let u = rng.canonical_f64();
+        self.cumulative.partition_point(|&c| c <= u)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    // Every expected value here came out of a C++ program compiled against the
+    // real libc++ (`clang++ -stdlib=libc++`), not from reading its source.
+
+    #[test]
+    fn mt19937_matches_libcxx_stream() {
+        let mut g = Mt19937::new(19_760_110);
+        let got: Vec<u32> = (0..10).map(|_| g.next_u32()).collect();
+        assert_eq!(
+            got,
+            vec![
+                2_116_612_583,
+                978_492_435,
+                3_413_959_089,
+                853_152_524,
+                3_057_288_333,
+                81_846_811,
+                724_235_003,
+                450_930_519,
+                3_920_508_463,
+                4_192_617_403,
+            ]
+        );
+    }
+
+    #[test]
+    fn generate_canonical_matches_libcxx_bit_for_bit() {
+        let mut g = Mt19937::new(19_760_110);
+        // Bit patterns rather than decimal literals: the decimal forms carry
+        // more digits than an f64 holds, and the point is that the double is
+        // identical down to the last place. A difference there changes which
+        // category a sample lands in.
+        let expected: [u64; 5] = [
+            0x3fcd_294e_09bf_1479, // 0.22782302356623399
+            0x3fc9_6d09_8665_be71, // 0.19864005148291455
+            0x3f93_8388_6ed8_ea12, // 0.019056445851882549
+            0x3fba_e0a7_572b_2af3, // 0.10499044302120435
+            0x3fef_3cc8_777d_35c7, // 0.97616980874743653
+        ];
+        for (i, &want) in expected.iter().enumerate() {
+            let got = g.canonical_f64();
+            assert_eq!(
+                got.to_bits(),
+                want,
+                "draw {i}: got {got:.17}, libc++ gives {:.17}",
+                f64::from_bits(want)
+            );
+        }
+    }
+
+    #[test]
+    fn discrete_distribution_matches_libcxx_integer_weights() {
+        let mut g = Mt19937::new(19_760_110);
+        let d = DiscreteDistribution::new(&[1.0, 2.0, 3.0, 4.0]);
+        let got: Vec<usize> = (0..20).map(|_| d.sample(&mut g)).collect();
+        assert_eq!(
+            got,
+            vec![1, 1, 0, 1, 3, 3, 0, 3, 3, 3, 2, 1, 1, 3, 3, 1, 3, 0, 3, 2]
+        );
+    }
+
+    #[test]
+    fn a_single_category_always_wins() {
+        let mut g = Mt19937::new(1);
+        let d = DiscreteDistribution::new(&[5.0]);
+        for _ in 0..8 {
+            assert_eq!(d.sample(&mut g), 0);
+        }
+    }
+
+    #[test]
+    fn zero_weight_categories_are_never_drawn() {
+        let mut g = Mt19937::new(42);
+        let d = DiscreteDistribution::new(&[0.0, 1.0, 0.0]);
+        for _ in 0..64 {
+            assert_eq!(d.sample(&mut g), 1);
+        }
+    }
+}

src/solo/mod.rs

@@ -11,6 +11,8 @@
 
 pub mod count;
 pub mod gene;
+pub mod libcxx_rng;
+pub mod sgt;
 pub mod smartseq;
 pub mod whitelist;