Correct entropy regional medians and singleton summaries

modkit-core/src/entropy/mod.rs

@@ -1386,8 +1386,13 @@ impl DescriptiveStats {
                 "measurements and n_reads should be the same length"
             );
             let mean_entropy = Self::mean(measurements);
-            let median_entropy =
-                percentile_linear_interp(measurements, 0.5f32)?;
+            let median_entropy = if measurements.len() == 1 {
+                measurements[0]
+            } else {
+                let mut sorted_measurements = measurements.to_vec();
+                sorted_measurements.sort_unstable_by(f32::total_cmp);
+                percentile_linear_interp(&sorted_measurements, 0.5f32)?
+            };
             // safe because of above check
             let (min_entropy, max_entropy) = match measurements.iter().minmax()
             {
@@ -1671,7 +1676,118 @@ impl BedRegion {
 
 #[cfg(test)]
 mod entropy_mod_tests {
-    use crate::entropy::BedRegion;
+    use crate::entropy::{BedRegion, DescriptiveStats};
+    use itertools::Itertools;
+
+    #[test]
+    fn singleton_entropy_summary_uses_its_value_as_the_median() {
+        let stats = DescriptiveStats::new(&[0.25], &[7], 2, 0, &(10..11))
+            .expect("a singleton is a valid regional entropy summary");
+
+        assert_eq!(stats.mean_entropy, 0.25);
+        assert_eq!(stats.median_entropy, 0.25);
+        assert_eq!(stats.min_entropy, 0.25);
+        assert_eq!(stats.max_entropy, 0.25);
+        assert_eq!(stats.mean_num_reads, 7.0);
+        assert_eq!(stats.min_num_reads, 7);
+        assert_eq!(stats.max_num_reads, 7);
+        assert_eq!(stats.successful_count, 1);
+        assert_eq!(stats.failed_count, 2);
+    }
+
+    #[test]
+    fn multi_window_entropy_summary_keeps_existing_exact_statistics() {
+        let stats = DescriptiveStats::new(
+            &[0.25, 0.5, 0.75],
+            &[2, 4, 6],
+            1,
+            0,
+            &(10..20),
+        )
+        .unwrap();
+
+        assert_eq!(stats.mean_entropy, 0.5);
+        assert_eq!(stats.median_entropy, 0.5);
+        assert_eq!(stats.min_entropy, 0.25);
+        assert_eq!(stats.max_entropy, 0.75);
+        assert_eq!(stats.mean_num_reads, 4.0);
+        assert_eq!(stats.min_num_reads, 2);
+        assert_eq!(stats.max_num_reads, 6);
+        assert_eq!(stats.successful_count, 3);
+        assert_eq!(stats.failed_count, 1);
+    }
+
+    #[test]
+    fn region_median_is_independent_of_window_encounter_order() {
+        let observed = DescriptiveStats::new(
+            &[0.9, 0.1, 0.4, 0.2],
+            &[9, 1, 4, 2],
+            0,
+            0,
+            &(10..20),
+        )
+        .unwrap();
+        let sorted = DescriptiveStats::new(
+            &[0.1, 0.2, 0.4, 0.9],
+            &[1, 2, 4, 9],
+            0,
+            0,
+            &(10..20),
+        )
+        .unwrap();
+
+        assert_eq!(observed.median_entropy, sorted.median_entropy);
+        assert_eq!(observed.median_entropy, 0.3);
+    }
+
+    #[test]
+    fn odd_region_median_is_independent_of_window_encounter_order() {
+        let observed = DescriptiveStats::new(
+            &[0.9, 0.1, 0.4],
+            &[9, 1, 4],
+            0,
+            0,
+            &(10..20),
+        )
+        .unwrap();
+        let sorted = DescriptiveStats::new(
+            &[0.1, 0.4, 0.9],
+            &[1, 4, 9],
+            0,
+            0,
+            &(10..20),
+        )
+        .unwrap();
+
+        assert_eq!(observed.median_entropy, sorted.median_entropy);
+        assert_eq!(observed.median_entropy, 0.4);
+    }
+
+    #[test]
+    fn regional_medians_are_stable_across_all_small_encounter_orders() {
+        for (values, expected_median) in [
+            (vec![0.75, 0.25, 0.5], 0.5),
+            (vec![0.75, 0.0, 0.25, 0.5], 0.375),
+        ] {
+            for measurements in
+                values.iter().copied().permutations(values.len())
+            {
+                let original_order = measurements.clone();
+                let reads = vec![1; measurements.len()];
+                let stats = DescriptiveStats::new(
+                    &measurements,
+                    &reads,
+                    0,
+                    0,
+                    &(10..20),
+                )
+                .unwrap();
+
+                assert_eq!(stats.median_entropy, expected_median);
+                assert_eq!(measurements, original_order);
+            }
+        }
+    }
 
     #[test]
     fn test_bed_region_parsing() {

modkit/tests/test_entropy_region_statistics.rs

@@ -0,0 +1,106 @@
+use rust_htslib::bam::{
+    self,
+    header::HeaderRecord,
+    record::{Aux, Cigar, CigarString},
+};
+use std::fs;
+use std::path::{Path, PathBuf};
+use std::process::Command;
+
+fn write_reference(root: &Path) -> PathBuf {
+    let reference = root.join("reference.fa");
+    fs::write(&reference, ">chr1\nAACGAA\n").unwrap();
+    fs::write(root.join("reference.fa.fai"), "chr1\t6\t6\t6\t7\n").unwrap();
+    reference
+}
+
+fn write_bam(root: &Path) -> PathBuf {
+    let bam_path = root.join("reads.bam");
+    let mut header = bam::Header::new();
+    let mut sq = HeaderRecord::new(b"SQ");
+    sq.push_tag(b"SN", "chr1").push_tag(b"LN", 6);
+    header.push_record(&sq);
+
+    let cigar = CigarString(vec![Cigar::Match(6)]);
+    let mut record = bam::Record::new();
+    record.set(b"read-0", Some(&cigar), b"AACGAA", &[30; 6]);
+    record.set_tid(0);
+    record.set_pos(0);
+    record.set_mapq(60);
+    record.push_aux(b"MM", Aux::String("C+m?,0;")).unwrap();
+    record.push_aux(b"ML", Aux::ArrayU8((&[255][..]).into())).unwrap();
+    record.push_aux(b"MN", Aux::U32(6)).unwrap();
+    record.push_aux(b"NM", Aux::U32(0)).unwrap();
+
+    let mut writer =
+        bam::Writer::from_path(&bam_path, &header, bam::Format::Bam).unwrap();
+    writer.write(&record).unwrap();
+    drop(writer);
+    bam::index::build(&bam_path, None, bam::index::Type::Bai, 1).unwrap();
+    bam_path
+}
+
+#[test]
+fn singleton_region_emits_exact_summary_statistics() {
+    let temp_dir = tempfile::tempdir().unwrap();
+    let reference = write_reference(temp_dir.path());
+    let bam = write_bam(temp_dir.path());
+    let regions = temp_dir.path().join("regions.bed");
+    fs::write(&regions, "chr1\t2\t3\tsingleton\n").unwrap();
+    let output_dir = temp_dir.path().join("output");
+
+    let result = Command::new(env!("CARGO_BIN_EXE_modkit"))
+        .args([
+            "entropy",
+            "--in-bam",
+            bam.to_str().unwrap(),
+            "--out-bed",
+            output_dir.to_str().unwrap(),
+            "--ref",
+            reference.to_str().unwrap(),
+            "--base",
+            "C",
+            "--num-positions",
+            "1",
+            "--window-size",
+            "1",
+            "--min-coverage",
+            "1",
+            "--max-filtered-positions",
+            "0",
+            "--filter-threshold",
+            "0",
+            "--threads",
+            "1",
+            "--io-threads",
+            "1",
+            "--regions",
+            regions.to_str().unwrap(),
+            "--prefix",
+            "singleton",
+            "--suppress-progress",
+        ])
+        .output()
+        .unwrap();
+
+    assert!(
+        result.status.success(),
+        "{}",
+        String::from_utf8_lossy(&result.stderr)
+    );
+    let windows = fs::read_to_string(
+        output_dir.join("singleton_windows.bedgraph"),
+    )
+    .unwrap();
+    let window_fields = windows.trim_end().split('\t').collect::<Vec<_>>();
+    // Motif/window interval normalization belongs to issue #681. This test
+    // isolates the regional summary and only requires one successful window
+    // with the expected biological value and coverage.
+    assert_eq!(window_fields[0], "chr1");
+    assert_eq!(window_fields[1], "2");
+    assert_eq!(&window_fields[3..], &["0", "+", "1"]);
+    assert_eq!(
+        fs::read(output_dir.join("singleton_regions.bed")).unwrap(),
+        b"chr1\t2\t3\tsingleton\t0\t+\t0\t0\t0\t1\t1\t1\t1\t0\n"
+    );
+}