.gitignore

@@ -87,3 +87,9 @@ unit_tests
 stuff
 .dirstamp
 configure.lineno
+
+# test coverage
+*.gcda
+*.gcno
+libbash_test.info
+testCoverage/

.travis.yml

@@ -13,14 +13,14 @@ matrix:
     - os: osx
       env: CXXFLAGS="-m32 -Werror"
 sudo: required
-dist: trusty
+dist: xenial
 before_install: .travis/setup.sh "$TRAVIS_OS_NAME" "$CXX" "$CXXFLAGS"
 before_script:
   - ./bootstrap
 script:
   - make unittests
-  - make scrm scrm_dbg && ./tests/test_binaries.sh
   - make algorithmtest
+  - ./tests/test_binaries.sh "$TRAVIS_OS_NAME"
 before_deploy:
   - ./.travis/build_src_pkg.sh
   - ./.travis/build_static_binaries.sh
@@ -40,3 +40,4 @@ deploy:
     condition: "$CC = gcc"
     condition: "$CXXFLAGS == -Werror"
     condition: "$TRAVIS_OS_NAME == linux"
+

.travis/test_coverage.sh

+#!/bin/bash
+
+os=$1
+cxx=$2
+cxxflags=$3
+
+echo "os: $os; cxx: $cxx; cxxflags: $cxxflags"
+
+if [ "$os" == "linux" ]; then
+  if [[ "$cxxflags" == *"-m32"* ]]; then
+    echo "On 32bit. Skipping coverage test."
+  else
+    coveralls --exclude lib --exclude tests --gcov-options '\-lp'
+  fi
+fi
+
+if [ "$os" == "osx" ]; then
+    echo "On OS X. Skipping coverage test."
+fi
+

Makefile.am

@@ -3,7 +3,7 @@ bin_PROGRAMS = scrm
 man_MANS = doc/scrm.1
 
 TESTS = unit_tests algorithm_tests
-check_PROGRAMS = unit_tests algorithm_tests scrm_dbg scrm_prof
+check_PROGRAMS = unit_tests algorithm_tests scrm_dbg scrm_asan scrm_prof
 PROG = SCRM
 
 dist-hook:
@@ -59,13 +59,15 @@ alg_test_src = tests/cppunit/test_runner.cc tests/algorithmtest/test_algorithm.c
 scrm_SOURCES = $(scrm_src) src/scrm.cc
 scrm_dbg_SOURCES = $(scrm_src) $(debug_src) src/scrm.cc
 scrm_prof_SOURCES = $(scrm_src) src/scrm.cc
+scrm_asan_SOURCES = $(scrm_src) src/scrm.cc
 unit_tests_SOURCES = $(scrm_src) $(debug_src) $(unit_test_src)
 algorithm_tests_SOURCES = $(scrm_src) $(alg_test_src)
 
 scrm_CXXFLAGS= -DNDEBUG @OPT_CXXFLAGS@
 scrm_dbg_CXXFLAGS= -g
 scrm_prof_CXXFLAGS= -pg -DNDEBUG
-unit_tests_CXXFLAGS = -g -DUNITTEST -DNDEBUG
+scrm_asan_CXXFLAGS= -g -DNDEBUG -fsanitize=undefined,address -fno-sanitize-recover
+unit_tests_CXXFLAGS = -g -DUNITTEST -DNDEBUG @TEST_CXXFLAGS@ 
 algorithm_tests_CXXFLAGS = -g -DNDEBUG
 unit_tests_LDADD= -L/opt/local/lib -lcppunit -ldl #link the cppunit unittest library in mac, cppunit was installed via macports
 algorithm_tests_LDADD= -L/opt/local/lib -lcppunit -ldl  #link the cppunit unittest library in mac, cppunit was installed via macports

NEWS.md

 scrm Version History
 ========================
 
