Add guards and determinism (sorting/emission-wise)

cpp/Averager.cpp

@@ -160,6 +160,18 @@ void Averager::compute_mean_coverage(std::vector<std::vector<BedGraphRow>>& all_
 
     std::cout << "[INFO] fastder is using " << nof_threads << " threads for computing mean coverage." << std::endl;
 
+    // Pre-build a vector of pointers to each chromosome's sample list so
+    // worker threads do read-only lookups by index. Calling chrom_samples[]
+    // from threads would be a non-const access on a shared map and could
+    // mutate it under concurrent reads.
+    std::vector<const std::vector<std::vector<BedGraphRow>>*> chrom_samples_ptrs;
+    chrom_samples_ptrs.reserve(chroms.size());
+    for (const auto& chrom : chroms)
+    {
+        const auto it = chrom_samples.find(chrom);
+        chrom_samples_ptrs.push_back(it == chrom_samples.end() ? nullptr : &it->second);
+    }
+
     for (unsigned int t = 0; t < static_cast<unsigned int>(nof_threads); ++t)
     {
         threads.emplace_back([&]()
@@ -169,7 +181,9 @@ void Averager::compute_mean_coverage(std::vector<std::vector<BedGraphRow>>& all_
                 unsigned int i = next_index++;
                 if (i >= chroms.size()) break;
                 const std::string& chrom = chroms[i];
-                std::vector<BedGraphRow> intervals = mean_for_chrom(chrom, chrom_samples[chrom], total_samples);
+                const auto* samples_ptr = chrom_samples_ptrs[i];
+                if (samples_ptr == nullptr) continue;
+                std::vector<BedGraphRow> intervals = mean_for_chrom(chrom, *samples_ptr, total_samples);
                 {
                     std::lock_guard<std::mutex> lock(map_mutex);
                     mean_intervals[chrom] = std::move(intervals);

cpp/Integrator.cpp

@@ -164,10 +164,17 @@ void Integrator::stitch_up(std::unordered_map<std::string, std::vector<BedGraphR
     //      with strand '.'.
     //   4. Sort the chromosome's StitchedERs by start so write_to_gtf reads
     //      them in genomic order even when chains came from different passes.
-    for (const auto& chrom_ers : expressed_regions)
+    //
+    // expressed_regions is an unordered_map, so its iteration order is
+    // implementation-defined. Sort the keys here to make cross-chromosome
+    // output order deterministic across runs and platforms.
+    std::vector<std::string> chroms_sorted;
+    chroms_sorted.reserve(expressed_regions.size());
+    for (const auto& chrom_ers : expressed_regions) chroms_sorted.push_back(chrom_ers.first);
+    std::sort(chroms_sorted.begin(), chroms_sorted.end());
+    for (const std::string& chrom : chroms_sorted)
     {
-        const std::string& chrom = chrom_ers.first;
-        const auto& ers = chrom_ers.second;
+        const auto& ers = expressed_regions.at(chrom);
         if (ers.empty()) continue;
 
         std::vector<uint32_t> plus_sjs;

cpp/Integrator.h

@@ -4,11 +4,8 @@
 
 #ifndef MLS_INTEGRATOR_H
 #define MLS_INTEGRATOR_H
-#include <unordered_map>
-
-#endif //MLS_INTEGRATOR_H
-
 
+#include <unordered_map>
 #include <SJRow.h>
 #include <BedGraphRow.h>
 #include <chrono>
@@ -49,4 +46,6 @@ class Integrator
     int position_tolerance = 5;
 
 
-};
+};
+
+#endif //MLS_INTEGRATOR_H

cpp/Parser.cpp

@@ -15,6 +15,7 @@
 #include <mutex>
 #include <charconv>
 #include <cstdint> // for library size which can be too large for unsigned int
+#include <cstdlib>
 
 #include "Parser.h"
 
@@ -377,11 +378,19 @@ void Parser::read_all_bedgraphs(std::vector<std::string> bedgraph_files, unsigne
                 std::vector<BedGraphRow> sample_bedgraph;
 
                 // pick the right reader by file extension. BigWig parsing is
-                // gated on libBigWig at compile time; if it is not built in,
-                // read_bigwig logs an error and returns an empty vector.
+                // gated on libBigWig at compile time; without it a .bw input
+                // is a hard error so we don't silently dilute the mean
+                // coverage by counting an empty sample toward total_samples.
                 if (filename.size() >= 3 && filename.substr(filename.size() - 3) == ".bw")
                 {
+#ifdef FASTDER_USE_LIBBIGWIG
                     sample_bedgraph = read_bigwig(filename, library_size);
+#else
+                    std::cerr << "[ERROR] BigWig input " << filename
+                              << " requires fastder built with -DFASTDER_USE_LIBBIGWIG=ON. "
+                              << "Reconfigure or convert the file to BedGraph." << std::endl;
+                    std::exit(1);
+#endif
                 }
                 else
                 {

cpp/StitchedER.h

@@ -5,8 +5,6 @@
 #ifndef MLS_STITCHEDER_H
 #define MLS_STITCHEDER_H
 
-#endif //MLS_STITCHEDER_H
-
 #include <vector>
 #include <BedGraphRow.h>
 #include <cstdint>
@@ -51,4 +49,6 @@ class StitchedER
         return os << stitched_er.across_er_coverage << "\t" << stitched_er.start << "\t" << stitched_er.end << "\t" << stitched_er.total_length << std::endl;
 
     }
-};
+};
+
+#endif //MLS_STITCHEDER_H