Merge pull request #54 from pachterlab/devel

Version 4.2.2

CMakeLists.txt

@@ -4,11 +4,11 @@ project(kallisto)
 
 set(EXT_PROJECTS_DIR ${PROJECT_SOURCE_DIR}/ext)
 
-#add_compile_options(-Wall -Wno-unused-function) 
+#add_compile_options(-Wall -Wno-unused-function)
 add_compile_options(-std=c++11)
 IF(CMAKE_BUILD_TYPE MATCHES Debug)
     message("debug mode")
-		add_compile_options(-O0)
+		add_compile_options(-O0 -g)
 ELSE(CMAKE_BUILD_TYPE MATCHES Debug)
 		add_compile_options(-O3)
 ENDIF(CMAKE_BUILD_TYPE MATCHES Debug)
@@ -17,7 +17,11 @@ if(CMAKE_BUILD_TYPE MATCHES Profile)
 	add_compile_options(-g)
 endif(CMAKE_BUILD_TYPE MATCHES Profile)
 
-find_package(HDF5)
+if(LINK MATCHES static)
+    message("static build")
+ELSE(LINK MATCHES shared)
+    message("shared build")
+ENDIF(LINK MATCHES static)
 
 # add_compile_options(-Wdeprecated-register)
 

src/Bootstrap.cpp

@@ -1,15 +1,93 @@
 #include "Bootstrap.h"
-#include "weights.h"
-#include "EMAlgorithm.h"
+// #include "weights.h"
+// #include "EMAlgorithm.h"
 