+scrm 1.7.3
+------------------------
+Released: 2018-11-18
+
+### Improvements
++ The labeling of internal nodes in the oriented forest output
+  no longer changes within the same tree topology. That makes
+  it easier to identify different topologies.
+
+### Bug Fixes
++ Now requires that subpopulation are defined via the `-I` argument
+  before any demographic options for all population (`-M`, `-G` or `-N`)
+  are given (#108).
+  Thanks to Jonathan Terhorst (@terhorst) for reporting this bug.
++ When the multiple of the merge/admixure arguemtns (`-ep`, `-eps` and `-ej`)
+  are used at the time, they are now executed in order in which
+  they are provided on the command line (#121). This behavior now
+  follows ms' implementation. 
+  Thanks to Scott T Small (@stsmall) for reporting this.
+
+
+
 scrm 1.7.2
 ------------------------
 Released: 2016-04-10

configure.ac

-AC_INIT([scrm], [1.7.2],[https://github.com/paulstaab/scrm/issues])
+AC_INIT([scrm], [1.7.3],[https://github.com/paulstaab/scrm/issues])
 AM_INIT_AUTOMAKE([subdir-objects -Wall -Werror foreign])
 
 # Use -O3 as default optimization level
@@ -33,6 +33,9 @@ AC_C_CONST
 AC_C_INLINE
 AC_TYPE_SIZE_T
 
+# Check if we reporting test coverage
+AC_SUBST(TEST_CXXFLAGS)
+
 # Enable Link-time optimization if supported (gcc only)
 if test x$CXX = xg++; then
   AX_CHECK_COMPILE_FLAG([-flto], [OPT_CXXFLAGS="$OPT_CXXFLAGS -flto"], [], [-Werror])

debian/changelog

+scrm (1.7.3-1) unstable; urgency=medium
+
+  * Team upload.
+  * New upstream version
+  * Standards-Version: 4.2.1
+
+ -- Andreas Tille <tille@debian.org>  Mon, 17 Dec 2018 13:26:38 +0100
+
 scrm (1.7.2-2) unstable; urgency=medium
 
   * Team upload.

debian/control

@@ -5,7 +5,7 @@ Section: science
 Priority: optional
 Build-Depends: debhelper (>= 11~),
                libcppunit-dev
-Standards-Version: 4.1.5
+Standards-Version: 4.2.1
 Vcs-Browser: https://salsa.debian.org/med-team/scrm
 Vcs-Git: https://salsa.debian.org/med-team/scrm.git
 Homepage: https://github.com/scrm/scrm

src/forest-debug.cc

@@ -317,7 +317,7 @@ bool Forest::checkTree(Node const* root) const {
     } 
   }
 
-  return child1*child2;
+  return child1 && child2;
 }
 
 

src/forest.cc

@@ -1103,39 +1103,29 @@ void Forest::implementMigration(const Event &event, const bool &recalculate, Tim
 void Forest::implementFixedTimeEvent(TimeIntervalIterator &ti) {
   dout << "* * Fixed time event" << std::endl;
   double sample;
-  bool migrated;
-  size_t chain_cnt, pop_number = model().population_number();
+  auto mig_events = model().single_mig_events();
 
   for (size_t i = 0; i < 2; ++i) {
     if (states_[i] != 1) continue;
-    chain_cnt = 0;
-    while (true) {
-      migrated = false;
+    
     sample = random_generator()->sample();
+    for (auto me : mig_events) {
+      if (active_node(i)->population() == me.source_pop) {
+        sample -= me.prob;
+      }
 
-      for (size_t j = 0; j < pop_number; ++j) {
-        sample -= model().single_mig_pop(active_node(i)->population(), j);
       if (sample < 0) {
+        dout << "* * * a" << i << ": Migration from " 
+             << me.source_pop << " to " 
+             << me.sink_pop << std::endl;
         tmp_event_ = Event((*ti).start_height());
-          tmp_event_.setToMigration(active_node(i), i, j);
+        tmp_event_.setToMigration(active_node(i), i, me.sink_pop);
         implementMigration(tmp_event_, false, ti);
-          migrated = true;
-          break;
+        sample = random_generator()->sample();
       }
     }
-
-      // Stop if no migration occurred
-      if (!migrated) {
-        dout << "* * No fixed time migration occurred" << std::endl;
-        break;
   }
 
-      // Resolve a maximum of 10k chained events for each node
-      if (chain_cnt == 10000)
-        throw std::logic_error("Circle detected when moving individuals between populations");
-      ++chain_cnt;
-    }
-  }
   assert( printTree() );
 }
 

src/model.cc

@@ -95,7 +95,7 @@ size_t Model::addChangeTime(double time, const bool &scaled) {
     growth_rates_list_.push_back(std::vector<double>());
     mig_rates_list_.push_back(std::vector<double>());
     total_mig_rates_list_.push_back(std::vector<double>());
-    single_mig_probs_list_.push_back(std::vector<double>());
+    single_mig_list_.push_back(std::vector<MigEvent>());
     return position;
   }
 
@@ -113,7 +113,7 @@ size_t Model::addChangeTime(double time, const bool &scaled) {
   growth_rates_list_.insert(growth_rates_list_.begin() + position, std::vector<double>());
   mig_rates_list_.insert(mig_rates_list_.begin() + position, std::vector<double>());
   total_mig_rates_list_.insert(total_mig_rates_list_.begin() + position, std::vector<double>());
-  single_mig_probs_list_.insert(single_mig_probs_list_.begin() + position, std::vector<double>());
+  single_mig_list_.insert(single_mig_list_.begin() + position, std::vector<MigEvent>());
   return position;
 }
 
@@ -433,18 +433,19 @@ void Model::addSingleMigrationEvent(const double time, const size_t source_pop,
                                     const bool &time_scaled) {
 
   size_t position = addChangeTime(time, time_scaled);
-  size_t popnr = population_number();
 
   if ( time < 0.0 ) throw std::invalid_argument("Single migration event: Negative time");
   if ( source_pop >= population_number() ) throw std::invalid_argument("Single migration event: Unknown population");
   if ( sink_pop >= population_number() ) throw std::invalid_argument("Single migration event: Unknown population");
   if ( fraction < 0.0 || fraction > 1.0 ) throw std::invalid_argument("Single migration event: Fraction out of range");
 
-  if ( single_mig_probs_list_.at(position).empty() ) {
-    single_mig_probs_list_.at(position) = std::vector<double>(popnr*popnr-popnr, 0.0);
+  if ( single_mig_list_.at(position).empty() ) {
+    single_mig_list_.at(position) = std::vector<MigEvent>(0);
   }
 
-  single_mig_probs_list_.at(position).at(getMigMatrixIndex(source_pop, sink_pop)) = fraction;
+  MigEvent migEvent = {source_pop, sink_pop, fraction};
+  single_mig_list_.at(position).push_back(migEvent);
+
   this->has_migration_ = true;
 }
 
@@ -486,13 +487,11 @@ std::ostream& operator<<(std::ostream& os, Model& model) {
       os << std::endl;
     }
 
-    for (size_t i = 0; i < n_pops; ++i) {
-      for (size_t j = 0; j < n_pops; ++j) {
-        if (model.single_mig_pop(i, j) != 0) {
-          os << " " << model.single_mig_pop(i, j) * 100 << "% of pop "
-             << i + 1 << " move to pop " << j + 1 << std::endl;
-        }
-      }
+    for (MigEvent me : model.single_mig_events()) {
+      os << " " 
+         << me.prob * 100 << "% of pop "
+         << me.source_pop + 1 << " move to pop " 
+         << me.sink_pop + 1 << std::endl;
     }
 
     if (idx < model.countChangeTimes() - 1) model.increaseTime();
@@ -632,7 +631,6 @@ void Model::addPopulation() {
 
   // Change Matrices
   addPopToMatrixList(mig_rates_list_, new_pop);
-  addPopToMatrixList(single_mig_probs_list_, new_pop, 0);
 }
 
 

src/model.h

@@ -50,6 +50,12 @@
 
 class Param;
 
+struct MigEvent {
+  size_t source_pop;
+  size_t sink_pop;
+  double prob;
+};
+
 enum SeqScale { relative, absolute, ms };
 
 class Model
@@ -67,7 +73,6 @@ class Model
    Model();
    Model(size_t sample_size);
    
-   
    // Default values;
    constexpr static double default_pop_size_ = 10000.0;
    constexpr static double default_growth_rate = 0.0;
@@ -223,10 +228,8 @@ class Model
     *
     * @return The probability/fraction of migration.
     */
-   double single_mig_pop(const size_t source, const size_t sink) const {
-    if (single_mig_probs_list_.at(current_time_idx_).empty()) return 0.0;
-    if (sink == source) return 0.0;
-    return single_mig_probs_list_.at(current_time_idx_).at( getMigMatrixIndex(source, sink) ); 
+   std::vector<MigEvent> single_mig_events() const {
+    return single_mig_list_.at(current_time_idx_);
    }
 
    void setMutationRate(double rate,
@@ -240,7 +243,7 @@ class Model
                              const double seq_position = 0);
 
    bool hasFixedTimeEvent(const double at_time) const {
-     if (single_mig_probs_list_.at(current_time_idx_).empty()) return false; 
+     if (single_mig_list_.at(current_time_idx_).empty()) return false; 
      if (getCurrentTime() != at_time) return false;
      return true;
    }
@@ -473,7 +476,8 @@ class Model
    std::vector<std::vector<double> > growth_rates_list_;
    std::vector<std::vector<double> > mig_rates_list_;
    std::vector<std::vector<double> > total_mig_rates_list_;
-   std::vector<std::vector<double> > single_mig_probs_list_;
+
+   std::vector<std::vector<MigEvent> > single_mig_list_;
 
    // Population sizes are saved as 1/(2N), where N is the actual population
    // size (do to fast multiplication rather than slow division in the

src/node_container.cc

@@ -241,6 +241,9 @@ void NodeContainer::clear() {
 
   // Clear free_slots_
   std::stack<Node*>().swap(free_slots_);
+  
+  // Delete nodes
+  for (std::vector<Node>* lane : node_lanes_) lane->clear();
 }
 
 

src/node_container.h

@@ -83,8 +83,9 @@ class NodeContainer {
         node_lanes_.push_back(new_lane);
       }
     }
-    (*node_lanes_.at(lane_counter_))[node_counter_] = Node(height, label);
-    return &*(node_lanes_[lane_counter_]->begin() + node_counter_++);
+    ++node_counter_;
+    node_lanes_.at(lane_counter_)->push_back(Node(height, label));
+    return &*(node_lanes_.at(lane_counter_)->end() - 1);
   }
 
   // Create Nodes
@@ -107,8 +108,9 @@ class NodeContainer {
         node_lanes_.push_back(new_lane);
       }
     }
-    (*node_lanes_.at(lane_counter_))[node_counter_] = Node(copiedNode);
-    return &*(node_lanes_[lane_counter_]->begin() + node_counter_++);
+    ++node_counter_;
+    node_lanes_.at(lane_counter_)->push_back(copiedNode);
+    return &*(node_lanes_.at(lane_counter_)->end() - 1);
   }
 
   void push_back(Node* node);

src/param.cc

@@ -53,6 +53,9 @@ Model Param::parse() {
        sfs = false,
        transpose = false;
 
+  // Tracks if demographic where added to the model.
+  // After the first demographic feature, defining substructure is no longer allowed.
+  bool has_demographic_feature = false;
 
   // The minimal time at which -eM, -eN, -eG, -eI, -ema and -es are allowed to happen. Is
   // increased by using -es.
@@ -115,6 +118,11 @@ Model Param::parse() {
     // ------------------------------------------------------------------
     // Set number of subpopulations and samples at time 0
     else if (*argv_i == "-I") {
+      // Check that -I is used immediately of the first two mandatory arguments
+      if (has_demographic_feature) {
+        throw std::invalid_argument("Option '-I' must be used before demographic '-M', '-N' or '-G' is used.");
+      }
+
       model.set_population_number(readNextInt());
       std::vector<size_t> sample_size;
       for (size_t i = 0; i < model.population_number(); ++i) {
@@ -147,6 +155,7 @@ Model Param::parse() {
     // Populations sizes 
     // ------------------------------------------------------------------
     else if (*argv_i == "-eN" || *argv_i == "-N") {
+      has_demographic_feature = true;
       if (*argv_i == "-eN") time = readNextInput<double>();
       else time = 0.0;
       if (time < min_time) {
@@ -169,6 +178,7 @@ Model Param::parse() {
     // Exponential Growth 
     // ------------------------------------------------------------------
     else if (*argv_i == "-G" || *argv_i == "-eG") {
+      has_demographic_feature = true;
       if (*argv_i == "-eG") time = readNextInput<double>();
       else time = 0.0;
       if (time < min_time) {
@@ -226,6 +236,7 @@ Model Param::parse() {
     }
 
     else if (*argv_i == "-M" || *argv_i == "-eM") {
+      has_demographic_feature = true;
       if (*argv_i == "-eM") {
         time = readNextInput<double>();
       }
@@ -234,6 +245,9 @@ Model Param::parse() {
         throw std::invalid_argument(std::string("If you use '-M' or '-eM' in a model with population merges ('-es'),") +
                                     std::string("then you need to sort both arguments by time."));
       }
+      if (model.population_number() == 1) {
+        throw std::invalid_argument("-M and -eM options can not be used if there is just one population");
+      }
       model.addSymmetricMigration(time, readNextInput<double>()/(model.population_number()-1), true, true);
     }
 

src/summary_statistics/oriented_forest.cc

@@ -66,10 +66,21 @@ void OrientedForest::generateTreeData(Node const* node, size_t &pos, int parent_
   parents_.at(pos) = parent_pos;
   parent_pos = pos--;
 
-  if (node->getLocalChild1() != NULL) {
-    generateTreeData(node->getLocalChild1(), pos, parent_pos+1, scaling_factor);
-    if (node->getLocalChild2() != NULL) {
-      generateTreeData(node->getLocalChild2(), pos, parent_pos+1, scaling_factor);
+  Node* local_child_1 = node->getLocalChild1();
+
+  if (local_child_1 != NULL) {
+    Node* local_child_2 = node->getLocalChild2();
+
+    if (local_child_2 != NULL) {
+      // Ensure that identical topologies lead to identical labels of nodes
+      if (local_child_2->height() > local_child_1->height()) {
+        Node* tmp = local_child_1;
+        local_child_1 = local_child_2;
+        local_child_2 = tmp;
       }
+      generateTreeData(local_child_2, pos, parent_pos+1, scaling_factor);
+    }
+
+    generateTreeData(local_child_1, pos, parent_pos+1, scaling_factor);
   }
 }

tests/algorithmtest/test_algorithm.cc

@@ -41,9 +41,9 @@ class TestAlgorithm : public CppUnit::TestCase {
     
     std::cout << "." << std::flush;
 
-    Forest forest = Forest(&model, this->rg);
     for (size_t i = 0; i < replicates; ++i) {
       // Check initial tree
+      Forest forest = Forest(&model, this->rg);
       forest.buildInitialTree();
       tmrca[0] += forest.getTMRCA(true);
       tree_length[0] += forest.getLocalTreeLength(true);
@@ -58,9 +58,6 @@ class TestAlgorithm : public CppUnit::TestCase {
           tree_length[j] += forest.getLocalTreeLength(true);
         }
       }
-
-      // Clear Forest
-      forest.clear();
     }
 
     // Allow an relative error of 2.5%. It would be nice to calculate

tests/test-coverage-local.sh

+#!/bin/bash
+
+# Generate html report
+lcov --base-directory . --directory . --zerocounters -q
+make check -mj
+lcov --base-directory . --directory . -c -o libbash_test.info
+# --rc lcov_branch_coverage=1 option will turn on branch check
+lcov --remove libbash_test.info "/usr*" -o libbash_test.info # remove output for external libraries
+rm -rf ./testCoverage
+genhtml -o ./testCoverage -t "libbash test coverage" --num-spaces 4 libbash_test.info --legend --function-coverage --demangle-cpp

tests/test_binaries.sh

 #!/bin/bash
-#
 # Test the binaries using scrms build in debug checks 
-# 
-# Author:   Paul R. Staab 
-# Email:    staab (at) bio.lmu.de
-# Licence:  GPLv3 or later
-#
+
+if [ "$1" == "osx" ]; then
+  echo "On macOS => Skipping valgrind checks"
+  valgrind=0
+else
+  valgrind=1
+fi
+
+make scrm scrm_dbg || exit 1
+make scrm_asan 2> /dev/null
+supports_sanitizers=$?
+if [ "$supports_sanitizers" == "0" ]; then
+  echo "Using address sanitizers"
+fi
 
 function test_scrm {
   echo -n " scrm $@ "
@@ -21,13 +29,24 @@ function test_scrm {
       exit 1
     fi
 
+    if [ "$supports_sanitizers" == "0" ]; then
+      ./scrm_asan $@ -seed $i > /dev/null
+      if [ $? -ne 0 ]; then
+        echo ""
+        echo "ASAN error in \"./scrm $@ -seed $i\""
+        exit 1
+      fi
+    fi
+
     # Test for memory leaks
+    if [ $valgrind -ne 0 ]; then
       valgrind --error-exitcode=1 --leak-check=full -q ./scrm $@ -seed $i > /dev/null
       if [ $? -ne 0 ]; then
         echo ""
         echo "Valgrind check of \"./scrm $@ -seed $i\" failed."
         exit 1
       fi
+    fi
 
   done
   echo " done."
@@ -61,13 +80,13 @@ echo ""
 echo "Testing Migration"
  test_scrm 5 5 -r 5 100 -I 2 3 2 1.2 || exit 1
  test_scrm 10 2 -r 20 200 -I 5 2 2 2 2 2 0.75 -l 5 || exit 1
- test_scrm 10 2 -r 10 100 -I 2 7 3 0.5 -eM 0.3 1.1 --print-model || exit 1
+ test_scrm 10 2 -r 10 100 -I 2 7 3 0.5 -eM 0.3 1.1 -O || exit 1
  test_scrm 10 2 -r 10 100 -I 2 7 3 -m 1 2 0.3 -em 0.5 2 1 0.6 -eM 2.0 1 || exit 1
  test_scrm 20 2 -I 3 2 2 2 1.0 -eI 1.0 2 2 2 -eI 2.0 2 3 3 || exit 1
 echo ""
 
 echo "Testing Size Change"
-  test_scrm 10 2 -r 1 100 -I 3 3 3 4 0.5 -eN 0.1 0.05 -eN 0.2 0.5 --print-model || exit 1 
+  test_scrm 10 2 -r 1 100 -I 3 3 3 4 0.5 -eN 0.1 0.05 -eN 0.2 0.5 -O || exit 1 
   test_scrm 10 2 -r 10 100 -I 3 3 3 4 0.5 -eN 0.1 0.05 -eN 0.2 0.5 -l 10 || exit 1 
 echo ""