-EMAlgorithm Bootstrap::run_em(const EMAlgorithm& em_start) {
+EMAlgorithm Bootstrap::run_em() {
     auto counts = mult_.sample();
-    auto weights = calc_weights(counts, ecmap_, eff_lens_);
-    EMAlgorithm em(ecmap_, counts, target_names_, eff_lens_, weights);
+    EMAlgorithm em(counts, index_, tc_, mean_fl);
 
     //em.set_start(em_start);
-    em.run(10000, 20, false);
+    em.run(10000, 50, false, false);
     /* em.compute_rho(); */
 
     return em;
 }
+
+BootstrapThreadPool::BootstrapThreadPool(
+    size_t n_threads,
+    std::vector<size_t> seeds,
+    const std::vector<int>& true_counts,
+    const KmerIndex& index,
+    const MinCollector& tc,
+    const std::vector<double>& eff_lens,
+    double mean,
+    const ProgramOptions& p_opts,
+    H5Writer& h5writer
+    ) :
+  n_threads_(n_threads),
+  seeds_(seeds),
+  n_complete_(0),
+  true_counts_(true_counts),
+  index_(index),
+  tc_(tc),
+  eff_lens_(eff_lens),
+  mean_fl_(mean),
+  opt_(p_opts),
+  writer_(h5writer)
+{
+  for (size_t i = 0; i < n_threads_; ++i) {
+    threads_.push_back( std::thread(BootstrapWorker(*this, i)) );
+  }
+}
+
+BootstrapThreadPool::~BootstrapThreadPool() {
+  for (size_t i = 0; i < n_threads_; ++i) {
+    threads_[i].join();
+  }
+}
+
+void BootstrapWorker::operator() (){
+  while (true) {
+    size_t cur_seed;
+    size_t cur_id;
+
+    // acquire a seed
+    {
+      std::unique_lock<std::mutex> lock(pool_.seeds_mutex_);
+
+      if (pool_.seeds_.empty()) {
+        // no more bootstraps to perform, this thread is done
+        return;
+      }
+
+      cur_id = pool_.seeds_.size() - 1;
+      cur_seed = pool_.seeds_.back();
+      pool_.seeds_.pop_back();
+      // std::cout << "cur seed from thread (" << thread_id_ << "): " <<
+      //   cur_seed <<  " id: " << cur_id << std::endl;
+    } // release lock
+
+    Bootstrap bs(pool_.true_counts_,
+        pool_.index_,
+        pool_.tc_,
+        pool_.eff_lens_,
+        pool_.mean_fl_,
+        cur_seed);
+
+    auto res = bs.run_em();
+
+    if (!pool_.opt_.plaintext) {
+      std::unique_lock<std::mutex> lock(pool_.write_lock_);
+      ++pool_.n_complete_;
+      std::cerr << "[bstrp] number of EM bootstraps complete: " << pool_.n_complete_ << "\r";
+      pool_.writer_.write_bootstrap(res, cur_id);
+      // release write lock
+    } else {
+      // can write out plaintext in parallel
+      plaintext_writer(pool_.opt_.output + "/bs_abundance_" +
+          std::to_string(cur_id) + ".tsv",
+          pool_.index_.target_names_, res.alpha_,
+          pool_.eff_lens_, pool_.index_.target_lens_);
+    }
+  }
+}

src/Bootstrap.h

@@ -1,41 +1,102 @@
 #ifndef KALLISTO_BOOTSTRAP_H
 #define KALLISTO_BOOTSTRAP_H
 
+#include <mutex>
+#include <thread>
+
 #include "KmerIndex.h"
+#include "MinCollector.h"
 #include "weights.h"
 #include "EMAlgorithm.h"
 #include "Multinomial.hpp"
+#include "H5Writer.h"
 
 class Bootstrap {
     // needs:
     // - "true" counts
     // - ecmap
     // - target_names
     // - eff_lens
-    public:
-        Bootstrap(const std::vector<int>& true_counts,
-                const EcMap& ecmap,
-                const std::vector<std::string>& target_names,
-                const std::vector<double>& eff_lens,
-                size_t seed) :
-            ecmap_(ecmap),
-            target_names_(target_names),
-            eff_lens_(eff_lens),
-            seed_(seed),
-            mult_(true_counts, seed_)
-        {}
-
-        // EM Algorithm generates a sample from the Multinomial, then returns
-        // an "EMAlgorithm" that has already run the EM as well as compute the
+public:
+  Bootstrap(const std::vector<int>& true_counts,
+            const KmerIndex& index,
+            const MinCollector& tc,
+            const std::vector<double>& eff_lens,
+            double mean,
+            size_t seed) :
+    index_(index),
+    tc_(tc),
+    eff_lens_(eff_lens),
+    mean_fl(mean),
+    seed_(seed),
+    mult_(true_counts, seed_)
+    {}
+
+  // EM Algorithm generates a sample from the Multinomial, then returns
+  // an "EMAlgorithm" that has already run the EM as well as compute the
         // rho values
-        EMAlgorithm run_em(const EMAlgorithm& em_start);
-
-    private:
-        const EcMap& ecmap_;
-        const std::vector<std::string>& target_names_;
-        const std::vector<double>& eff_lens_;
-        size_t seed_;
-        Multinomial mult_;
+  EMAlgorithm run_em();
+
+private:
+  const KmerIndex& index_;
+  const MinCollector& tc_;
+  const std::vector<double>& eff_lens_;
+  double mean_fl;
+  size_t seed_;
+  Multinomial mult_;
+};
+
+class BootstrapThreadPool {
+  friend class BootstrapWorker;
+
+  public:
+    BootstrapThreadPool(
+        size_t n_threads,
+        std::vector<size_t> seeds,
+        const std::vector<int>& true_counts,
+        const KmerIndex& index,
+        const MinCollector& tc,
+        const std::vector<double>& eff_lens,
+        double mean,
+        const ProgramOptions& p_opts,
+        H5Writer& h5writer
+        );
+
+    size_t num_threads() {return n_threads_;}
+
+    ~BootstrapThreadPool();
+  private:
+    std::vector<size_t> seeds_;
+    size_t n_threads_;
+
+    std::vector<std::thread> threads_;
+    std::mutex seeds_mutex_;
+    std::mutex write_lock_;
+
+    size_t n_complete_;
+
+    // things to run bootstrap
+    const std::vector<int> true_counts_;
+    const KmerIndex& index_;
+    const MinCollector& tc_;
+    const std::vector<double>& eff_lens_;
+    double mean_fl_;
+    const ProgramOptions& opt_;
+    H5Writer& writer_;
+};
+
+class BootstrapWorker {
+  public:
+    BootstrapWorker(BootstrapThreadPool& pool, size_t thread_id) :
+      pool_(pool),
+      thread_id_(thread_id)
+    {}
+
+  void operator()();
+
+  private:
+    BootstrapThreadPool& pool_;
+    size_t thread_id_;
 };
 
 #endif

src/CMakeLists.txt

@@ -4,30 +4,57 @@ file(GLOB headers *.h *.hpp)
 list(REMOVE_ITEM sources main.cpp)
 
 include_directories(../ext/)
+
 add_library(kallisto_core ${sources} ${headers})
 target_include_directories(kallisto_core PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
 
 add_executable(kallisto main.cpp)
 
+find_package( Threads REQUIRED )
+target_link_libraries(kallisto kallisto_core pthread)
+
+if(LINK MATCHES static)
+    set(BUILD_SHARED_LIBS OFF)
+    set(HDF5_USE_STATIC_LIBRARIES 1)
+
+    if (UNIX AND NOT APPLE)
+        #set(CMAKE_EXE_LINKER_FLAGS  "-static -static-libgcc -static-libstdc++")
+        set(CMAKE_EXE_LINKER_FLAGS  "-static -static-libstdc++")
+        SET(CMAKE_FIND_LIBRARY_SUFFIXES ".a")
+        set(CMAKE_EXE_LINKER_FLAGS  "-static -static-libgcc -static-libstdc++")
+    endif(UNIX AND NOT APPLE)
+
+    SET_TARGET_PROPERTIES(kallisto kallisto_core PROPERTIES LINK_SEARCH_END_STATIC 1)
+endif(LINK MATCHES static)
+
+find_package( HDF5 REQUIRED )
 find_package( ZLIB REQUIRED )
+
 if ( ZLIB_FOUND )
     include_directories( ${ZLIB_INCLUDE_DIRS} )
-    target_link_libraries( kallisto ${ZLIB_LIBRARIES} )
 else()
-    message(FATA_ERROR "zlib not found. Required for to output files" )
+    message(FATAL_ERROR "zlib not found. Required for to output files" )
 endif( ZLIB_FOUND )
 
 if(HDF5_FOUND)
     include_directories( ${HDF5_INCLUDE_DIR} )
+    target_link_libraries( kallisto_core ${HDF5_LIBRARIES} )
     target_link_libraries( kallisto ${HDF5_LIBRARIES} )
 else()
     message(FATAL_ERROR "HDF5 not found. Required to output files")
 endif()
 
-target_link_libraries(kallisto kallisto_core pthread)
 
-if (UNIX AND NOT APPLE)
-    target_link_libraries(kallisto rt)
-endif()
+if(LINK MATCHES static)
+    if (UNIX AND NOT APPLE)
+        target_link_libraries(kallisto librt.a)
+    endif()
+else()
+    if (UNIX AND NOT APPLE)
+        target_link_libraries(kallisto rt)
+    endif()
+endif(LINK MATCHES static)
+
+
 
 install(TARGETS kallisto DESTINATION bin)