[med-svn] [scrm] 01/04: Imported Upstream version 1.7.0
Kevin Murray
daube-guest at moszumanska.debian.org
Wed Feb 10 01:51:06 UTC 2016
This is an automated email from the git hooks/post-receive script.
daube-guest pushed a commit to branch master
in repository scrm.
commit 41e78d9a19481ba47dfb09676c1f2a5eeb23fb58
Author: Kevin Murray <spam at kdmurray.id.au>
Date: Tue Feb 9 17:46:16 2016 -0800
Imported Upstream version 1.7.0
---
.gitignore | 1 +
.travis.yml | 36 ++--
.travis/build_win_binaries.sh | 2 +
.travis/setup.sh | 24 +++
Makefile.am | 32 +--
NEWS.md | 26 +++
README.md | 9 +-
configure.ac | 12 +-
src/contemporaries_container.h | 116 +++++------
src/forest-debug.cc | 1 -
src/forest.cc | 241 +++++++++++++----------
src/forest.h | 52 ++---
src/model.cc | 212 ++++++++++----------
src/model.h | 9 +-
src/node_container.cc | 46 +++--
src/node_container.h | 130 +++++++-----
src/param.cc | 10 +-
src/random/mersenne_twister.cc | 7 -
src/random/mersenne_twister.h | 15 +-
src/random/random_generator.cc | 2 +-
src/random/random_generator.h | 27 +--
src/scrm.cc | 12 +-
src/summary_statistics/newick_tree.cc | 6 +-
src/summary_statistics/newick_tree.h | 1 -
src/summary_statistics/oriented_forest.h | 1 -
src/summary_statistics/seg_sites.cc | 36 +++-
src/summary_statistics/seg_sites.h | 7 +-
src/summary_statistics/summary_statistic.h | 4 +-
src/summary_statistics/tmrca.cc | 11 +-
src/time_interval.h | 9 +-
tests/algorithmtest/test_algorithm.cc | 35 +---
tests/test_binaries.sh | 46 ++---
tests/unittests/test_contemporaries_container.cc | 21 +-
tests/unittests/test_forest.cc | 171 ++++++++--------
tests/unittests/test_model.cc | 168 ++++++++--------
tests/unittests/test_node.cc | 8 +-
tests/unittests/test_node_container.cc | 24 +--
tests/unittests/test_param.cc | 89 +++++----
tests/unittests/test_random_generator.cc | 16 +-
tests/unittests/test_summary_statistics.cc | 45 ++++-
tests/unittests/test_time_interval.cc | 63 +++---
tools/compare_versions.sh | 6 +-
42 files changed, 980 insertions(+), 809 deletions(-)
diff --git a/.gitignore b/.gitignore
index a6a5ea2..37cbfd9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -86,3 +86,4 @@ algorithm_tests
unit_tests
stuff
.dirstamp
+configure.lineno
diff --git a/.travis.yml b/.travis.yml
index bec3864..5da88cd 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -1,27 +1,30 @@
language: cpp
os:
-- linux
-- osx
+ - linux
+ - osx
compiler:
-- gcc
-- clang
+ - gcc
+ - clang
+env:
+ - CXXFLAGS="-Werror"
+ - CXXFLAGS="-m32 -Werror"
+matrix:
+ exclude:
+ - os: osx
+ env: CXXFLAGS="-m32 -Werror"
sudo: required
-before_install:
-- if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo add-apt-repository -y ppa:dns/gnu; fi
-- if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo add-apt-repository -y ppa:ubuntu-toolchain-r/test; fi
-- if [ "$TRAVIS_OS_NAME" == "linux" ]; then sudo apt-get update -qq && sudo apt-get install -qqy automake libcppunit-dev valgrind; fi
-- if [ "$TRAVIS_OS_NAME" == "linux" -a "$CC" = "gcc" ]; then sudo apt-get install -qqy g++-5 && sudo update-alternatives --install /usr/bin/g++ g++ /usr/bin/g++-5 50; fi
-- if [ "$TRAVIS_OS_NAME" == "osx" ]; then brew update && brew install cppunit valgrind; fi
+dist: trusty
+before_install: .travis/setup.sh "$TRAVIS_OS_NAME" "$CXX" "$CXXFLAGS"
before_script:
-- ./bootstrap
+ - ./bootstrap
script:
-- make scrm scrm_dbg
-- make unittests
-- ./tests/test_binaries.sh
-- make algorithmtest
+ - make unittests
+ - make scrm scrm_dbg && ./tests/test_binaries.sh
+ - make algorithmtest
before_deploy:
- ./.travis/build_src_pkg.sh
- ./.travis/build_static_binaries.sh
+ - ./.travis/build_win_binaries.sh
deploy:
provider: releases
api_key:
@@ -29,7 +32,10 @@ deploy:
file:
- "scrm-src.tar.gz"
- "scrm-x64-static.tar.gz"
+ - "scrm-win32.zip"
+ - "scrm-win64.zip"
on:
tags: true
condition: "$CC = gcc"
+ condition: "$CXXFLAGS == -Werror"
condition: "$TRAVIS_OS_NAME == linux"
diff --git a/.travis/build_win_binaries.sh b/.travis/build_win_binaries.sh
index 0892066..c36587b 100755
--- a/.travis/build_win_binaries.sh
+++ b/.travis/build_win_binaries.sh
@@ -1,5 +1,7 @@
#!/bin/bash
+sudo apt-get install -qqy mingw-w64 || exit 1
+
# 64 bit
CXX=i686-w64-mingw32-g++ CXXFLAGS='-O3' LDFLAGS='-static-libgcc -static -lpthread' \
./configure --host=i686-w64-mingw32 || exit 1
diff --git a/.travis/setup.sh b/.travis/setup.sh
new file mode 100755
index 0000000..2eaefa9
--- /dev/null
+++ b/.travis/setup.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+os=$1
+cxx=$2
+cxxflags=$3
+
+echo "os: $os; cxx: $cxx; cxxflags: $cxxflags"
+
+if [ "$os" == "linux" ]; then
+ if [[ "$cxxflags" == *"-m32"* ]]; then
+ sudo dpkg --add-architecture i386
+ sudo apt-get update -qq
+ sudo apt-get install -qqy automake valgrind g++-multilib libc6-dbg:i386 libcppunit-dev:i386
+ else
+ sudo apt-get update -qq
+ sudo apt-get install -qqy automake libcppunit-dev valgrind;
+ fi
+fi
+
+if [ "$os" == "osx" ]; then
+ brew update
+ brew install cppunit valgrind
+fi
+
diff --git a/Makefile.am b/Makefile.am
index eb95f96..806abfd 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,22 +1,22 @@
-EXTRA_DIST = doc/create-manual.sh doc/knitr.css doc/scrm.1
-bin_PROGRAMS = scrm
+EXTRA_DIST = doc/create-manual.sh doc/knitr.css doc/scrm.1
+bin_PROGRAMS = scrm
man_MANS = doc/scrm.1
-TESTS = unit_tests algorithm_tests
+TESTS = unit_tests algorithm_tests
check_PROGRAMS = unit_tests algorithm_tests scrm_dbg scrm_prof
PROG = SCRM
-dist-hook:
+dist-hook:
chmod u+w $(distdir)/doc
cd $(distdir); ./doc/create-manual.sh $(VERSION)
test: unittests
unittests: unit_tests
- ./unit_tests
+ ./unit_tests
algorithmtest: algorithm_tests
- ./algorithm_tests
+ ./algorithm_tests
nrml_src = src/param.cc src/forest.cc src/node.cc src/node_container.cc src/time_interval.cc \
src/model.cc src/tree_point.cc \
@@ -27,7 +27,7 @@ nrml_src = src/param.cc src/forest.cc src/node.cc src/node_container.cc src/time
random_src = src/random/random_generator.cc src/random/mersenne_twister.cc \
src/random/fastfunc.cc \
src/random/random_generator.h src/random/mersenne_twister.h \
- src/random/fastfunc.h
+ src/random/fastfunc.h
sumstat_src = src/summary_statistics/tmrca.cc \
src/summary_statistics/seg_sites.cc \
@@ -40,7 +40,7 @@ sumstat_src = src/summary_statistics/tmrca.cc \
src/summary_statistics/summary_statistic.h \
src/summary_statistics/oriented_forest.cc \
src/summary_statistics/oriented_forest.h
-
+
scrm_src = $(nrml_src) $(random_src) $(sumstat_src)
debug_src = src/random/constant_generator.cc src/random/constant_generator.h \
@@ -48,12 +48,12 @@ debug_src = src/random/constant_generator.cc src/random/constant_generator.h \
unit_test_src = tests/unittests/test_forest.cc tests/unittests/test_model.cc\
tests/unittests/test_node.cc tests/unittests/test_node_container.cc\
- tests/cppunit/test_runner.cc tests/unittests/test_time_interval.cc\
- tests/unittests/test_fastfunc.cc tests/unittests/test_param.cc\
+ tests/cppunit/test_runner.cc tests/unittests/test_time_interval.cc\
+ tests/unittests/test_fastfunc.cc tests/unittests/test_param.cc\
tests/unittests/test_random_generator.cc tests/unittests/test_summary_statistics.cc\
- tests/unittests/test_contemporaries_container.cc
+ tests/unittests/test_contemporaries_container.cc
-alg_test_src = tests/cppunit/test_runner.cc tests/algorithmtest/test_algorithm.cc
+alg_test_src = tests/cppunit/test_runner.cc tests/algorithmtest/test_algorithm.cc
scrm_SOURCES = $(scrm_src) src/scrm.cc
@@ -61,11 +61,11 @@ scrm_dbg_SOURCES = $(scrm_src) $(debug_src) src/scrm.cc
scrm_prof_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=
+scrm_dbg_CXXFLAGS= -g
scrm_prof_CXXFLAGS= -pg -DNDEBUG
-unit_tests_CXXFLAGS = -DUNITTEST -DNDEBUG
-algorithm_tests_CXXFLAGS = -DNDEBUG
+unit_tests_CXXFLAGS = -g -DUNITTEST -DNDEBUG
+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
diff --git a/NEWS.md b/NEWS.md
index 5082373..f655922 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,32 @@
scrm Version History
========================
+scrm 1.7.0
+------------------------
+Released: 2016-02-07
+
+### New Features
++ New command line option "--transpose-segsites" that prints the segregating
+ sites matrix with rows repesenting mutations instead of individuals. This
+ also adds the time at which the mutation occurred as additional information
+ (#85).
+
+### Improvements
++ scrm now throws an error when mutation or recombination rate changes
+ have invalid sequence positions (#82).
++ The length of segments between recombinations is now always reported
+ in non-scientific notation (#81).
++ Improved the error message when lines could not coalescence because
+ of missing migration or negative population growth (#87).
+
+### Bug Fixes
++ `Forest.coalescence_finished` has not set to `true` after a pairwise
+ coalescence event (#89). This had no effect on the command line version
+ of scrm.
++ Fixes scrm's tests suite on 32 bit systems (#98).
+
+
+
scrm 1.6.1
------------------------
Released: 2015-07-09
diff --git a/README.md b/README.md
index 8c3a885..4d2ad57 100644
--- a/README.md
+++ b/README.md
@@ -1,10 +1,11 @@
scrm
====
-_scrm_ is a coalescent simulator for biological sequences. Different to similar programs,
-it can approximate the coalescent with recombination as closely as needed, but still has
-only linear runtime cost for long sequences. It allows you to rapidly simulate chromosome
-scale sequences with an essentially correct genetic linkage structure.
+_scrm_ simulates the evolution of genetic sequences. It takes a neutral evolutionary model as input,
+and generates random sequences that evolved under the model. As coalescent simulator, it traces
+the ancestry of the sampled sequences backwards in time and is therefore extremely efficient. Compared to
+other coalescent simulators, it can simulate chromosome-scale sequences without a measureable reduction of
+genetic linkage between different sites.
## Installation
diff --git a/configure.ac b/configure.ac
index 98fdd15..73009df 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,4 +1,4 @@
-AC_INIT([scrm], [1.6.1],[https://github.com/paulstaab/scrm/issues])
+AC_INIT([scrm], [1.7.0],[https://github.com/paulstaab/scrm/issues])
AM_INIT_AUTOMAKE([subdir-objects -Wall -Werror foreign])
# Use -O3 as default optimization level
@@ -35,12 +35,20 @@ AC_TYPE_SIZE_T
# Enable Link-time optimization if supported (gcc only)
if test x$CXX = xg++; then
- AX_CHECK_COMPILE_FLAG([-flto], [OPT_CXXFLAGS="-flto"], [], [-Werror])
+ AX_CHECK_COMPILE_FLAG([-flto], [OPT_CXXFLAGS="$OPT_CXXFLAGS -flto"], [], [-Werror])
AC_SUBST(OPT_CXXFLAGS)
fi
# Enable Warnings
AX_CXXFLAGS_WARN_ALL
+AX_CHECK_COMPILE_FLAG([-Wextra], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Wextra"], [], [-Werror])
+# AX_CHECK_COMPILE_FLAG([-Wshadow], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Wshadow"], [], [-Werror])
+AX_CHECK_COMPILE_FLAG([-Wnon-virtual-dtor], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Wnon-virtual-dtor"], [], [-Werror])
+#AX_CHECK_COMPILE_FLAG([-Wold-style-cast], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Wold-style-cast"], [], [-Werror])
+AX_CHECK_COMPILE_FLAG([-Wcast-align], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Wcast-align"], [], [-Werror])
+AX_CHECK_COMPILE_FLAG([-Wunused], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Wunused"], [], [-Werror])
+AX_CHECK_COMPILE_FLAG([-Woverloaded-virtual], [OPT_CXXFLAGS="$OPT_CXXFLAGS -Woverloaded-virtual"], [], [-Werror])
+AX_CHECK_COMPILE_FLAG([-pedantic], [OPT_CXXFLAGS="$OPT_CXXFLAGS -pedantic"], [], [-Werror])
AC_CONFIG_FILES([Makefile])
AC_OUTPUT
diff --git a/src/contemporaries_container.h b/src/contemporaries_container.h
index 1281de1..75ff36c 100644
--- a/src/contemporaries_container.h
+++ b/src/contemporaries_container.h
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -35,39 +35,22 @@
class ContemporariesIterator {
public:
- ContemporariesIterator(std::unordered_set<Node*>::iterator it) {
- it_set_ = it;
- use_set_ = true;
- };
-
- ContemporariesIterator(std::vector<Node*>::iterator it) {
- it_vec_ = it;
- use_set_ = false;
- };
+ ContemporariesIterator(std::unordered_set<Node*>::iterator it) : it_set_(it), use_set_(true) { };
+ ContemporariesIterator(std::vector<Node*>::iterator it) : it_vec_(it), use_set_(false) { };
Node* operator*() const {
if (use_set_) return *it_set_;
else return *it_vec_;
}
- Node* operator++() {
- if (use_set_) {
- ++it_set_;
- return *it_set_;
- } else {
- ++it_vec_;
- return *it_vec_;
- }
+ ContemporariesIterator& operator++() {
+ if (use_set_) (++it_set_);
+ else ++it_vec_;
+ return *this;
}
- Node* operator++(int) {
- Node* node = **this;
- ++*this;
- return node;
- }
-
bool operator==(const ContemporariesIterator &other) const {
- if (use_set_) return (it_set_ == other.it_set_);
+ if (use_set_) return (it_set_ == other.it_set_);
else return (it_vec_ == other.it_vec_);
}
@@ -76,6 +59,8 @@ class ContemporariesIterator {
}
private:
+ ContemporariesIterator();
+
std::unordered_set<Node*>::iterator it_set_;
std::vector<Node*>::iterator it_vec_;
bool use_set_;
@@ -84,39 +69,23 @@ class ContemporariesIterator {
class ContemporariesConstIterator {
public:
- ContemporariesConstIterator(std::unordered_set<Node*>::const_iterator it) {
- it_set_ = it;
- use_set_ = true;
- };
-
- ContemporariesConstIterator(std::vector<Node*>::const_iterator it) {
- it_vec_ = it;
- use_set_ = false;
- };
+ ContemporariesConstIterator(std::unordered_set<Node*>::const_iterator it) : it_set_(it), use_set_(true) { };
+ ContemporariesConstIterator(std::vector<Node*>::const_iterator it) : it_vec_(it), use_set_(false) { };
Node const* operator*() const {
if (use_set_) return *it_set_;
else return *it_vec_;
}
- Node const* operator++() {
- if (use_set_) {
- ++it_set_;
- return *it_set_;
- } else {
- ++it_vec_;
- return *it_vec_;
- }
+ ContemporariesConstIterator& operator++() {
+ if (use_set_) ++it_set_;
+ else ++it_vec_;
+ return *this;
}
- Node const* operator++(int) {
- Node const* node = **this;
- ++*this;
- return node;
- }
-
bool operator==(const ContemporariesConstIterator &other) const {
- if (use_set_) return (it_set_ == other.it_set_);
+ if (use_set_ != other.use_set_) return false;
+ if (use_set_) return (it_set_ == other.it_set_);
else return (it_vec_ == other.it_vec_);
}
@@ -133,9 +102,10 @@ class ContemporariesConstIterator {
class ContemporariesContainer {
public:
ContemporariesContainer();
- ContemporariesContainer(const size_t pop_number,
+ ContemporariesContainer(const size_t pop_number,
const size_t sample_number,
RandomGenerator *rg);
+ ~ContemporariesContainer(){};
void add(Node* node);
void remove(Node *node);
@@ -151,30 +121,30 @@ class ContemporariesContainer {
bool use_set() const { return use_set_; };
// Create Iterators
- ContemporariesConstIterator begin(const size_t pop) const {
- if (use_set_) return ContemporariesConstIterator(contemporaries_set().at(pop).cbegin());
+ ContemporariesConstIterator begin(const size_t pop) const {
+ if (use_set_) return ContemporariesConstIterator(contemporaries_set().at(pop).cbegin());
else return ContemporariesConstIterator(contemporaries_vector().at(pop).cbegin());
}
- ContemporariesConstIterator end(const size_t pop) const {
- if (use_set_) return ContemporariesConstIterator(contemporaries_set().at(pop).cend());
+ ContemporariesConstIterator end(const size_t pop) const {
+ if (use_set_) return ContemporariesConstIterator(contemporaries_set().at(pop).cend());
else return ContemporariesConstIterator(contemporaries_vector().at(pop).cend());
}
// Create Iterators for the buffer
- ContemporariesConstIterator buffer_begin(const size_t pop) const {
- if (use_set_) return ContemporariesConstIterator(buffer_set().at(pop).cbegin());
+ ContemporariesConstIterator buffer_begin(const size_t pop) const {
+ if (use_set_) return ContemporariesConstIterator(buffer_set().at(pop).cbegin());
else return ContemporariesConstIterator(buffer_vector().at(pop).cbegin());
}
- ContemporariesConstIterator buffer_end(const size_t pop) const {
- if (use_set_) return ContemporariesConstIterator(buffer_set().at(pop).cend());
+ ContemporariesConstIterator buffer_end(const size_t pop) const {
+ if (use_set_) return ContemporariesConstIterator(buffer_set().at(pop).cend());
else return ContemporariesConstIterator(buffer_vector().at(pop).cend());
}
- ContemporariesIterator buffer_begin(const size_t pop) {
- if (use_set_) return ContemporariesIterator(buffer_set().at(pop).begin());
+ ContemporariesIterator buffer_begin(const size_t pop) {
+ if (use_set_) return ContemporariesIterator(buffer_set().at(pop).begin());
else return ContemporariesIterator(buffer_vector().at(pop).begin());
}
- ContemporariesIterator buffer_end(const size_t pop) {
- if (use_set_) return ContemporariesIterator(buffer_set().at(pop).end());
+ ContemporariesIterator buffer_end(const size_t pop) {
+ if (use_set_) return ContemporariesIterator(buffer_set().at(pop).end());
else return ContemporariesIterator(buffer_vector().at(pop).end());
}
@@ -223,6 +193,7 @@ class ContemporariesContainer {
RandomGenerator* rg_;
};
+
inline ContemporariesContainer::ContemporariesContainer() {
contemporaries_vec1_ = std::vector<std::vector<Node*> >(1, std::vector<Node*>(100));
contemporaries_vec2_ = std::vector<std::vector<Node*> >(1, std::vector<Node*>(100));
@@ -233,13 +204,14 @@ inline ContemporariesContainer::ContemporariesContainer() {
use_set_ = false;
}
-inline ContemporariesContainer::ContemporariesContainer(const size_t pop_number,
+
+inline ContemporariesContainer::ContemporariesContainer(const size_t pop_number,
const size_t sample_number,
RandomGenerator* rg) {
// Use vectors for the storage if the number of samples is below 750.
// This threshold is mostly arbitrary, with simulation supporting it in
- // special situations.
+ // special situations.
if (sample_number <= 750) {
contemporaries_vec1_ = std::vector<std::vector<Node*> >(pop_number);
for ( auto it : contemporaries_vec1_ ) it.reserve(sample_number + 200);
@@ -279,7 +251,7 @@ inline void ContemporariesContainer::remove(Node* node) {
}
inline void ContemporariesContainer::replaceChildren(Node *add_node) {
- replace(add_node, add_node->first_child(), add_node->second_child());
+ replace(add_node, add_node->first_child(), add_node->second_child());
}
inline void ContemporariesContainer::replace(Node *add_node, Node *del_node_1, Node *del_node_2) {
@@ -308,18 +280,18 @@ inline size_t ContemporariesContainer::size(const size_t pop) const {
/**
* @brief Function that buffers the current state of contemporaries for later
- * use.
+ * use.
*
* Be careful not to change the tree at the time of buffer, as this
* will not be reflected in the buffered contemporaries.
*
* @param current_time The time for with the current state is valid.
*
- * @return
+ * @return
*/
inline void ContemporariesContainer::buffer(const double current_time) {
buffer_time_ = current_time;
- use_first_ = 1 - use_first_;
+ use_first_ = 1 - use_first_;
this->clear();
}
@@ -335,7 +307,7 @@ inline Node* ContemporariesContainer::sample(const size_t pop) const {
if (use_set_) {
// Sample the position of the Node we return
- for (auto it = contemporaries_set().at(pop).begin();
+ for (auto it = contemporaries_set().at(pop).begin();
it != contemporaries_set().at(pop).end(); ++it) {
assert( *it != NULL );
if ( sample == 0 ) return (*it);
diff --git a/src/forest-debug.cc b/src/forest-debug.cc
index 421f95f..ab83be2 100644
--- a/src/forest-debug.cc
+++ b/src/forest-debug.cc
@@ -32,7 +32,6 @@ void Forest::createExampleTree() {
this->rec_bases_.push_back(5.0);
this->current_rec_ = 1;
- //this->rec_bases_.push_back(105.0);
Node* leaf1 = nodes()->createNode(0, 1);
Node* leaf2 = nodes()->createNode(0, 2);
diff --git a/src/forest.cc b/src/forest.cc
index 436353f..efaaef1 100644
--- a/src/forest.cc
+++ b/src/forest.cc
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -32,7 +32,7 @@ Forest::Forest(Model* model, RandomGenerator* random_generator) {
}
// Sets member variable to default values
-void Forest::initialize(Model* model,
+void Forest::initialize(Model* model,
RandomGenerator* rg) {
model->resetTime();
@@ -49,9 +49,10 @@ void Forest::initialize(Model* model,
this->coalescence_finished_ = true;
- this->contemporaries_ = ContemporariesContainer(model->population_number(),
+ this->contemporaries_ = ContemporariesContainer(model->population_number(),
model->sample_size(),
rg);
+
tmp_event_time_ = -1;
}
@@ -59,7 +60,7 @@ void Forest::initialize(Model* model,
* @brief Copy constructor for forest
*
* This creates a copy of an Forest. It only produces valid results
- * when there is no ongoing coalescence event (e.g. when we are not
+ * when there is no ongoing coalescence event (e.g. when we are not
* inside of sampleNextGenealogy()).
*
* Also, both copies of it share the model and the random generator,
@@ -68,7 +69,7 @@ void Forest::initialize(Model* model,
*
* @param current_forest Forest that needs to be duplicated
*/
-Forest::Forest(const Forest ¤t_forest) {
+Forest::Forest(const Forest ¤t_forest) {
if (!current_forest.coalescence_finished_) {
throw std::logic_error("Can not copy forest during an ongoing coalescence");
}
@@ -93,10 +94,12 @@ Forest::Forest(const Forest ¤t_forest) {
}
// Set initial values for temporary variables
- this->contemporaries_ = ContemporariesContainer(model().population_number(),
+ this->contemporaries_ = ContemporariesContainer(model().population_number(),
model().sample_size(),
random_generator());
- this->tmp_event_time_ = -1;
+
+ this->tmp_event_time_ = this->getTMRCA(false); // Disable buffer for next genealogy.
+
this->coalescence_finished_ = true;
dout<<" #################### check copied forest ###############"<<std::endl;
@@ -108,12 +111,12 @@ Forest::Forest(const Forest ¤t_forest) {
}
-/**
+/**
* function that cuts a subtree out of a tree of the forest and reinserts it as
* a separate tree.
*
* This is primarily used to cut the subtree below an recombination
- * away.
+ * away.
*
* \param cut_point A TreePoint marking the top of the subtree to cut.
* \return The root of the now separated subtree.
@@ -125,7 +128,7 @@ Node* Forest::cut(const TreePoint &cut_point) {
//The new end of the old branch after the cut
Node* new_leaf = nodes()->createNode(cut_point.height());
-
+
if ( !cut_point.base_node()->local() )
new_leaf->make_nonlocal(cut_point.base_node()->last_update());
else
@@ -156,7 +159,7 @@ Node* Forest::cut(const TreePoint &cut_point) {
new_root->set_first_child(cut_point.base_node());
// Set invariants of new root
- new_root->set_length_below(cut_point.base_node()->length_below() +
+ new_root->set_length_below(cut_point.base_node()->length_below() +
cut_point.relative_height() );
new_root->set_samples_below(cut_point.base_node()->samples_below() );
@@ -177,7 +180,7 @@ Node* Forest::cut(const TreePoint &cut_point) {
/**
- * Function to update the invariants (local, samples_below, length_below)
+ * Function to update the invariants (local, samples_below, length_below)
* of a 'node' and all of its (grand-)parents. Also sets local_root_ if it
* encounters it. Never makes non-local nodes local, only the other way round.
*
@@ -192,11 +195,11 @@ Node* Forest::cut(const TreePoint &cut_point) {
* \param invariants_only If true, it only updates the nodes invariants, but
* does not make nodes non-local and change the local root.
*/
-void Forest::updateAbove(Node* node, bool above_local_root,
+void Forest::updateAbove(Node* node, bool above_local_root,
const bool &recursive, const bool &invariants_only) {
//dout << "Updating: " << node << " above_local_root: " << above_local_root << std::endl;
-
+
// Fast forward above local root because this part is fairly straight forward
if (above_local_root) {
// Assure that everything is non-local
@@ -208,7 +211,7 @@ void Forest::updateAbove(Node* node, bool above_local_root,
set_primary_root(node);
}
return;
- }
+ }
if ( recursive ) updateAbove(node->parent(), true, true);
return;
}
@@ -237,7 +240,7 @@ void Forest::updateAbove(Node* node, bool above_local_root,
}
assert( length_below >= 0 );
- // Update whether the node is local or not
+ // Update whether the node is local or not
if (!invariants_only) {
if (samples_below == 0) {
if ( node->local() ) node->make_nonlocal(current_rec());
@@ -246,7 +249,7 @@ void Forest::updateAbove(Node* node, bool above_local_root,
if ( node->local() ) node->make_nonlocal(current_rec());
// Are we the local root?
- if (node->countChildren() == 2 &&
+ if (node->countChildren() == 2 &&
l_child->samples_below() > 0 && h_child->samples_below() > 0) {
set_local_root(node);
}
@@ -258,7 +261,7 @@ void Forest::updateAbove(Node* node, bool above_local_root,
}
// If nothing changed, we also don't need to update the tree further above...
- if (samples_below == node->samples_below() &&
+ if (samples_below == node->samples_below() &&
areSame(length_below, node->length_below()) ) {
return;
}
@@ -318,6 +321,7 @@ void Forest::buildInitialTree() {
assert(this->checkLeafsOnLocalTree());
assert(this->printTree());
assert(this->printNodes());
+ assert(this->coalescence_finished());
}
this->sampleNextBase();
dout << "Next Sequence position: " << this->next_base() << std::endl;
@@ -385,12 +389,12 @@ TreePoint Forest::samplePoint(Node* node, double length_left) const {
double tmp = node->first_child()->height_above() + node->first_child()->length_below();
if ( length_left <= tmp )
return samplePoint(node->first_child(), length_left);
- else
+ else
return samplePoint(node->second_child(), length_left - tmp);
}
/* Alternative inefficient implementation
TreePoint Forest::samplePoint(Node* node, double length_left) {
- length_left = random_generator()->sample() * local_tree_length();
+ length_left = random_generator()->sample() * local_tree_length();
for (auto ni = nodes()->iterator(); ni.good(); ++ni) {
if (!(*ni)->local()) continue;
if (length_left < (*ni)->height_above()) return TreePoint(*ni, length_left, true);
@@ -402,7 +406,7 @@ TreePoint Forest::samplePoint(Node* node, double length_left) {
-/**
+/**
* Function to modify the tree after we encountered a recombination on the
* sequence. Also samples a place for this recombination on the tree, marks the
* branch above as non-local (and updates invariants) if needed, cuts the
@@ -452,6 +456,7 @@ void Forest::sampleNextGenealogy() {
assert( this->checkTree() );
assert( this->printTree() );
assert( this->printNodes() );
+ assert( this->coalescence_finished() );
this->sampleNextBase();
dout << "Next Sequence position: " << this->next_base() << std::endl;
@@ -459,7 +464,7 @@ void Forest::sampleNextGenealogy() {
}
-/**
+/**
* Function for doing a coalescence.
*
* \param start_node The node at which the coalescence starts. Must be the root
@@ -477,8 +482,8 @@ void Forest::sampleCoalescences(Node *start_node) {
coalescence_finished_ = false;
// This assertion needs an exception for building the initial tree
- assert ( current_rec() == 1 ||
- active_node(1)->in_sample() ||
+ assert ( current_rec() == 1 ||
+ active_node(1)->in_sample() ||
start_node->height() <= active_node(1)->height() );
// Only prune every second round
@@ -487,10 +492,10 @@ void Forest::sampleCoalescences(Node *start_node) {
dout << "* * Time interval: " << (*ti).start_height() << " - "
<< (*ti).end_height() << " (Last event at " << tmp_event_.time() << ")" << std::endl;
- // Assert that we don't accidentally jump in time
+ // Assert that we don't accidentally jump in time
assert( tmp_event_.time() < 0 || tmp_event_.time() == (*ti).start_height() );
- // Update States & Rates (see their declaration for explanation);
+ // Update States & Rates (see their declaration for explanation);
states_[0] = getNodeState(active_node(0), (*ti).start_height());
states_[1] = getNodeState(active_node(1), (*ti).start_height());
@@ -503,10 +508,10 @@ void Forest::sampleCoalescences(Node *start_node) {
// Some debug checks
dout << "* * * Active Nodes: a0:" << active_node(0) << ":s" << states_[0]
- << "(p" << active_node(0)->population() << ")"
- << " a1:" << active_node(1) << ":s" << states_[1]
+ << "(p" << active_node(0)->population() << ")"
+ << " a1:" << active_node(1) << ":s" << states_[1]
<< "(p" << active_node(1)->population() << ")" << std::endl
- << "* * * Total Rates: " << rates_[0] << " "
+ << "* * * Total Rates: " << rates_[0] << " "
<< rates_[1] << " " << rates_[2] << std::endl;
assert( active_node(0) != active_node(1) );
@@ -517,7 +522,7 @@ void Forest::sampleCoalescences(Node *start_node) {
assert( states_[1] == 1 || active_node(1)->parent_height() >= tmp_event_.time() );
assert( states_[0] != 2 || !active_node(0)->local() );
assert( states_[1] != 2 || !active_node(1)->local() );
-
+
assert( active_node(0)->first_child() == NULL || active_node(0)->first_child()->local() ||
active_node(0)->second_child() == NULL || active_node(0)->second_child()->local() );
assert( active_node(1)->first_child() == NULL || active_node(1)->first_child()->local() ||
@@ -530,16 +535,14 @@ void Forest::sampleCoalescences(Node *start_node) {
assert( tmp_event_.isNoEvent() || tmp_event_.time() <= (*ti).end_height() );
- // Go on if nothing happens in this time interval
+ // Implement the event
if ( tmp_event_.isNoEvent() ) {
this->implementNoEvent(*ti, coalescence_finished_);
- if (coalescence_finished_) return;
}
- // First take care of pairwise coalescence
else if ( tmp_event_.isPwCoalescence() ) {
this->implementPwCoalescence(active_node(0), active_node(1), tmp_event_.time());
- return;
+ this->coalescence_finished_ = true;
}
else if ( tmp_event_.isRecombination() ) {
@@ -548,17 +551,17 @@ void Forest::sampleCoalescences(Node *start_node) {
else if ( tmp_event_.isMigration() ) {
this->implementMigration(tmp_event_, true, ti);
- assert( this->printTree() );
}
else if ( tmp_event_.isCoalescence() ) {
this->implementCoalescence(tmp_event_, ti);
assert( checkInvariants(tmp_event_.node()) );
- if (coalescence_finished_) return;
assert( this->printTree() );
}
+
+ if (coalescence_finished()) return;
}
-}
+}
void Forest::calcRates(const TimeInterval &ti) {
@@ -572,8 +575,8 @@ void Forest::calcRates(const TimeInterval &ti) {
if (states_[0] == 1) {
// coalescing or migrating
rates_[0] += model().total_migration_rate(active_node(0)->population());
- if (model().growth_rate(active_node(0)->population()) == 0.0)
- rates_[0] += calcCoalescenceRate(active_node(0)->population(), ti);
+ if (model().growth_rate(active_node(0)->population()) == 0.0)
+ rates_[0] += calcCoalescenceRate(active_node(0)->population(), ti);
else {
// exponential growth -- assign this node to timeline 1
rates_[1] += calcCoalescenceRate(active_node(0)->population(), ti);
@@ -582,10 +585,10 @@ void Forest::calcRates(const TimeInterval &ti) {
}
else if (states_[0] == 2) {
// recombining
- rates_[0] += calcRecombinationRate(active_node(0));
+ rates_[0] += calcRecombinationRate(active_node(0));
}
- // The second node is a bit more complicated
+ // The second node is a bit more complicated
if (states_[1] == 1) {
// coalescing or migrating
rates_[0] += model().total_migration_rate(active_node(1)->population());
@@ -593,9 +596,9 @@ void Forest::calcRates(const TimeInterval &ti) {
// No Growth => Normal time
rates_[0] += calcCoalescenceRate(active_node(1)->population(), ti);
- if (states_[0] == 1 && active_node(0)->population() == active_node(1)->population()) {
+ if (states_[0] == 1 && active_node(0)->population() == active_node(1)->population()) {
// Also add rates for pw coalescence
- rates_[0] += calcPwCoalescenceRate(active_node(1)->population(), ti);
+ rates_[0] += calcPwCoalescenceRate(active_node(1)->population(), ti);
}
}
else {
@@ -604,9 +607,9 @@ void Forest::calcRates(const TimeInterval &ti) {
// Coalescing or migrating; and we can use the timeline of the first node
rates_[1] += calcCoalescenceRate(active_node(1)->population(), ti);
// And must add pw coalescence again
- rates_[1] += calcPwCoalescenceRate(active_node(1)->population(), ti);
+ rates_[1] += calcPwCoalescenceRate(active_node(1)->population(), ti);
active_nodes_timelines_[1] = 1;
- }
+ }
else {
// No chance of a pairwise coalescence, but there is growth.
// We might need our own timeline (This could be made more efficient if both populations have
@@ -624,15 +627,15 @@ void Forest::calcRates(const TimeInterval &ti) {
assert(rates_[0] >= 0);
assert(rates_[1] >= 0);
assert(rates_[2] >= 0);
-}
+}
/**
* Samples if an event (coalescence, recombination or migration of active nodes)
- * happens in the current TimeInterval or not.
+ * happens in the current TimeInterval or not.
*
* In particular requires that the 'temporary' forest members samples_, rates_
- * and active_nodes_ are set correctly beforehand.
+ * and active_nodes_ are set correctly beforehand.
*
* \param ti The current time interval
* \returns the event that has happened (can also be a "NoEvent" event)
@@ -644,13 +647,13 @@ void Forest::sampleEvent(const TimeInterval &ti, double &event_time, Event &retu
// Sample on which time and time line the event happens (if any)
for (size_t i = 0; i < 3; ++i) {
if (rates_[i] == 0.0) continue;
- selectFirstTime(random_generator()->sampleExpoExpoLimit(rates_[i], getTimeLineGrowth(i), ti.length()),
+ selectFirstTime(random_generator()->sampleExpoExpoLimit(rates_[i], getTimeLineGrowth(i), ti.length()),
i, event_time, event_line);
}
- // Correct the time from relative to the time interval to absolute
+ // Correct the time from relative to the time interval to absolute
if (event_time != -1) event_time += ti.start_height();
- assert( (event_time == -1) ||
+ assert( (event_time == -1) ||
(ti.start_height() <= event_time && event_time <= ti.end_height()) );
assert( (event_line == -1 && event_time == -1) || (event_line != -1 && event_time != -1));
@@ -660,12 +663,12 @@ void Forest::sampleEvent(const TimeInterval &ti, double &event_time, Event &retu
/**
- * Given that an event has happened, this function samples the events type.
- *
- * In particular requires that the 'temporary' forest members samples_, rates_,
- * active_nodes_, and nodes_timelines_ are set correctly beforehand.
+ * Given that an event has happened, this function samples the events type.
+ *
+ * In particular requires that the 'temporary' forest members samples_, rates_,
+ * active_nodes_, and nodes_timelines_ are set correctly beforehand.
*/
-void Forest::sampleEventType(const double time, const size_t time_line,
+void Forest::sampleEventType(const double time, const size_t time_line,
const TimeInterval &ti, Event &event) const {
event = Event(time);
@@ -683,13 +686,13 @@ void Forest::sampleEventType(const double time, const size_t time_line,
// Only Nodes in state 1 or 2 can do something
if ( states_[i] == 0 ) continue;
- // Coalescence can occur on all time lines
+ // Coalescence can occur on all time lines
if (states_[i] == 1 && active_nodes_timelines_[i] == time_line) {
sample -= calcCoalescenceRate(active_node(i)->population(), ti);
- if (sample <= 0.0) return event.setToCoalescence(active_node(i), i);
+ if (sample <= 0.0) return event.setToCoalescence(active_node(i), i);
}
- // Migration and Recombination only on time line 0
+ // Migration and Recombination only on time line 0
if (time_line != 0) continue;
// Recombination
@@ -705,7 +708,7 @@ void Forest::sampleEventType(const double time, const size_t time_line,
for ( size_t j = 0; j < model().population_number(); ++j) {
sample -= model().migration_rate(active_node(i)->population(), j);
if ( sample <= 0.0 ) return event.setToMigration(active_node(i), i, j);
- }
+ }
throw std::logic_error("Error Sampling Type of Event");
}
sample -= model().total_migration_rate(active_node(i)->population());
@@ -714,7 +717,7 @@ void Forest::sampleEventType(const double time, const size_t time_line,
// If we are here, than we should have sampled a pw coalescence...
assert( states_[0] == 1 && states_[1] == 1 );
assert( active_nodes_[0]->population() == active_nodes_[1]->population() );
- assert( sample <= calcPwCoalescenceRate(active_nodes_[0]->population(), ti) );
+ assert( sample <= calcPwCoalescenceRate(active_nodes_[0]->population(), ti) );
return event.setToPwCoalescence();
}
@@ -725,12 +728,12 @@ void Forest::sampleEventType(const double time, const size_t time_line,
* the time_line counter.
*
* \param new_time An ExpoLimit Sample
- * \param time_line The timeline that the sample was from
- * \param current_time The variable that save the time of the nearest event
- * \param time_line The variable that saves the timeline of the nearest event
- * \return Nothing, but updates current_time and current_time_line
+ * \param time_line The timeline that the sample was from
+ * \param current_time The variable that save the time of the nearest event
+ * \param time_line The variable that saves the timeline of the nearest event
+ * \return Nothing, but updates current_time and current_time_line
*/
-void Forest::selectFirstTime(const double new_time, const size_t time_line,
+void Forest::selectFirstTime(const double new_time, const size_t time_line,
double ¤t_time, size_t ¤t_time_line) const {
if (new_time == -1) return;
if (current_time == -1 || new_time < current_time) {
@@ -755,8 +758,8 @@ size_t Forest::getNodeState(Node const *node, const double current_time) const {
if (node->is_root()) return(1);
if (!node->local()) return(2);
dout << "Error getting node state." << std::endl;
- dout << "Height: " << node->height()
- << " current time: " << current_time
+ dout << "Height: " << node->height()
+ << " current time: " << current_time
<< " diff: " << node->height() - current_time << std::endl;
dout << "Node local: " << node->local() << std::endl;
dout << "Node root: " << node->is_root() << std::endl;
@@ -788,8 +791,8 @@ double Forest::calcRecombinationRate(Node const* node) const {
} else {
// Rec rate may change. Accumulate the total rate.
- double rate = model().recombination_rate() *
- (this->current_base() - model().getCurrentSequencePosition());
+ double rate = model().recombination_rate() *
+ (this->current_base() - model().getCurrentSequencePosition());
size_t idx = model().get_position_index() - 1;
while (model().change_position(idx) > last_update_pos) {
@@ -804,7 +807,7 @@ double Forest::calcRecombinationRate(Node const* node) const {
}
-/**
+/**
* Even if no event occurred in a time interval, there is some stuff that we
* might have to do. Mainly moving the "active" flag upwards if a active node
* was looking for recombinations, but none occurred. In that case, we also
@@ -812,11 +815,27 @@ double Forest::calcRecombinationRate(Node const* node) const {
*
* \param ti The current time interval
* \param coalescences_finished temp variable to pass the information that the
- * coalescence has finished.
+ * coalescence has finished.
*/
void Forest::implementNoEvent(const TimeInterval &ti, bool &coalescence_finished) {
- if (ti.end_height() == DBL_MAX)
- throw std::logic_error("Lines did not coalescence. If you use an negative growth parameter (population rapidly declining forward in time), you need to set it to a non-negative value at some time.");
+ if (ti.end_height() == DBL_MAX) {
+ std::stringstream message;
+ message << "Lines did not coalescence." << std::endl;
+ if (active_node(0)->population() != active_node(1)->population()) {
+ message << "The lines were in populations " << active_node(0)->population() + 1
+ << " and " << active_node(1)->population() + 1 << "." << std::endl
+ << "You should add on opportunity for migration between these populations."
+ << std::endl;
+ }
+
+ else if (model().growth_rate(active_node(0)->population())) {
+ message << "Population " << active_node(0)->population() + 1
+ << " has a negative growth factor for infinite time." << std::endl
+ << "This can prevent coalescence. " << std::endl;
+ }
+
+ throw std::logic_error(message.str());
+ }
if (states_[0] == 2) {
set_active_node(0, possiblyMoveUpwards(active_node(0), ti));
if (active_node(0)->local()) {
@@ -825,8 +844,8 @@ void Forest::implementNoEvent(const TimeInterval &ti, bool &coalescence_finished
updateAbove(active_node(0));
coalescence_finished = true;
tmp_event_time_ = active_node(0)->height();
- contemporaries_.replace(active_node(0),
- active_node(0)->first_child(),
+ contemporaries_.replace(active_node(0),
+ active_node(0)->first_child(),
active_node(0)->second_child());
return;
}
@@ -837,7 +856,7 @@ void Forest::implementNoEvent(const TimeInterval &ti, bool &coalescence_finished
if (states_[1] == 2) set_active_node(1, possiblyMoveUpwards(active_node(1), ti));
if (active_node(0) == active_node(1)) {
- dout << "* * * Active Nodes hit each other in " << active_node(0)
+ dout << "* * * Active Nodes hit each other in " << active_node(0)
<< ". Done." << std::endl;
updateAbove(active_node(0));
coalescence_finished = true;
@@ -865,7 +884,7 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
Node* target = contemporaries_.sample(coal_node->population());
dout << "* * * Above node " << target << std::endl;
- assert( target->height() < event.time() );
+ assert( target->height() < event.time() );
assert( coal_node->population() == target->population() );
assert( getEventNode() != NULL );
assert( getOtherNode() != NULL );
@@ -876,7 +895,7 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
// Actually implement the coalescence
// ---------------------------------------------
- // Look if we can reuse the root that coalesced for marking the point of
+ // Look if we can reuse the root that coalesced for marking the point of
// coalescence
if ( coal_node->countChildren() == 1 && !coal_node->is_migrating() ){
assert( coal_node->local() );
@@ -894,7 +913,7 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
// Now we have:
// target: Node in the target tree under the coalescence
- // coal_node: Root of the coalescing tree
+ // coal_node: Root of the coalescing tree
// new_node: New parent of 'target' and 'coal_node'
// Update new_node
@@ -914,7 +933,7 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
// And update
coal_node->make_local();
- updateAbove(coal_node, false, false);
+ updateAbove(coal_node, false, false);
set_active_node(event.active_node_nr(), new_node);
@@ -924,7 +943,7 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
// ---------------------------------------------
if ( getOtherNodesState() == 2 ) {
- // If the coalescing node coalesced into the branch directly above
+ // If the coalescing node coalesced into the branch directly above
// a recombining node, we are done.
if ( getOtherNode()->parent() == getEventNode() ) {
dout << "* * * Recombining Node moved into coalesced node. Done." << std::endl;
@@ -940,15 +959,15 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
}
if ( target->local() ) {
- // Only active_node(0) can coalescence into local nodes. active_node(1) is
- // at least the local root and hence above all local nodes.
+ // Only active_node(0) can coalescence into local nodes. active_node(1) is
+ // at least the local root and hence above all local nodes.
// If active_node(0) coalescences into the local tree, there are no more
- // active nodes and we are done.
+ // active nodes and we are done.
assert( event.active_node_nr() == 0 );
assert( states_[1] == 0 );
dout << "* * * We hit the local tree. Done." << std::endl;
- updateAbove(getEventNode());
+ updateAbove(getEventNode());
coalescence_finished_ = true;
contemporaries_.replace(new_node, coal_node, target);
return;
@@ -962,9 +981,9 @@ void Forest::implementCoalescence(const Event &event, TimeIntervalIterator &tii)
-/**
+/**
* @brief Modifies the forest to reflect that two coalescing nodes coalesced together.
- *
+ *
* @param root_1 The first coalescing node
* @param root_2 The second coalescing node
* @param time The time at which the coalescence happens
@@ -998,7 +1017,7 @@ void Forest::implementPwCoalescence(Node* root_1, Node* root_2, const double tim
new_root = root_1;
root_1 = root_1->first_child();
assert( root_1 != NULL );
- }
+ }
else if (root_2->countChildren() == 1 && !root_2->is_migrating()) {
// only root_2 has a single branch => use as new root
this->nodes()->move(root_2, time);
@@ -1044,8 +1063,8 @@ void Forest::implementMigration(const Event &event, const bool &recalculate, Tim
dout << "* * Implementing migration event... " << std::flush;
assert( event.node()->is_root() );
- // There is only little to do if we can reuse the event.node()
- if ( event.node()->is_unimportant() ||
+ // There is only little to do if we can reuse the event.node()
+ if ( event.node()->is_unimportant() ||
( event.node()->height() == event.time() && event.node()->is_migrating() ) ) {
dout << "Reusing: " << event.node() << "... " << std::flush;
nodes()->move(event.node(), event.time());
@@ -1054,7 +1073,7 @@ void Forest::implementMigration(const Event &event, const bool &recalculate, Tim
} else {
// Otherwise create a new node that marks the migration event,
Node* mig_node = nodes()->createNode(event.time());
- dout << "Marker: " << mig_node << "... " << std::flush;
+ dout << "Marker: " << mig_node << "... " << std::flush;
nodes()->add(mig_node, event.node());
mig_node->set_population(event.mig_pop());
@@ -1104,10 +1123,13 @@ void Forest::implementFixedTimeEvent(TimeIntervalIterator &ti) {
}
// Stop if no migration occurred
- if (!migrated) break;
+ 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)
+ if (chain_cnt == 10000)
throw std::logic_error("Circle detected when moving individuals between populations");
++chain_cnt;
}
@@ -1115,7 +1137,7 @@ void Forest::implementFixedTimeEvent(TimeIntervalIterator &ti) {
assert( printTree() );
}
-/**
+/**
* Helper function for doing a coalescence.
* Moves the 'active' flag (i.e. the node stored in root_1 or root_2 in sampleCoalescence)
* from a node to it's parent if the branch above the node
@@ -1130,7 +1152,7 @@ void Forest::implementFixedTimeEvent(TimeIntervalIterator &ti) {
*
* \param node An active node
* \param time_interval The time interval the coalescence is currently in.
- *
+ *
* \return Either the parent of 'node' if we need to move upwards or 'node'
* itself
*/
@@ -1159,7 +1181,7 @@ bool Forest::pruneNodeIfNeeded(Node* node, const bool prune_orphans) {
// of the current coalescence.
if (node == primary_root()) set_primary_root(NULL);
nodes()->remove(node);
- return true;
+ return true;
}
// Other roots stay
return false;
@@ -1171,20 +1193,20 @@ bool Forest::pruneNodeIfNeeded(Node* node, const bool prune_orphans) {
assert(!node->is_root());
node->parent()->change_child(node, NULL);
- if (node->countChildren() == 0) nodes()->remove(node);
- else {
+ if (node->countChildren() == 0) nodes()->remove(node);
+ else {
// Remove link between `node` and its parent,
- // which effectively removes the branch,
+ // which effectively removes the branch,
// and separates the subtree below from the main tree
Node* parent = node->parent();
node->set_parent(NULL);
updateAbove(parent, false, true, true);
}
return true;
- }
+ }
else if (node->countChildren() == 1 && !node->is_migrating()) {
- // Unneeded nodes
+ // Unneeded nodes
dout << "* * * PRUNING: Removing node " << node << " from tree (unneeded)" << std::endl;
assert(!node->is_migrating());
assert(node->first_child()->last_update() == node->last_update());
@@ -1192,7 +1214,7 @@ bool Forest::pruneNodeIfNeeded(Node* node, const bool prune_orphans) {
Node* child = node->first_child();
child->set_parent(node->parent());
- node->parent()->change_child(node, child);
+ node->parent()->change_child(node, child);
nodes()->remove(node);
return true;
}
@@ -1245,7 +1267,7 @@ void Forest::clear() {
// Clear Summary Statistics
this->clearSumStats();
-
+
// Reset Model
writable_model()->resetTime();
writable_model()->resetSequencePosition();
@@ -1340,5 +1362,8 @@ void Forest::sampleNextBase() {
// Recombination in the sequence segment
set_next_base(current_base() + length);
}
-}
+
+ assert(next_base() > current_base());
+ assert(next_base() <= model().loci_length());
+}
diff --git a/src/forest.h b/src/forest.h
index 7dec984..db2a867 100644
--- a/src/forest.h
+++ b/src/forest.h
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -95,9 +95,9 @@ class Forest
Node* primary_root() const { return primary_root_; }
void set_primary_root(Node* primary_root) { primary_root_ = primary_root; };
- size_t sample_size() const {
- if (sample_size_ == 0) return model().sample_size();
- return this->sample_size_;
+ size_t sample_size() const {
+ if (sample_size_ == 0) return model().sample_size();
+ return this->sample_size_;
}
void set_sample_size(const size_t size ) { sample_size_ = size; }
@@ -111,7 +111,7 @@ class Forest
*
* @param finite_sites If 'true', the length is measured in number of bases
* (e.g. integer sequence positions) for which the tree is valid. Otherwise,
- * the length is measured real valued number on a continuous chromosome.
+ * the length is measured real valued number on a continuous chromosome.
*
* @return The length of the current segment (see above for its unit)
*/
@@ -124,7 +124,7 @@ class Forest
}
void set_random_generator(RandomGenerator *rg) {
- this->random_generator_ = rg;
+ this->random_generator_ = rg;
}
RandomGenerator* random_generator() const { return this->random_generator_; }
@@ -175,15 +175,15 @@ class Forest
}
//derived class from Forest
- virtual void record_Recombevent(size_t pop_i,
- double opportunity,
+ virtual void record_Recombevent(size_t pop_i,
+ double opportunity,
eventCode event_code){
(void)pop_i;
(void)opportunity;
- (void)event_code;
+ (void)event_code;
}
virtual void record_all_event( TimeInterval const &ti){
- (void) ti;
+ (void) ti;
}
// Calc & Print Summary Statistics
@@ -193,24 +193,25 @@ class Forest
void printSegmentSumStats(std::ostream &output) const;
double get_rec_base(const size_t idx) const {
- assert(idx < rec_bases_.size());
- return rec_bases_[idx];
+ return rec_bases_.at(idx);
}
double current_base() const { return get_rec_base(current_rec_); }
double next_base() const { return get_rec_base(current_rec_ + 1); }
- void set_current_base(double const base) { rec_bases_[current_rec_] = base; };
- void set_next_base(double const base) { rec_bases_.push_back(base); };
+ void set_current_base(double const base) { rec_bases_[current_rec_] = base; };
+ void set_next_base(double const base) { rec_bases_.push_back(base); };
size_t current_rec() const { return current_rec_; };
+ bool coalescence_finished() const { return this->coalescence_finished_; }
+
private:
Forest() { this->initialize(); }
//Operations on the Tree
Node* cut(const TreePoint &cut_point);
- void updateAbove(Node* node,
+ void updateAbove(Node* node,
bool above_local_root = false,
const bool &recursive = true,
const bool &invariants_only = false);
@@ -218,7 +219,7 @@ class Forest
// Tools for doing coalescence & recombination
void sampleCoalescences(Node *start_node);
size_t getNodeState(Node const *node, const double current_time) const;
- Node* updateBranchBelowEvent(Node* node, const TreePoint &event_point);
+ Node* updateBranchBelowEvent(Node* node, const TreePoint &event_point);
Node* possiblyMoveUpwards(Node* node, const TimeInterval &event);
// Implementation of the different events
@@ -233,11 +234,11 @@ class Forest
bool nodeIsOld(Node const* node) const {
if ( node->local() ) return false;
if ( node->is_root() ) return false;
- if ( model().has_window_rec() &&
+ if ( model().has_window_rec() &&
segment_count() - node->last_update() > model().window_length_rec()) {
return true;
}
- if ( model().has_window_seq() &&
+ if ( model().has_window_seq() &&
current_base() - get_rec_base(node->last_update()) > model().window_length_seq()) {
return true;
}
@@ -258,10 +259,10 @@ class Forest
void sampleEvent(const TimeInterval &ti, double &event_time, Event &return_event) const;
- void sampleEventType(const double time, const size_t time_line,
+ void sampleEventType(const double time, const size_t time_line,
const TimeInterval &ti, Event &return_event) const;
- void selectFirstTime(const double new_time, const size_t time_line,
+ void selectFirstTime(const double new_time, const size_t time_line,
double ¤t_time, size_t ¤t_time_line) const;
double getTimeLineGrowth(const size_t time_line) const {
@@ -271,6 +272,7 @@ class Forest
else throw std::out_of_range("Trying to get growthrate of unknown time line.");
}
+
// Private Members
NodeContainer nodes_; // The nodes of the Tree/Forest
@@ -287,7 +289,7 @@ class Forest
size_t sample_size_; // The number of sampled nodes (changes while building the initial tree)
size_t current_rec_; // A counter for recombinations
- std::vector<double> rec_bases_; // Genetic positions of the recombinations
+ std::vector<double> rec_bases_; // Genetic positions of the recombinations
Model* model_;
RandomGenerator* random_generator_;
@@ -299,7 +301,7 @@ class Forest
// Temporarily used when doing the coalescence
- // Rates:
+ // Rates:
double rates_[3];
// States: Each (branch above an) active node can either be in state
diff --git a/src/model.cc b/src/model.cc
index 815e3bb..aa3127c 100644
--- a/src/model.cc
+++ b/src/model.cc
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
@@ -22,9 +22,8 @@
#include "model.h"
-Model::Model() {
+Model::Model() {
default_pop_size = 10000;
- default_loci_length = 100000;
default_growth_rate = 0.0;
default_mig_rate = 0.0;
scaling_factor_ = 1.0 / (4 * default_pop_size);
@@ -32,15 +31,13 @@ Model::Model() {
has_migration_ = false;
has_recombination_ = false;
+ this->set_loci_number(1);
+ this->setLocusLength(1);
this->addChangeTime(0.0);
this->addChangePosition(0.0);
this->set_population_number(1);
- this->set_loci_number(1);
- this->loci_length_ = this->default_loci_length;
-
- this->setLocusLength(default_loci_length);
this->setMutationRate(0.0);
this->setRecombinationRate(0.0);
@@ -56,6 +53,7 @@ Model::Model() {
Model::Model(size_t sample_size) : Model() {
this->addSampleSizes(0.0, std::vector<size_t>(1, sample_size));
+ this->setLocusLength(1000);
this->resetTime();
}
@@ -71,7 +69,7 @@ void Model::reset() {
}
*/
-/**
+/**
* Function to add a new change time to the model.
*
* It preserves the relation between the times and the *param*_list_ containers.
@@ -80,7 +78,7 @@ void Model::reset() {
*
* @param time The time that is added
* @param scaled set to TRUE if the time is in units of 4N0 generations, and
- * FALSE if it is in units of generations.
+ * FALSE if it is in units of generations.
*
* @returns The position the time has now in the vector
*/
@@ -98,16 +96,16 @@ size_t Model::addChangeTime(double time, const bool &scaled) {
return position;
}
- std::vector<double>::iterator ti;
+ std::vector<double>::iterator ti;
for (ti = change_times_.begin(); ti != change_times_.end(); ++ti) {
if ( *ti == time ) return position;
- if ( *ti > time ) break;
+ if ( *ti > time ) break;
++position;
}
change_times_.insert(ti, time);
-
- // Add Null at the right position in all parameter vectors
+
+ // Add Null at the right position in all parameter vectors
pop_sizes_list_.insert(pop_sizes_list_.begin() + position, std::vector<double>());
growth_rates_list_.insert(growth_rates_list_.begin() + position, std::vector<double>());
mig_rates_list_.insert(mig_rates_list_.begin() + position, std::vector<double>());
@@ -117,21 +115,23 @@ size_t Model::addChangeTime(double time, const bool &scaled) {
}
-/**
- * Function to add a new change time to the model.
+/**
+ * Adds a new change position to the model.
*
- * It preserves the relation between the times and the *param*_list_ containers.
- * If the same time is added multiple times, it is just added once to the model,
- * but this should not make a difference when using this function.
+ * Change position are sequence positions where mutation or recombination rates
+ * change. This creates a new position, but does not add the new rates.
*
- * @param time The time that is added
- * @param scaled set to TRUE if the time is in units of 4N0 generations, and
- * FALSE if it is in units of generations.
+ * @param position The sequence position add which a change is added
*
- * @returns The position the time has now in the vector
+ * @returns The index of the new rates in the recombination_rates_ and
+ * mutation_rates vectors.
*/
size_t Model::addChangePosition(const double position) {
- size_t idx = 0;
+ if (position < 0 || position > loci_length()) {
+ std::stringstream ss;
+ ss << "Rate change position " << position << " is outside of the simulated sequence.";
+ throw std::invalid_argument(ss.str());
+ }
if ( change_position_.size() == 0 ) {
change_position_ = std::vector<double>(1, position);
@@ -140,16 +140,17 @@ size_t Model::addChangePosition(const double position) {
return 0;
}
- std::vector<double>::iterator ti;
+ size_t idx = 0;
+ std::vector<double>::iterator ti;
for (ti = change_position_.begin(); ti != change_position_.end(); ++ti) {
if ( *ti == position ) return idx;
if ( *ti > position ) break;
- ++idx;
+ ++idx;
}
change_position_.insert(ti, position);
-
- // Add Null at the right position in all parameter vectors
+
+ // Add Null at the right position in all parameter vectors
recombination_rates_.insert(recombination_rates_.begin() + idx, -1);
mutation_rates_.insert(mutation_rates_.begin() + idx, -1);
@@ -171,9 +172,9 @@ void Model::addSampleSizes(double time, const std::vector<size_t> &samples_sizes
/**
* @brief This changes the size of all populations to the given values at a
- * specific point in time.
+ * specific point in time.
*
- * The sizes apply for with point on backwards in time.
+ * The sizes apply for with point on backwards in time.
*
* @param time The time at which the population changes their sizes.
* @param pop_sizes A vector of new population sizes. Can either be given as
@@ -183,10 +184,10 @@ void Model::addSampleSizes(double time, const std::vector<size_t> &samples_sizes
* @param relative set to TRUE, if the population sizes are given relative to
* N0, or to FALSE if they are absolute values.
*/
-void Model::addPopulationSizes(double time, const std::vector<double> &pop_sizes,
+void Model::addPopulationSizes(double time, const std::vector<double> &pop_sizes,
const bool &time_scaled, const bool &relative) {
- if ( pop_sizes.size() != population_number() )
+ if ( pop_sizes.size() != population_number() )
throw std::logic_error("Population size values do not meet the number of populations");
size_t position = addChangeTime(time, time_scaled);
@@ -199,8 +200,8 @@ void Model::addPopulationSizes(double time, const std::vector<double> &pop_sizes
// Save inverse double value
if (pop_size <= 0.0) throw std::invalid_argument("population size <= 0");
- pop_size = 1.0 / (2 * pop_size);
- }
+ pop_size = 1.0 / (2 * pop_size);
+ }
pop_sizes_list_[position].push_back(pop_size);
}
}
@@ -208,9 +209,9 @@ void Model::addPopulationSizes(double time, const std::vector<double> &pop_sizes
/**
* @brief This changes the size of all populations to a given value at a
- * specific point in time.
+ * specific point in time.
*
- * The sizes apply for with point on backwards in time.
+ * The sizes apply for with point on backwards in time.
*
* @param time The time at which the population changes their sizes.
* @param pop_sizes The size to which we set all populations. Can either be given as
@@ -220,17 +221,17 @@ void Model::addPopulationSizes(double time, const std::vector<double> &pop_sizes
* @param relative set to TRUE, if the population sizes are given relative to
* N0, or to FALSE if they are absolute values.
*/
-void Model::addPopulationSizes(const double time, const double pop_size,
+void Model::addPopulationSizes(const double time, const double pop_size,
const bool &time_scaled, const bool &relative) {
addPopulationSizes(time, std::vector<double>(population_number(), pop_size), time_scaled, relative);
}
-
+
/**
* @brief This changes the size of a single populations to a given value at a
- * specific point in time.
+ * specific point in time.
*
- * The sizes apply for with point on backwards in time.
+ * The sizes apply for with point on backwards in time.
* Requires Model.finalization() to be called after the model is set up.
*
* @param time The time at which the population change its size.
@@ -258,7 +259,7 @@ void Model::addPopulationSize(const double time, const size_t pop, double popula
* @brief Set the population size growth rates at a certain time point.
*
* The population growth or shrinks exponentially from that time point on
- * backwards in time.
+ * backwards in time.
* Requires Model.finalization() to be called after the model is set up.
*
* @param time The time at which to set the growth rates
@@ -268,7 +269,7 @@ void Model::addPopulationSize(const double time, const size_t pop, double popula
*/
void Model::addGrowthRates(const double time, const std::vector<double> &growth_rates,
const bool &time_scaled, const bool &rate_scaled) {
- if ( growth_rates.size() != population_number() )
+ if ( growth_rates.size() != population_number() )
throw std::logic_error("Growth rates values do not meet the number of populations");
size_t position = addChangeTime(time, time_scaled);
@@ -284,7 +285,7 @@ void Model::addGrowthRates(const double time, const std::vector<double> &growth_
* @brief Set the population size growth rates at a certain time point.
*
* The population growth or shrinks exponentially from that time point on
- * backwards in time.
+ * backwards in time.
* Requires Model.finalization() to be called after the model is set up.
*
* @param time The time at which to set the growth rates
@@ -302,7 +303,7 @@ void Model::addGrowthRates(const double time, const double growth_rate,
* @brief Set the population size growth rates of a population at a certain time point.
*
* The population growth or shrinks exponentially from that time point on
- * backwards in time.
+ * backwards in time.
* Requires Model.finalization() to be called after the model is set up.
*
* @param time The time at which to set the growth rates
@@ -311,19 +312,19 @@ void Model::addGrowthRates(const double time, const double growth_rate,
* @param time_scaled Set to true if the time is given in units of 4*N0
* generations, or to false if the time is given in units of generations.
*/
-void Model::addGrowthRate(const double time, const size_t population,
+void Model::addGrowthRate(const double time, const size_t population,
double growth_rate, const bool &time_scaled, const bool &rate_scaled) {
checkPopulation(population);
size_t position = addChangeTime(time, time_scaled);
if (rate_scaled) growth_rate *= scaling_factor();
- if (growth_rates_list_.at(position).empty()) addGrowthRates(time, nan("number to replace"), time_scaled);
+ if (growth_rates_list_.at(position).empty()) addGrowthRates(time, nan("number to replace"), time_scaled);
growth_rates_list_.at(position).at(population) = growth_rate;
}
/**
* @brief Sets a migration rate form a specific population to another starting from a
- * certain time point (going backwards in time);
+ * certain time point (going backwards in time);
*
* This requires model finalization, e.g. call model.finalize() after you set up
* the model completely.
@@ -334,8 +335,8 @@ void Model::addGrowthRate(const double time, const size_t population,
* looking backwards in time. Is the sink population when looking forward.
* @param sink The population to which the individuals migrate to (also
* when looking backwards in time)
- * @param mig_rate The backwards scaled migration rate M_ij = 4N0 * m_ij,
- * where m_ij is the fraction for population i = source that migrates
+ * @param mig_rate The backwards scaled migration rate M_ij = 4N0 * m_ij,
+ * where m_ij is the fraction for population i = source that migrates
* to population j = sink (again, when looking backwards in time).
* @param scaled_time Set to true if the time is given in units of 4*N0
* generations, or to false if the time is given in units of generations.
@@ -349,16 +350,16 @@ void Model::addMigrationRate(double time, size_t source, size_t sink, double mig
checkPopulation(sink);
size_t position = addChangeTime(time, scaled_time);
if (scaled_rates) mig_rate *= scaling_factor();
- if (mig_rates_list_.at(position).empty()) {
- addSymmetricMigration(time, nan("value to replace"), scaled_time);
+ if (mig_rates_list_.at(position).empty()) {
+ addSymmetricMigration(time, nan("value to replace"), scaled_time);
}
- mig_rates_list_.at(position).at(getMigMatrixIndex(source, sink)) = mig_rate;
+ mig_rates_list_.at(position).at(getMigMatrixIndex(source, sink)) = mig_rate;
}
-/**
+/**
* @brief Sets the migration matrix to the given values for the time following at
- * certain time point (backwards in time).
+ * certain time point (backwards in time).
*
* This requires model finalization, e.g. call model.finalize() after you set up
* the model completely.
@@ -370,7 +371,7 @@ void Model::addMigrationRate(double time, size_t source, size_t sink, double mig
* M_ij = 4N0 * m_ij and m_ij is the faction of population i that
* migrates to population j (viewed backwards in time; forwards the
* migration is from population j to i). The diagonal elements of the
- * matrix are ignored and can be set to "x" for better readability.
+ * matrix are ignored and can be set to "x" for better readability.
* @param scaled_time Set to true if the time is given in units of 4*N0
* generations, or to false if the time is given in units of generations.
* @param scaled_rate Set to true if the rate is given as M = 4*N0*m and to
@@ -381,7 +382,7 @@ void Model::addMigrationRates(double time, const std::vector<double> &mig_rates,
double popnr = population_number();
double scaling = 1;
if (scaled_rates) scaling = scaling_factor();
- if ( mig_rates.size() != population_number()*population_number() )
+ if ( mig_rates.size() != population_number()*population_number() )
throw std::logic_error("Migration rates values do not meet the number of populations");
size_t position = addChangeTime(time, scaled_time);
@@ -390,7 +391,7 @@ void Model::addMigrationRates(double time, const std::vector<double> &mig_rates,
for (size_t i = 0; i < popnr; ++i) {
for (size_t j = 0; j < popnr; ++j) {
if (i == j) continue;
- mig_rates_list_[position].push_back(mig_rates.at(i*popnr+j) * scaling);
+ mig_rates_list_[position].push_back(mig_rates.at(i*popnr+j) * scaling);
}
}
}
@@ -410,21 +411,21 @@ void Model::addMigrationRates(double time, const std::vector<double> &mig_rates,
* @param time The time at which the migration is set to the given value.
* It applies backwards in time until it is changed again.
* @param mig_rate The scaled migration rate M_ij = 4N0 * m_ij that is used
- * between all populations i and j. m_ij is the fraction of population i
+ * between all populations i and j. m_ij is the fraction of population i
* that migrates to population j.
* @param time_scaled Set to true if the time is given in units of 4*N0
* generations, or to false if the time is given in units of generations.
* @param rate_scaled Set to true if the rate is given as M = 4*N0*m and to
* false if it is given as m.
*/
-void Model::addSymmetricMigration(const double time, const double mig_rate,
+void Model::addSymmetricMigration(const double time, const double mig_rate,
const bool &time_scaled, const bool &rate_scaled) {
std::vector<double> mig_rates = std::vector<double>(population_number()*population_number(), mig_rate);
this->addMigrationRates(time, mig_rates, time_scaled, rate_scaled);
}
-void Model::addSingleMigrationEvent(const double time, const size_t source_pop,
+void Model::addSingleMigrationEvent(const double time, const size_t source_pop,
const size_t sink_pop, const double fraction,
const bool &time_scaled) {
@@ -440,37 +441,37 @@ void Model::addSingleMigrationEvent(const double time, const size_t source_pop,
single_mig_probs_list_.at(position) = std::vector<double>(popnr*popnr-popnr, 0.0);
}
- single_mig_probs_list_.at(position).at(getMigMatrixIndex(source_pop, sink_pop)) = fraction;
+ single_mig_probs_list_.at(position).at(getMigMatrixIndex(source_pop, sink_pop)) = fraction;
this->has_migration_ = true;
-}
+}
std::ostream& operator<<(std::ostream& os, Model& model) {
size_t n_pops = model.population_number();
os << "---- Model: ------------------------" << std::endl;
- os << "Total Sample Size: " << model.sample_size() << std::endl;
- os << "N0 is assumed to be " << model.default_pop_size << std::endl;
+ os << "Total Sample Size: " << model.sample_size() << std::endl;
+ os << "N0 is assumed to be " << model.default_pop_size << std::endl;
model.resetSequencePosition();
for (size_t idx = 0; idx < model.countChangePositions(); ++idx) {
- os << std::endl << "At Position " << model.getCurrentSequencePosition() << ":" << std::endl;
- os << " Mutation Rate: " << model.mutation_rate() << std::endl;
- os << " Recombination Rate: " << model.recombination_rate() << std::endl;
+ os << std::endl << "At Position " << model.getCurrentSequencePosition() << ":" << std::endl;
+ os << " Mutation Rate: " << model.mutation_rate() << std::endl;
+ os << " Recombination Rate: " << model.recombination_rate() << std::endl;
model.increaseSequencePosition();
}
model.resetSequencePosition();
-
+
model.resetTime();
- for (size_t idx = 0; idx < model.countChangeTimes(); ++idx) {
- os << std::endl << "At Time " << model.getCurrentTime() << ":" << std::endl;
-
+ for (size_t idx = 0; idx < model.countChangeTimes(); ++idx) {
+ os << std::endl << "At Time " << model.getCurrentTime() << ":" << std::endl;
+
os << " Population Sizes: ";
- for (size_t pop = 0; pop < n_pops; ++pop)
+ for (size_t pop = 0; pop < n_pops; ++pop)
os << std::setw(10) << std::right << model.population_size(pop, model.getCurrentTime());
os << std::endl;
os << " Growth Rates: ";
- for (size_t pop = 0; pop < n_pops; ++pop)
+ for (size_t pop = 0; pop < n_pops; ++pop)
os << std::setw(10) << std::right << model.growth_rate(pop);
os << std::endl;
@@ -481,11 +482,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 "
+ os << " " << model.single_mig_pop(i, j) * 100 << "% of pop "
<< i + 1 << " move to pop " << j + 1 << std::endl;
}
}
@@ -504,7 +505,7 @@ void Model::updateTotalMigRates(const size_t position) {
total_mig_rates_list_.at(position) = std::vector<double>(population_number(), 0.0);
}
- std::vector<double>* mig_rates = &(total_mig_rates_list_.at(position));
+ std::vector<double>* mig_rates = &(total_mig_rates_list_.at(position));
for (size_t i = 0; i < population_number(); ++i) {
for (size_t j = 0; j < population_number(); ++j) {
@@ -524,7 +525,7 @@ void Model::finalize() {
for (size_t j = 0; j < mig_rates_list_.size(); ++j) {
if (mig_rates_list_.at(j).empty()) continue;
updateTotalMigRates(j);
- }
+ }
// Fill in missing recombination rates
assert( mutation_rates_.at(0) != -1 );
@@ -547,7 +548,7 @@ void Model::finalize() {
void Model::calcPopSizes() {
// Set initial population sizes
- if (pop_sizes_list_.at(0).empty()) addPopulationSizes(0, default_pop_size);
+ if (pop_sizes_list_.at(0).empty()) addPopulationSizes(0, default_pop_size);
else {
// Replace values not set by the user with the default size
for (size_t pop = 0; pop < population_number(); ++pop) {
@@ -564,28 +565,28 @@ void Model::calcPopSizes() {
// Make sure we always have a pop sizes vector
if (pop_sizes_list_.at(i).empty()) {
addPopulationSizes(change_times_.at(i), nan("value to replace"));
- assert( ! pop_sizes_list_.at(i).empty() );
+ assert( ! pop_sizes_list_.at(i).empty() );
}
- // Calculate the effective duration of a growth period if it ends here
+ // Calculate the effective duration of a growth period if it ends here
duration = change_times_.at(i) - change_times_.at(i - 1);
- // Calculate the population sizes:
+ // Calculate the population sizes:
for (size_t pop = 0; pop < population_number(); ++pop) {
// If the user explicitly gave a value => always use this value
- if ( !std::isnan(pop_sizes_list_.at(i).at(pop)) ) continue;
-
- assert(!std::isnan(pop_sizes_list_.at(i - 1).at(pop)));
+ if ( !std::isnan(pop_sizes_list_.at(i).at(pop)) ) continue;
+
+ assert(!std::isnan(pop_sizes_list_.at(i - 1).at(pop)));
// Otherwise use last value
- pop_sizes_list_.at(i).at(pop) = pop_sizes_list_.at(i - 1).at(pop);
+ pop_sizes_list_.at(i).at(pop) = pop_sizes_list_.at(i - 1).at(pop);
- // ... and scale it if there was growth
+ // ... and scale it if there was growth
if (last_growth != -1) {
- pop_sizes_list_.at(i).at(pop) *=
+ pop_sizes_list_.at(i).at(pop) *=
std::exp((growth_rates_list_.at(last_growth).at(pop)) * duration);
}
}
- }
+ }
}
@@ -595,13 +596,13 @@ void Model::check() {
// Structure without migration?
if (population_number() > 1 && !has_migration())
- throw std::invalid_argument("Model has multiple populations but no migration. Coalescence impossible");
+ throw std::invalid_argument("Model has multiple populations but no migration. Coalescence impossible");
}
void Model::fillVectorList(std::vector<std::vector<double> > &vector_list, const double default_value) {
- std::vector<double>* last = NULL;
- std::vector<double>* current = NULL;
+ std::vector<double>* last = NULL;
+ std::vector<double>* current = NULL;
for (size_t j = 0; j < vector_list.size(); ++j) {
current = &(vector_list.at(j));
if (current->empty()) continue;
@@ -610,7 +611,7 @@ void Model::fillVectorList(std::vector<std::vector<double> > &vector_list, const
if ( !std::isnan(current->at(i)) ) continue;
if (last == NULL) current->at(i) = default_value;
- else current->at(i) = last->at(i);
+ else current->at(i) = last->at(i);
}
last = current;
}
@@ -656,9 +657,9 @@ void Model::addPopToVectorList(std::vector<std::vector<double> > &vector_list) {
/**
* @brief Sets the recombination rate
*
- * @param rate The recombination rate per sequence length per time unit.
+ * @param rate The recombination rate per sequence length per time unit.
* The sequence length can be either per locus or per base pair and the time
- * can be given in units of 4N0 generations ("scaled") or per generation.
+ * can be given in units of 4N0 generations ("scaled") or per generation.
*
* @param loci_length The length of every loci.
* @param per_locus Set to TRUE, if the rate is given per_locus, and to FALSE
@@ -666,16 +667,16 @@ void Model::addPopToVectorList(std::vector<std::vector<double> > &vector_list) {
* @param scaled Set to TRUE is the rate is scaled with 4N0, or to FALSE if
* it isn't
*/
-void Model::setRecombinationRate(double rate,
- const bool &per_locus,
+void Model::setRecombinationRate(double rate,
+ const bool &per_locus,
const bool &scaled,
- const double seq_position) {
+ const double seq_position) {
- if (rate < 0.0) {
+ if (rate < 0.0) {
throw std::invalid_argument("Recombination rate must be non-negative");
}
- if (scaled) rate /= 4.0 * default_pop_size;
+ if (scaled) rate /= 4.0 * default_pop_size;
if (per_locus) {
if (loci_length() <= 1) {
throw std::invalid_argument("Locus length must be at least two");
@@ -684,7 +685,7 @@ void Model::setRecombinationRate(double rate,
}
if (rate > 0.0) has_recombination_ = true;
- recombination_rates_[addChangePosition(seq_position)] = rate;
+ recombination_rates_[addChangePosition(seq_position)] = rate;
}
@@ -697,9 +698,9 @@ void Model::setRecombinationRate(double rate,
* @param scaled Set this to TRUE if you give the mutation rate in scaled
* units (e.g. as theta rather than as u).
*/
-void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled,
- const double seq_position) {
- if (scaled) rate /= 4.0 * default_pop_size;
+void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled,
+ const double seq_position) {
+ if (scaled) rate /= 4.0 * default_pop_size;
size_t idx = addChangePosition(seq_position);
if (per_locus) {
@@ -708,3 +709,4 @@ void Model::setMutationRate(double rate, const bool &per_locus, const bool &scal
mutation_rates_.at(idx) = rate;
}
}
+
diff --git a/src/model.h b/src/model.h
index f2f9e15..b358732 100644
--- a/src/model.h
+++ b/src/model.h
@@ -44,6 +44,7 @@
#include <cmath>
#include <iomanip>
#include <memory>
+#include <sstream>
#include "summary_statistics/summary_statistic.h"
@@ -69,7 +70,6 @@ class Model
// Default values;
double default_pop_size;
- size_t default_loci_length;
double default_growth_rate;
double default_mig_rate;
@@ -402,9 +402,7 @@ class Model
SeqScale getSequenceScaling() const { return seq_scale_; }
void setSequenceScaling(SeqScale seq_scale) { seq_scale_ = seq_scale; };
-
- private:
- std::vector<double> change_times_;
+
void setLocusLength(const size_t length) {
// Rescale the rates that are per base pair
for (size_t i = 0; i < change_position_.size(); ++i) {
@@ -414,6 +412,9 @@ class Model
loci_length_ = length;
}
+ private:
+ std::vector<double> change_times_;
+
double change_position(size_t idx) const {
return this->change_position_.at(idx);
}
diff --git a/src/node_container.cc b/src/node_container.cc
index 7d996df..a8510ae 100644
--- a/src/node_container.cc
+++ b/src/node_container.cc
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -33,7 +33,7 @@ NodeContainer::NodeContainer() {
size_ = 0;
node_counter_ = 0;
- lane_counter_ = 0;
+ lane_counter_ = 0;
std::vector<Node>* new_lane = new std::vector<Node>();
new_lane->reserve(10000);
node_lanes_.push_back(new_lane);
@@ -43,17 +43,23 @@ NodeContainer::NodeContainer(const NodeContainer &nc) {
size_ = 0;
set_first(NULL);
set_last(NULL);
+
+ node_counter_ = 0;
+ lane_counter_ = 0;
+ std::vector<Node>* new_lane = new std::vector<Node>();
+ new_lane->reserve(10000);
+ node_lanes_.push_back(new_lane);
+
this->unsorted_node_ = NULL;
- this->last_node_ = NULL;
- this->first_node_ = NULL;
- std::map<Node const*, Node*> node_mapping;
+ std::map<Node const*, Node*> node_mapping;
node_mapping[NULL] = NULL;
for (auto it = nc.iterator(); it.good(); ++it) {
- Node *node = new Node(**it);
+ Node *node = this->createNode(**it);
+ //Node *node = new Node(**it);
node_mapping[*it] = node;
- add(node);
+ this->add(node);
}
assert( this->sorted() );
@@ -62,7 +68,6 @@ NodeContainer::NodeContainer(const NodeContainer &nc) {
(*it)->set_first_child(node_mapping[(*it)->first_child()]);
(*it)->set_second_child(node_mapping[(*it)->second_child()]);
}
-
unsorted_node_ = node_mapping[nc.unsorted_node_];
}
@@ -78,7 +83,7 @@ Node* NodeContainer::at(size_t nr) const {
current = current->next();
}
- if ( current == NULL ) throw std::out_of_range("NodeContainer out of range");
+ if ( current == NULL ) throw std::out_of_range("NodeContainer out of range");
return current;
}
@@ -162,14 +167,14 @@ void NodeContainer::add(Node* node, Node* after_node) {
while ( current->height() <= node->height() ) {
assert( !current->is_last() );
if ( !current->is_root() ) {
- if ( current->parent_height() < node->height() ) {
+ if ( current->parent_height() < node->height() ) {
current = current->parent();
continue;
}
}
current = current->next();
}
-
+
// And add the node;
this->add_before(node, current);
assert( this->sorted() );
@@ -177,7 +182,7 @@ void NodeContainer::add(Node* node, Node* after_node) {
void NodeContainer::remove(Node *node, const bool &del) {
- --size_;
+ --size_;
if ( node->is_first() && node->is_last() ) {
this->set_first(NULL);
this->set_last(NULL);
@@ -194,7 +199,7 @@ void NodeContainer::remove(Node *node, const bool &del) {
node->previous()->set_next(node->next());
node->next()->set_previous(node->previous());
}
-
+
if (del) free_slots_.push(node);
assert( this->sorted() );
}
@@ -238,6 +243,7 @@ void NodeContainer::clear() {
std::stack<Node*>().swap(free_slots_);
}
+
void NodeContainer::add_before(Node* add, Node* next_node){
add->set_next(next_node);
add->set_previous(next_node->previous());
@@ -271,7 +277,7 @@ bool NodeContainer::sorted() const {
return 0;
}
}
-
+
if ( !current->is_last() ) {
dout << "NodeContainer: Last Node not last" << std::endl;
return 0;
@@ -287,13 +293,17 @@ void swap(NodeContainer& first, NodeContainer& second) {
swap(first.last_node_, second.last_node_);
swap(first.size_, second.size_);
swap(first.unsorted_node_, second.unsorted_node_);
+ swap(first.node_counter_, second.node_counter_);
+ swap(first.lane_counter_, second.lane_counter_);
+ swap(first.node_lanes_, second.node_lanes_);
+ swap(first.free_slots_, second.free_slots_);
}
std::ostream& operator<< (std::ostream& stream, const NodeContainer& nc) {
for ( ConstNodeIterator it = nc.iterator(); it.good(); ++it ) {
stream << *it << "(" << (*it)->height() << ")";
- if (*it != nc.last()) stream << " <--> ";
+ if (*it != nc.last()) stream << " <--> ";
}
return stream;
}
diff --git a/src/node_container.h b/src/node_container.h
index 5a5d2d1..6d8146b 100644
--- a/src/node_container.h
+++ b/src/node_container.h
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -44,19 +44,19 @@ class ReverseConstNodeIterator;
class NodeContainer {
public:
NodeContainer();
- ~NodeContainer() {
- clear();
+ ~NodeContainer() {
+ clear();
for (std::vector<Node>* lane : node_lanes_) delete lane;
};
NodeContainer& operator=(NodeContainer nc) {
- swap(*this, nc);
+ swap(*this, nc);
return(*this);
};
NodeContainer(const NodeContainer &nc);
-
- NodeIterator iterator();
+
+ NodeIterator iterator();
NodeIterator iterator(Node* node);
ConstNodeIterator iterator() const;
ConstNodeIterator iterator(Node* node) const;
@@ -64,7 +64,7 @@ class NodeContainer {
ReverseConstNodeIterator reverse_iterator(Node* node) const;
// Create Nodes
- Node* createNode(double height, size_t label = 0) {
+ Node* createNode(double height, size_t label = 0) {
// Use the slot of a previously deleted node if possible
if (free_slots_.size() > 0) {
Node* node = free_slots_.top();
@@ -82,9 +82,33 @@ class NodeContainer {
new_lane->reserve(10000);
node_lanes_.push_back(new_lane);
}
- }
+ }
(*node_lanes_.at(lane_counter_))[node_counter_] = Node(height, label);
- return &*(node_lanes_[lane_counter_]->begin() + node_counter_++);
+ return &*(node_lanes_[lane_counter_]->begin() + node_counter_++);
+ }
+
+ // Create Nodes
+ Node* createNode(const Node copiedNode) {
+ // Use the slot of a previously deleted node if possible
+ if (free_slots_.size() > 0) {
+ Node* node = free_slots_.top();
+ free_slots_.pop();
+ *node = Node(copiedNode);
+ return node;
+ }
+
+ // Otherwise, use a new slot
+ if (node_counter_ >= 10000) {
+ ++lane_counter_;
+ node_counter_ = 0;
+ if (lane_counter_ == node_lanes_.size()) {
+ std::vector<Node>* new_lane = new std::vector<Node>();
+ new_lane->reserve(10000);
+ node_lanes_.push_back(new_lane);
+ }
+ }
+ (*node_lanes_.at(lane_counter_))[node_counter_] = Node(copiedNode);
+ return &*(node_lanes_[lane_counter_]->begin() + node_counter_++);
}
void push_back(Node* node);
@@ -94,14 +118,14 @@ class NodeContainer {
void move(Node *node, const double new_height);
void clear();
- Node* at(size_t nr) const;
- Node const* get(size_t nr) const { return at(nr); };
+ Node* at(size_t nr) const;
+ Node const* get(size_t nr) const { return at(nr); };
Node* first() const { return first_node_; };
Node* last() const { return last_node_; };
-
- size_t size() const { return size_; };
- bool sorted() const;
+
+ size_t size() const { return size_; };
+ bool sorted() const;
#ifdef UNITTEST
friend class TestNodeContainer;
@@ -119,7 +143,7 @@ class NodeContainer {
Node* first_node_;
Node* last_node_;
-
+
void set_first(Node* node) { first_node_ = node; }
void set_last(Node* node) { last_node_ = node; }
@@ -141,14 +165,14 @@ class NodeIterator {
NodeIterator(Node* node) { current_node_ = node; };
~NodeIterator() {};
- Node* operator*() {
- if (current_node_ == NULL)
+ Node* operator*() {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
- return current_node_;
+ return current_node_;
}
-
+
Node* operator++() {
- if (current_node_ == NULL)
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
if ( current_node_->is_last() ) current_node_ = NULL;
@@ -157,25 +181,25 @@ class NodeIterator {
}
Node* operator--() {
- if (current_node_ == NULL)
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
if ( current_node_->is_first() ) current_node_ = NULL;
else current_node_ = current_node_->previous();
return current_node_;
}
-
+
Node* operator++(int) {
- if (current_node_ == NULL)
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
-
+
Node* ret = current_node_;
if ( current_node_->is_last() ) current_node_ = NULL;
else current_node_ = current_node_->next();
return ret;
}
- bool good() const {
+ bool good() const {
return current_node_ != NULL;
}
@@ -183,7 +207,7 @@ class NodeIterator {
if ( good() ) return current_node_->height();
else return DBL_MAX;
}
-
+
#ifdef UNITTEST
friend class TestNodeContainer;
#endif
@@ -200,14 +224,14 @@ class ConstNodeIterator {
ConstNodeIterator(Node const* node) { current_node_ = node; };
~ConstNodeIterator() {};
- Node const* operator*() {
- if (current_node_ == NULL)
+ Node const* operator*() {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
- return current_node_;
+ return current_node_;
}
-
- Node const* operator++() {
- if (current_node_ == NULL)
+
+ Node const* operator++() {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
if ( current_node_->is_last() ) current_node_ = NULL;
@@ -215,27 +239,27 @@ class ConstNodeIterator {
return current_node_;
}
-
+
Node const* operator--() {
- if (current_node_ == NULL)
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
if ( current_node_->is_first() ) current_node_ = NULL;
else current_node_ = current_node_->previous();
return current_node_;
}
-
- Node const* operator++(int) {
- if (current_node_ == NULL)
+
+ Node const* operator++(int) {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
-
+
Node const* ret = current_node_;
if ( current_node_->is_last() ) current_node_ = NULL;
else current_node_ = current_node_->next();
return ret;
}
-
- bool good() const {
+
+ bool good() const {
return current_node_ != NULL;
}
@@ -256,14 +280,14 @@ class ReverseConstNodeIterator {
ReverseConstNodeIterator(Node const* node) { current_node_ = node; };
~ReverseConstNodeIterator() {};
- Node const* operator*() {
- if (current_node_ == NULL)
+ Node const* operator*() {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
- return current_node_;
+ return current_node_;
}
- Node const* operator++() {
- if (current_node_ == NULL)
+ Node const* operator++() {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
if ( current_node_->is_first() ) current_node_ = NULL;
@@ -271,23 +295,23 @@ class ReverseConstNodeIterator {
return current_node_;
}
- Node const* operator++(int) {
- if (current_node_ == NULL)
+ Node const* operator++(int) {
+ if (current_node_ == NULL)
throw std::out_of_range("Node iterator out of range");
-
+
Node const* ret = current_node_;
if ( current_node_->is_first() ) current_node_ = NULL;
else current_node_ = current_node_->previous();
return ret;
}
- bool good() const {
+ bool good() const {
return current_node_ != NULL;
}
private:
Node const* current_node_;
-};
+};
inline NodeIterator NodeContainer::iterator() { return NodeIterator(*this); }
diff --git a/src/param.cc b/src/param.cc
index 0004518..141f582 100644
--- a/src/param.cc
+++ b/src/param.cc
@@ -50,7 +50,8 @@ Model Param::parse() {
bool tmrca = false,
newick_trees = false,
orientedForest = false,
- sfs = false;
+ sfs = false,
+ transpose = false;
// The minimal time at which -eM, -eN, -eG, -eI, -ema and -es are allowed to happen. Is
@@ -394,7 +395,11 @@ Model Param::parse() {
else if (*argv_i == "-print-model" || *argv_i == "--print-model") {
this->set_print_model(true);
- }
+ }
+
+ else if (*argv_i == "-transpose-segsites" || *argv_i == "--transpose-segsites") {
+ transpose = true;
+ }
// ------------------------------------------------------------------
@@ -432,6 +437,7 @@ Model Param::parse() {
}
if (tmrca) model.addSummaryStatistic(std::make_shared<TMRCA>());
if (seg_sites.get() != NULL) model.addSummaryStatistic(seg_sites);
+ if (seg_sites.get() != NULL && transpose ) seg_sites->set_transpose(true);
if (sfs) {
if (seg_sites == NULL)
throw std::invalid_argument("You need to give a mutation rate ('-t') to simulate a SFS");
diff --git a/src/random/mersenne_twister.cc b/src/random/mersenne_twister.cc
index e934e2f..e75a011 100644
--- a/src/random/mersenne_twister.cc
+++ b/src/random/mersenne_twister.cc
@@ -40,13 +40,6 @@ MersenneTwister::MersenneTwister(const bool use_seed, size_t seed){
this->construct_common(seed);
}
-MersenneTwister::MersenneTwister(FastFunc* ff):RandomGenerator(ff) {
- this->construct_common(generateRandomSeed());
-}
-
-MersenneTwister::MersenneTwister(const size_t seed, FastFunc* ff):RandomGenerator(ff) {
- this->construct_common(seed);
-}
/**
* @brief Generates a random seed using entropy provided by the operating
diff --git a/src/random/mersenne_twister.h b/src/random/mersenne_twister.h
index 64de12d..00e0684 100644
--- a/src/random/mersenne_twister.h
+++ b/src/random/mersenne_twister.h
@@ -24,19 +24,26 @@
#define scrm_src_random_mersenne_twister
#include <random>
+#include <memory>
#include "random_generator.h"
class MersenneTwister : public RandomGenerator
{
public:
MersenneTwister();
- MersenneTwister(FastFunc* ff);
MersenneTwister(const size_t seed);
MersenneTwister(const bool use_seed, size_t seed);
- MersenneTwister(const size_t seed, FastFunc* ff);
- virtual ~MersenneTwister() {};
- void initialize() {};
+ MersenneTwister(std::shared_ptr<FastFunc> ff):RandomGenerator(ff) {
+ this->construct_common(generateRandomSeed());
+ }
+
+ MersenneTwister(const size_t seed, std::shared_ptr<FastFunc> ff):RandomGenerator(ff) {
+ this->construct_common(seed);
+ }
+
+ ~MersenneTwister() {};
+
void set_seed(const size_t seed);
void construct_common(const size_t seed);
diff --git a/src/random/random_generator.cc b/src/random/random_generator.cc
index 3ed3dc9..8f3dd34 100644
--- a/src/random/random_generator.cc
+++ b/src/random/random_generator.cc
@@ -74,5 +74,5 @@ double RandomGenerator::sampleExpoExpoLimit(double b, double c, double limit){
assert( result > 0 );
return result;
}
- }
+ }
}
diff --git a/src/random/random_generator.h b/src/random/random_generator.h
index 30fdb2c..96374d7 100644
--- a/src/random/random_generator.h
+++ b/src/random/random_generator.h
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
@@ -27,6 +27,7 @@
#include <cassert>
#include <cmath>
+#include <memory>
#include "fastfunc.h"
@@ -35,10 +36,10 @@ class RandomGenerator
{
friend class MersenneTwister;
public:
- RandomGenerator() { this->ff_ = new FastFunc(); };
- RandomGenerator( FastFunc* ff ) { this->ff_ = ff; };
+ RandomGenerator() : ff_(std::make_shared<FastFunc>()) { };
+ RandomGenerator(std::shared_ptr<FastFunc> ff) : ff_(ff) { };
+
virtual ~RandomGenerator() {}
- void clearFastFunc() { delete this->ff_; }
//Getters & Setters
size_t seed() const { return seed_; }
@@ -47,7 +48,6 @@ class RandomGenerator
this->seed_ = seed;
}
- virtual void initialize() =0;
virtual double sample() =0;
// Base class methods
@@ -57,13 +57,13 @@ class RandomGenerator
}
// Uniformly samples a number out of 0, ..., range-1
- // Unit tested
+ // Unit tested
int sampleInt(const int max_value) {
return(static_cast<int>(this->sample()*max_value));
}
- // Samples from an exponential distribution
- // Unit tested
+ // Samples from an exponential distribution
+ // Unit tested
double sampleExpo(const double lambda) {
return sampleUnitExponential() / lambda;
}
@@ -82,7 +82,8 @@ class RandomGenerator
#endif
// fast functions
- FastFunc *ff() { return this->ff_; }
+ std::shared_ptr<FastFunc> ff() { return this->ff_; }
+
protected:
// Sample from a unit exponential distribution
@@ -99,7 +100,7 @@ class RandomGenerator
size_t seed_;
// cache for a unit-exponentially distributed variable
double unit_exponential_;
- FastFunc *ff_;
+ std::shared_ptr<FastFunc> ff_;
};
#endif
diff --git a/src/scrm.cc b/src/scrm.cc
index 27fe937..6f85604 100644
--- a/src/scrm.cc
+++ b/src/scrm.cc
@@ -64,7 +64,6 @@ int main(int argc, char *argv[]){
// Loop over the independent loci/chromosomes
for (size_t rep_i=0; rep_i < model.loci_number(); ++rep_i) {
-
// Mark the start of a new independent sample
*output << std::endl << "//" << std::endl;
@@ -79,20 +78,21 @@ int main(int argc, char *argv[]){
forest.sampleNextGenealogy();
forest.printSegmentSumStats(*output);
}
-
+ assert(forest.next_base() == model.loci_length());
+
forest.printLocusSumStats(*output);
forest.clear();
}
- // Clean-up and exit
- rg.clearFastFunc();
return EXIT_SUCCESS;
}
- catch (const std::exception &e)
- {
+
+ catch (const std::exception &e) {
std::cerr << "Error: " << e.what() << std::endl;
std::cerr << "Try 'scrm --help' for more information." << std::endl;
return EXIT_FAILURE;
}
}
+
#endif
+
diff --git a/src/summary_statistics/newick_tree.cc b/src/summary_statistics/newick_tree.cc
index 2ab3ce1..015d98e 100644
--- a/src/summary_statistics/newick_tree.cc
+++ b/src/summary_statistics/newick_tree.cc
@@ -30,7 +30,11 @@ void NewickTree::calculate(const Forest &forest) {
void NewickTree::printSegmentOutput(std::ostream &output) const {
if (segment_length_ == 0.0) return;
- if (has_rec_) output << "[" << segment_length_ << "]";
+ if (has_rec_) {
+ double intpart; // dummy variable for modf
+ if (modf(segment_length_, &intpart) == 0.0) output << "[" << (size_t)segment_length_ << "]";
+ else output << "[" << segment_length_ << "]";
+ }
output << tree_ << ";" << std::endl;
}
diff --git a/src/summary_statistics/newick_tree.h b/src/summary_statistics/newick_tree.h
index e37ce41..4d58248 100644
--- a/src/summary_statistics/newick_tree.h
+++ b/src/summary_statistics/newick_tree.h
@@ -56,7 +56,6 @@ class NewickTree : public SummaryStatistic
//Virtual methods
void calculate(const Forest &forest);
void printSegmentOutput(std::ostream &output) const;
- void printLocusOutput(std::ostream &output) const {};
NewickTree* clone() const { return new NewickTree(precision_, has_rec_); };
diff --git a/src/summary_statistics/oriented_forest.h b/src/summary_statistics/oriented_forest.h
index 17b08e9..67f17d8 100644
--- a/src/summary_statistics/oriented_forest.h
+++ b/src/summary_statistics/oriented_forest.h
@@ -40,7 +40,6 @@ class OrientedForest : public SummaryStatistic
//Virtual methods
void calculate(const Forest &forest);
- void printLocusOutput(std::ostream &output) const {};
void printSegmentOutput(std::ostream &output) const;
void clear() { }
diff --git a/src/summary_statistics/seg_sites.cc b/src/summary_statistics/seg_sites.cc
index 970d5eb..1025b40 100644
--- a/src/summary_statistics/seg_sites.cc
+++ b/src/summary_statistics/seg_sites.cc
@@ -32,6 +32,7 @@ void SegSites::calculate(const Forest &forest) {
while (position_at < forest.next_base()) {
TreePoint mutation = forest.samplePoint();
+ heights_.push_back(mutation.height() / (4 * forest.model().default_pop_size));
haplotypes_.push_back(getHaplotypes(mutation, forest));
if (forest.model().getSequenceScaling() == absolute) {
positions_.push_back(position_at);
@@ -39,23 +40,40 @@ void SegSites::calculate(const Forest &forest) {
positions_.push_back(position_at / forest.model().loci_length());
}
position_at += forest.random_generator()->sampleExpo(forest.getLocalTreeLength() * forest.model().mutation_rate());
- }
+ }
set_position(forest.next_base());
}
void SegSites::printLocusOutput(std::ostream &output) const {
- output << "segsites: "<< countMutations() << std::endl;
- if ( countMutations() == 0 ) return;
-
- output << "positions: " << positions_ << std::endl;
-
- for (size_t i = 0; i < haplotypes_.at(0).size(); i++){
- for (size_t j = 0; j < haplotypes_.size(); j++){
- output << haplotypes_[j][i];
+ if ( transpose_ ) {
+ output << "transposed segsites: " << countMutations() << std::endl;
+ if ( countMutations() == 0 ) return;
+ output << "position time";
+ for (size_t i = 0; i < haplotypes_.at(0).size(); i++){
+ output << " " << i+1;
}
output <<"\n";
+
+ for (size_t j = 0; j < haplotypes_.size(); j++){
+ output << positions_[j] << " " << heights_[j];
+ for (size_t i = 0; i < haplotypes_.at(0).size(); i++){
+ output << " " << haplotypes_[j][i];
+ }
+ output <<"\n";
+ }
+ } else {
+ output << "segsites: " << countMutations() << std::endl;
+ if ( countMutations() == 0 ) return;
+ output << "positions: " << positions_ << std::endl;
+
+ for (size_t i = 0; i < haplotypes_.at(0).size(); i++){
+ for (size_t j = 0; j < haplotypes_.size(); j++){
+ output << haplotypes_[j][i];
+ }
+ output <<"\n";
+ }
}
}
diff --git a/src/summary_statistics/seg_sites.h b/src/summary_statistics/seg_sites.h
index 6b1b35e..0fd0bc9 100644
--- a/src/summary_statistics/seg_sites.h
+++ b/src/summary_statistics/seg_sites.h
@@ -34,7 +34,7 @@
class SegSites : public SummaryStatistic
{
public:
- SegSites() { set_position(0.0); }
+ SegSites( ) { set_position(0.0); set_transpose(false);}
~SegSites() {}
#ifdef UNITTEST
@@ -44,6 +44,7 @@ class SegSites : public SummaryStatistic
//Virtual methods
void calculate(const Forest &forest);
void printLocusOutput(std::ostream &output) const;
+
SegSites* clone() const { return new SegSites(*this); }
void clear() {
@@ -60,16 +61,20 @@ class SegSites : public SummaryStatistic
std::valarray<bool> const* getHaplotype(const size_t mutation) const {
return &(haplotypes_.at(mutation));
}
+ void set_transpose(const bool transpose) { transpose_ = transpose; };
+ bool get_transpose() const { return transpose_; }
private:
std::valarray<bool> getHaplotypes(TreePoint mutation, const Forest &Forest);
std::vector<double> positions_;
+ std::vector<double> heights_;
std::vector<std::valarray<bool>> haplotypes_;
void traversal(Node const* node, std::valarray<bool> &haplotype) const;
void set_position(const double position) { position_ = position; };
double position_;
+ bool transpose_;
};
#endif
diff --git a/src/summary_statistics/summary_statistic.h b/src/summary_statistics/summary_statistic.h
index e557a95..07f8c63 100644
--- a/src/summary_statistics/summary_statistic.h
+++ b/src/summary_statistics/summary_statistic.h
@@ -35,12 +35,12 @@ class SummaryStatistic
// Virtual methods
virtual void calculate(const Forest &forest) =0;
- virtual void printLocusOutput(std::ostream &output) const =0;
virtual void clear() =0;
virtual SummaryStatistic* clone() const =0;
// Optional methods
- virtual void printSegmentOutput(std::ostream &output) const { };
+ virtual void printLocusOutput(std::ostream &output) const { (void) output; };
+ virtual void printSegmentOutput(std::ostream &output) const { (void) output; };
};
#endif
diff --git a/src/summary_statistics/tmrca.cc b/src/summary_statistics/tmrca.cc
index 1b16c4a..cdddc57 100644
--- a/src/summary_statistics/tmrca.cc
+++ b/src/summary_statistics/tmrca.cc
@@ -1,10 +1,10 @@
/*
* scrm is an implementation of the Sequential-Coalescent-with-Recombination Model.
- *
+ *
* Copyright (C) 2013, 2014 Paul R. Staab, Sha (Joe) Zhu, Dirk Metzler and Gerton Lunter
- *
+ *
* This file is part of scrm.
- *
+ *
* scrm is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
@@ -14,7 +14,7 @@
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
- *
+ *
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
@@ -27,8 +27,9 @@ void TMRCA::calculate(const Forest &forest) {
tree_length_.push_back(forest.getLocalTreeLength(true));
}
+
void TMRCA::printLocusOutput(std::ostream &output) const {
for (size_t i = 0; i < tmrca_.size(); ++i) {
- output << "time:\t" << tmrca_.at(i) << " \t" << tree_length_.at(i) << "\n";
+ output << "time:\t" << tmrca_.at(i) << " \t" << tree_length_.at(i) << "\n";
}
}
diff --git a/src/time_interval.h b/src/time_interval.h
index c53cad0..3cb866c 100644
--- a/src/time_interval.h
+++ b/src/time_interval.h
@@ -136,14 +136,7 @@ class TimeIntervalIterator {
* @return Nothing. Nodes are saved in contemporaries_.
*/
inline void TimeIntervalIterator::searchContemporaries(Node *node) {
- // Prefer top-down search if the target is above .8 coalescence units in
- // sample space. This is relatively high, but the iterative bottom-up approach
- // is faster than the recursion.
- if (node->height() >= contemporaries()->buffer_time()) {
- searchContemporariesBottomUp(node, true);
- } else {
- searchContemporariesBottomUp(node);
- }
+ searchContemporariesBottomUp(node, node->height() >= contemporaries()->buffer_time());
}
#endif
diff --git a/tests/algorithmtest/test_algorithm.cc b/tests/algorithmtest/test_algorithm.cc
index 3d51734..1b1a624 100644
--- a/tests/algorithmtest/test_algorithm.cc
+++ b/tests/algorithmtest/test_algorithm.cc
@@ -94,7 +94,6 @@ class TestAlgorithm : public CppUnit::TestCase {
}
void tearDown() {
- rg->clearFastFunc();
delete rg;
}
@@ -179,9 +178,7 @@ class TestAlgorithm : public CppUnit::TestCase {
void testSizeChange() {
Model model = Model(0);
- char *argv[] = { "scrm", "10", "1", "-I", "3", "3", "3", "4", "0.5", "-eN", "0.1", "0.05", "-eN", "0.2", "0.5" };
- model = Param(15, argv).parse();
- model.setRecombinationRate(1, false, true);
+ model = Param("10 1 -r 1 100 -I 3 3 3 4 0.5 -eN 0.1 0.05 -eN 0.2 0.5").parse();
testTree(model, 1000, 3.75, 2.68, 9.31, 5.67);
model.set_window_length_seq(5);
@@ -190,27 +187,19 @@ class TestAlgorithm : public CppUnit::TestCase {
void testGrowth() {
- Model model = Model(10);
- model.setRecombinationRate(1, false, true);
- model.addGrowthRates(0.0, 5, true, true);
- model.finalize();
+ Model model = Param("10 30 -r 50 50 -G 5").parse();
testTree(model, 2500, 0.321, 0.089, 1.31, 0.28);
- char *argv0[] = { "scrm", "10", "30", "-G", "-0.5", "-eG", "0.75", "2" };
- model = Param(8, argv0).parse();
- model.setRecombinationRate(1, false, true);
- model.set_window_length_seq(5);
+ model = Param("10 30 -r 50 50 -G -0.5 -eG 0.75 2 -l 5").parse();
testTree(model, 2500, 0.918, 0.38, 2.95, 1.00);
- char *argv[] = { "scrm", "4", "30", "-G", "-2.5", "-eN", "1", "0.25",
- "-eG", "2", "0.0", "-eN", "2.5", "0.25" };
- model = Param(14, argv).parse();
- model.setRecombinationRate(1, false, true);
+ model = Param("4 30 -r 50 50 -G -2.5 -eN 1 0.25 -eG 2 0.0 -eN 2.5 0.25 -l 5").parse();
testTree(model, 1000, 0.964, 0.34, 2.48, 0.97);
}
void testSplit() {
Model model = Model(0);
+ model.setLocusLength(100);
model.setRecombinationRate(1, false, true);
model.set_population_number(2);
std::vector<size_t> sample_size;
@@ -223,9 +212,7 @@ class TestAlgorithm : public CppUnit::TestCase {
model.finalize();
testTree(model, 2500, 1.51, 0.55, 5.20, 1.44);
- char *argv[] = { "scrm", "15", "5", "-I", "3", "7", "3", "5", "0.5", "-ej", "0.5", "3", "2", "-ej", "1.0", "2", "1" };
- model = Param(17, argv).parse();
- model.setRecombinationRate(1, false, true);
+ model = Param("15 5 -r 50 50 -I 3 7 3 5 0.5 -ej 0.5 3 2 -ej 1.0 2 1").parse();
testTree(model, 2500, 1.60, 0.54, 7.29, 1.48);
model.set_window_length_seq(5);
@@ -233,18 +220,14 @@ class TestAlgorithm : public CppUnit::TestCase {
}
void testMerge() {
- Model model;
- char *argv1[] = { "scrm", "20", "10", "-I", "2", "10", "10", "1.5", "-es", "1.6", "2", "0.5", "-eM", "2.0", "1" };
- model = Param(15, argv1).parse();
- model.setRecombinationRate(1, false, true);
+ Model model = Param("20 10 -r 50 50 -I 2 10 10 1.5 -es 1.6 2 0.5 -eM 2.0 1").parse();
testTree(model, 1000, 2.88, 2.26, 9.36, 4.87);
model.set_window_length_seq(5);
testTree(model, 1000, 2.88, 2.26, 9.36, 4.87);
- char *argv2[] = { "scrm", "20", "10", "-I", "2", "10", "10", "1.5", "-es", "0.9", "1", "0.8", "-es", "1.6", "2", "0.5", "-eM", "2.0", "1" };
- model = Param(19, argv2).parse();
- model.setRecombinationRate(1, false, true);
+
+ model = Param("20 10 -r 50 50 -I 2 10 10 1.5 -es 0.9 1 0.8 -es 1.6 2 0.5 -eM 2.0 1").parse();
testTree(model, 1000, 3.82, 3.32, 11.39, 7.09);
model.set_window_length_seq(5);
diff --git a/tests/test_binaries.sh b/tests/test_binaries.sh
index e44c40c..30c5960 100755
--- a/tests/test_binaries.sh
+++ b/tests/test_binaries.sh
@@ -35,31 +35,31 @@ function test_scrm {
echo "Testing Initial Tree"
- test_scrm 5 1000 -t 5 || exit 1
- test_scrm 3 1000 -t 5 -L -oSFS || exit 1
- test_scrm 100 20 -O || exit 1
- test_scrm 250 1 -T || exit 1
+ test_scrm 5 100 -t 5 || exit 1
+ test_scrm 3 100 -t 5 -L -oSFS -transpose-segsites || exit 1
+ test_scrm 50 20 -O || exit 1
+ test_scrm 25 1 -T || exit 1
echo ""
echo "Testing Recombinations"
- test_scrm 4 500 -r 5 100 || exit 1
- test_scrm 6 500 -r 1 100 -t 5 -L -T || exit 1
- test_scrm 8 500 -r 1 100 -t 5 -oSFS -O || exit 1
+ test_scrm 4 10 -r 5 100 || exit 1
+ test_scrm 6 10 -r 1 100 -t 5 -L -T -transpose-segsites || exit 1
+ test_scrm 8 10 -r 1 100 -t 5 -oSFS -O || exit 1
echo ""
echo "Testing Pruning"
- test_scrm 10 200 -r 10 500 -l 10 -t 5 -oSFS || exit 1
- test_scrm 3 500 -r 10 500 -l 0 -t 1 -oSFS || exit 1
- test_scrm 10 200 -r 10 500 -l 2r -L || exit 1
- test_scrm 10 200 -r 1 500 -l -1 -T || exit 1
+ test_scrm 10 20 -r 10 500 -l 10 -t 5 -oSFS || exit 1
+ test_scrm 3 50 -r 10 500 -l 0 -t 1 -oSFS || exit 1
+ test_scrm 10 20 -r 10 500 -l 2r -L || exit 1
+ test_scrm 10 20 -r 1 500 -l -1 -T || exit 1
echo ""
echo "Testing Migration"
- test_scrm 5 50 -r 5 100 -I 2 3 2 1.2 || exit 1
- test_scrm 10 1 -r 20 200 -I 5 2 2 2 2 2 0.75 -l 5 || exit 1
- test_scrm 10 1 -r 10 100 -I 2 7 3 0.5 -eM 0.3 1.1 --print-model || exit 1
- test_scrm 10 1 -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 100 -I 3 2 2 2 1.0 -eI 1.0 2 2 2 -eI 2.0 2 3 3 || exit 1
+ 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 -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"
@@ -68,25 +68,25 @@ echo "Testing Size Change"
echo ""
echo "Testing Splits & Merges"
- test_scrm 5 30 -r 20 200 -I 2 3 2 0.4 -ej 1.1 2 1 -l 25 || exit 1
- test_scrm 6 30 -r 20 200 -I 3 2 2 2 -ej 0.2 2 1 -ej 0.25 3 1 -l 25 || exit 1
+ test_scrm 5 3 -r 20 200 -I 2 3 2 0.4 -ej 1.1 2 1 -l 25 || exit 1
+ test_scrm 6 3 -r 20 200 -I 3 2 2 2 -ej 0.2 2 1 -ej 0.25 3 1 -l 25 || exit 1
test_scrm 20 2 -r 5 200 -I 2 10 10 1.5 -es 1.6 2 0.5 -eM 2.0 1 -l 25 || exit 1
test_scrm 20 2 -r 5 200 -I 2 10 10 1.5 -es 0.9 1 0.8 -es 1.6 2 0.5 -eM 2.0 1 -l 50 || exit 1
test_scrm 20 2 -r 5 200 -I 2 10 10 -eps 1.0 2 1 0.0 || exit 1
echo ""
echo "Testing Growth"
- test_scrm 5 100 -r 20 200 -G 1.5 -l 25 || exit 1
- test_scrm 5 60 -r 2 200 -G -2.5 -eG 1 0.0 -l 25 || exit 1
- test_scrm 4 30 -r 2 200 -G -2.5 -eN 1 0.25 -eG 2 0.0 -eN 2.5 0.25 -l 25 || exit 1
+ test_scrm 5 10 -r 20 200 -G 1.5 -l 25 || exit 1
+ test_scrm 5 6 -r 2 200 -G -2.5 -eG 1 0.0 -l 25 || exit 1
+ test_scrm 4 3 -r 2 200 -G -2.5 -eN 1 0.25 -eG 2 0.0 -eN 2.5 0.25 -l 25 || exit 1
echo ""
echo "Testing Variable Rates"
- test_scrm 3 200 -r 2 100 -t 5 -st 10 10 -sr 20 5 -st 30 1 -sr 40 0 -st 50 20 -T || exit 1
+ test_scrm 3 2 -r 2 100 -t 5 -st 10 10 -sr 20 5 -st 30 1 -sr 40 0 -st 50 20 -T || exit 1
echo ""
echo "Various Edge Cases"
test_scrm 6 1 -I 2 3 3 0.5 -r 1 100 -es 0 2 0.5 -ej 1 3 1 -t 1 || exit 1
- test_scrm 10 10 -es 1.0 1 0.5 -ej 1.0 2 1 || exit 1
+ test_scrm 10 1 -es 1.0 1 0.5 -ej 1.0 2 1 || exit 1
echo ""
diff --git a/tests/unittests/test_contemporaries_container.cc b/tests/unittests/test_contemporaries_container.cc
index d41ec0c..49ae2f5 100644
--- a/tests/unittests/test_contemporaries_container.cc
+++ b/tests/unittests/test_contemporaries_container.cc
@@ -38,10 +38,11 @@ class TestContemporariesContainer : public CppUnit::TestCase {
}
void tearDown() {
- delete rg, nc;
+ delete rg;
+ delete nc;
}
- void add() {
+ void add() {
// Vector
ContemporariesContainer cc = ContemporariesContainer(3, 10, rg);
CPPUNIT_ASSERT_EQUAL( (size_t)0, cc.size(0) );
@@ -168,7 +169,7 @@ class TestContemporariesContainer : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( (size_t)0, cc.size(0) );
CPPUNIT_ASSERT_EQUAL( (size_t)0, cc.size(1) );
CPPUNIT_ASSERT_EQUAL( (size_t)0, cc.size(2) );
-
+
cc = ContemporariesContainer(3, 1000, rg);
cc.add(node1);
cc.add(node2);
@@ -246,16 +247,16 @@ class TestContemporariesContainer : public CppUnit::TestCase {
for (size_t i = 0; i < 1000; ++i) {
CPPUNIT_ASSERT_EQUAL(node1, cc.sample(0));
CPPUNIT_ASSERT_EQUAL(node2, cc.sample(1));
- }
+ }
double count = 0;
for (size_t i = 0; i < 10000; ++i) {
Node* node = cc.sample(2);
CPPUNIT_ASSERT( node == node3 || node == node4 );
if (node == node3) ++count;
- }
+ }
count /= 10000;
- CPPUNIT_ASSERT( 0.49 < count && count < 0.51 );
+ CPPUNIT_ASSERT( 0.49 < count && count < 0.51 );
// Set
cc = ContemporariesContainer(3, 1000, rg);
@@ -267,16 +268,16 @@ class TestContemporariesContainer : public CppUnit::TestCase {
for (size_t i = 0; i < 1000; ++i) {
CPPUNIT_ASSERT_EQUAL(node1, cc.sample(0));
CPPUNIT_ASSERT_EQUAL(node2, cc.sample(1));
- }
+ }
count = 0;
for (size_t i = 0; i < 10000; ++i) {
Node* node = cc.sample(2);
CPPUNIT_ASSERT( node == node3 || node == node4 );
if (node == node3) ++count;
- }
+ }
count /= 10000;
- CPPUNIT_ASSERT( 0.49 < count && count < 0.51 );
+ CPPUNIT_ASSERT( 0.49 < count && count < 0.51 );
}
void buffer() {
@@ -368,7 +369,7 @@ class TestContemporariesContainer : public CppUnit::TestCase {
}
void empty() {
- // Vector
+ // Vector
ContemporariesContainer cc = ContemporariesContainer(3, 10, rg);
CPPUNIT_ASSERT( cc.empty() );
cc.add(node1);
diff --git a/tests/unittests/test_forest.cc b/tests/unittests/test_forest.cc
index dad9c2c..65e9f24 100644
--- a/tests/unittests/test_forest.cc
+++ b/tests/unittests/test_forest.cc
@@ -24,18 +24,18 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_TEST( testSelectFirstTime );
CPPUNIT_TEST( testSampleEventType );
CPPUNIT_TEST( testSampleEvent );
- CPPUNIT_TEST( testGetNodeState );
+ CPPUNIT_TEST( testGetNodeState );
CPPUNIT_TEST( testCut );
CPPUNIT_TEST( testImplementCoalescence );
- CPPUNIT_TEST( testBuildInitialTree );
- CPPUNIT_TEST( testImplementRecombination );
- CPPUNIT_TEST( testImplementFixTimeEvent );
+ CPPUNIT_TEST( testBuildInitialTree );
+ CPPUNIT_TEST( testImplementRecombination );
+ CPPUNIT_TEST( testImplementFixTimeEvent );
CPPUNIT_TEST( testPrintTree );
CPPUNIT_TEST( testCopyConstructor );
CPPUNIT_TEST( testCheckForNodeAtHeight );
CPPUNIT_TEST( testPrintLocusSumStats );
- CPPUNIT_TEST( testSampleNextPosition );
- CPPUNIT_TEST( testClear );
+ CPPUNIT_TEST( testSampleNextPosition );
+ CPPUNIT_TEST( testClear );
CPPUNIT_TEST_SUITE_END();
@@ -57,8 +57,10 @@ class TestForest : public CppUnit::TestCase {
}
void tearDown() {
- delete forest, forest_2pop;
- delete model, model_2pop;
+ delete forest;
+ delete model;
+ delete forest_2pop;
+ delete model_2pop;
delete rg;
}
@@ -84,7 +86,7 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT( forest->primary_root() == forest->nodes()->get(8) );
CPPUNIT_ASSERT( forest->getLocalTreeLength() == 24 );
CPPUNIT_ASSERT( forest->checkTree() );
-
+
forest->createExampleTree();
CPPUNIT_ASSERT_EQUAL((size_t)9, forest->nodes()->size());
CPPUNIT_ASSERT( forest->local_root() == forest->nodes()->get(8) );
@@ -100,7 +102,7 @@ class TestForest : public CppUnit::TestCase {
void testCalcRate() {
TimeIntervalIterator tii(forest, forest->nodes()->at(0));
- size_t pop_size = 2*forest->model().population_size(0);
+ double pop_size = 2*forest->model().population_size(0);
Node *node1 = new Node(0.1);
Node *node2 = new Node(0.2);
@@ -109,23 +111,23 @@ class TestForest : public CppUnit::TestCase {
forest->states_[0] = 0;
forest->states_[1] = 0;
CPPUNIT_ASSERT_NO_THROW( forest->calcRates(*tii) );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[0] );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[0] );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
forest->states_[0] = 1;
CPPUNIT_ASSERT_NO_THROW( forest->calcRates(*tii) );
- CPPUNIT_ASSERT_EQUAL( 4.0/pop_size, forest->rates_[0] );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
+ CPPUNIT_ASSERT( areSame(4.0/pop_size, forest->rates_[0]) );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
forest->states_[1] = 1;
CPPUNIT_ASSERT_NO_THROW( forest->calcRates(*tii) );
- CPPUNIT_ASSERT( areSame(9.0/pop_size, forest->rates_[0]) );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
- CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
+ CPPUNIT_ASSERT( areSame(9.0/pop_size, forest->rates_[0]) );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
+ CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
- // Coalescence with structure
+ // Coalescence with structure
forest_2pop->set_active_node(0, node1);
forest_2pop->set_active_node(1, node2);
forest_2pop->states_[0] = 1;
@@ -149,33 +151,33 @@ class TestForest : public CppUnit::TestCase {
TimeIntervalIterator tii3(forest_2pop, forest_2pop->nodes()->at(0));
CPPUNIT_ASSERT_NO_THROW( forest_2pop->calcRates(*tii3) );
- CPPUNIT_ASSERT( areSame(3.0/pop_size, forest_2pop->rates_[0]) );
- CPPUNIT_ASSERT( areSame(1.0/pop_size, forest_2pop->rates_[1]) );
- CPPUNIT_ASSERT( areSame(0.0/pop_size, forest_2pop->rates_[2]) );
- CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[0] );
- CPPUNIT_ASSERT_EQUAL( (size_t)0, forest_2pop->active_nodes_timelines_[1] );
+ CPPUNIT_ASSERT( areSame(3.0/pop_size, forest_2pop->rates_[0]) );
+ CPPUNIT_ASSERT( areSame(1.0/pop_size, forest_2pop->rates_[1]) );
+ CPPUNIT_ASSERT( areSame(0.0/pop_size, forest_2pop->rates_[2]) );
+ CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[0] );
+ CPPUNIT_ASSERT_EQUAL( (size_t)0, forest_2pop->active_nodes_timelines_[1] );
forest_2pop->writable_model()->increaseTime();
CPPUNIT_ASSERT_NO_THROW( forest_2pop->calcRates(*tii3) );
- //CPPUNIT_ASSERT_EQUAL( 0.0/pop_size, forest_2pop->rates_[0] );
- //CPPUNIT_ASSERT_EQUAL( 1.0/pop_size, forest_2pop->rates_[1] );
- //CPPUNIT_ASSERT_EQUAL( 3.0/pop_size, forest_2pop->rates_[2] );
- CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[0] );
- CPPUNIT_ASSERT_EQUAL( (size_t)2, forest_2pop->active_nodes_timelines_[1] );
+ //CPPUNIT_ASSERT_EQUAL( 0.0/pop_size, forest_2pop->rates_[0] );
+ //CPPUNIT_ASSERT_EQUAL( 1.0/pop_size, forest_2pop->rates_[1] );
+ //CPPUNIT_ASSERT_EQUAL( 3.0/pop_size, forest_2pop->rates_[2] );
+ CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[0] );
+ CPPUNIT_ASSERT_EQUAL( (size_t)2, forest_2pop->active_nodes_timelines_[1] );
node2->set_population(1);
CPPUNIT_ASSERT_NO_THROW( forest_2pop->calcRates(*tii3) );
- //CPPUNIT_ASSERT_EQUAL( 0.0/pop_size, forest_2pop->rates_[0] );
- //CPPUNIT_ASSERT( areSame(3.0/pop_size, forest_2pop->rates_[1]) );
- //CPPUNIT_ASSERT_EQUAL( 0.0/pop_size, forest_2pop->rates_[2] );
- CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[0] );
- CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[1] );
+ //CPPUNIT_ASSERT_EQUAL( 0.0/pop_size, forest_2pop->rates_[0] );
+ //CPPUNIT_ASSERT( areSame(3.0/pop_size, forest_2pop->rates_[1]) );
+ //CPPUNIT_ASSERT_EQUAL( 0.0/pop_size, forest_2pop->rates_[2] );
+ CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[0] );
+ CPPUNIT_ASSERT_EQUAL( (size_t)1, forest_2pop->active_nodes_timelines_[1] );
delete node1;
delete node2;
}
- void testGetNodeState() {
+ void testGetNodeState() {
CPPUNIT_ASSERT( forest->getNodeState(forest->getNodes()->get(8), 11) == 1 );
CPPUNIT_ASSERT( forest->getNodeState(forest->getNodes()->get(6), 11) == 2 );
CPPUNIT_ASSERT( forest->getNodeState(forest->getNodes()->get(7), 11) == 1 );
@@ -211,7 +213,7 @@ class TestForest : public CppUnit::TestCase {
void testSampleEventType() {
Event event;
- Forest *forest2 = forest_2pop;
+ Forest *forest2 = forest_2pop;
forest2->nodes()->at(0)->set_population(1);
forest2->nodes()->at(1)->set_population(1);
@@ -247,9 +249,9 @@ class TestForest : public CppUnit::TestCase {
};
// Coalescence of two nodes
- // active_node 0: Pop 1, 2 Contemporaries
+ // active_node 0: Pop 1, 2 Contemporaries
// active_node 1: Pop 0, 2 Contemporaries
- // => 50% each
+ // => 50% each
forest2->states_[1] = 1;
forest2->calcRates(*tii);
size_t count = 0, count2 = 0;
@@ -261,7 +263,7 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT( 4900 < count && count < 5100 ); // ~5000
// Test with Pw Coalescence
- // active_node 0: Pop 1, 2 Contemporaries
+ // active_node 0: Pop 1, 2 Contemporaries
// active_node 1: Pop 1, 2 Contemporaries
// 4/5 Prob of Coalescence, 1/5 of Pw coalescence
forest2->writable_model()->resetTime();
@@ -273,14 +275,14 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT( event.isCoalescence() || event.isPwCoalescence() );
count += event.isCoalescence();
};
- CPPUNIT_ASSERT( 7900 < count && count < 8100 );
+ CPPUNIT_ASSERT( 7900 < count && count < 8100 );
// Test with migration
- // active_node 0: Pop 1, 2 Contemporaries
+ // active_node 0: Pop 1, 2 Contemporaries
// active_node 1: Pop 1, 2 Contemporaries
// => Coal: 4/2Ne, Pw Coal: 1/2Ne
// Migration: 2 * 5/4Ne
- // => 40% Coal, 10% Pw Coal, 50% Mig
+ // => 40% Coal, 10% Pw Coal, 50% Mig
forest2->writable_model()->increaseTime();
forest2->writable_model()->increaseTime();
forest2->calcRates(*tii);
@@ -296,10 +298,10 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT( 3900 < count2 && count2 < 4100 );
- // Coalescence and Recombination
- // active_node 0: Pop 1, 2 Contemporaries => Coal rate: 2 / 2 * Ne = 1/Ne
- // active_node 1: Pop 1, Recombination => Rec rate: 10 Bases * 0.4 / 4Ne = 1/Ne
- // => 50% Recombination, 50% Coalescence
+ // Coalescence and Recombination
+ // active_node 0: Pop 1, 2 Contemporaries => Coal rate: 2 / 2 * Ne = 1/Ne
+ // active_node 1: Pop 1, Recombination => Rec rate: 10 Bases * 0.4 / 4Ne = 1/Ne
+ // => 50% Recombination, 50% Coalescence
forest2->writable_model()->setLocusLength(101);
forest2->writable_model()->setRecombinationRate(0.4, false, true);
forest2->set_current_base(10);
@@ -345,7 +347,7 @@ class TestForest : public CppUnit::TestCase {
forest2->set_current_base(10.0);
forest2->set_next_base(15.0);
forest2->set_next_base(20.0);
-
+
forest2->current_rec_ += 2;
forest2->active_node(0)->make_nonlocal(forest2->current_rec_ - 1);
forest2->active_node(1)->make_nonlocal(forest2->current_rec_ - 2);
@@ -373,7 +375,7 @@ class TestForest : public CppUnit::TestCase {
};
// Combinations With Growth
- forest2->writable_model()->resetTime();
+ forest2->writable_model()->resetTime();
forest2->writable_model()->increaseTime();
forest2->writable_model()->increaseTime();
forest2->writable_model()->increaseTime();
@@ -403,9 +405,6 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT_NO_THROW( forest2->sampleEventType(0.5, 1, *tii, event) );
CPPUNIT_ASSERT( event.isCoalescence() || event.isPwCoalescence() );
};
-
- delete forest2->writable_model();
- delete forest2;
}
void testCalcRecombinationRate() {
@@ -454,7 +453,7 @@ class TestForest : public CppUnit::TestCase {
Event event;
double tmp_event_time = 0.0;
for (size_t i = 0; i < 1000; ++i) {
- forest2.sampleEvent(*tii, tmp_event_time, event);
+ forest2.sampleEvent(*tii, tmp_event_time, event);
CPPUNIT_ASSERT( event.isNoEvent() || ( 0 <= event.time() && event.time() < forest2.nodes()->at(4)->height() ) );
CPPUNIT_ASSERT( event.isNoEvent() || event.isCoalescence() );
}
@@ -462,7 +461,7 @@ class TestForest : public CppUnit::TestCase {
++tii;
forest2.calcRates(*tii);
for (size_t i = 0; i < 1000; ++i) {
- forest2.sampleEvent(*tii, tmp_event_time, event);
+ forest2.sampleEvent(*tii, tmp_event_time, event);
CPPUNIT_ASSERT( event.isNoEvent() || ( forest2.nodes()->at(4)->height() <= event.time() && event.time() < forest2.nodes()->at(5)->height() ) );
CPPUNIT_ASSERT( event.isNoEvent() || event.isCoalescence() );
}
@@ -544,10 +543,10 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT( forest->checkTree() == 1 );
// In-Between Nodes should be pruned, iff they are of same age
- Node *parent = new Node(20),
- *inbetween1 = new Node(19),
- *inbetween2 = new Node(18),
- *child = new Node(17);
+ Node *parent = forest->nodes()->createNode(20),
+ *inbetween1 = forest->nodes()->createNode(19),
+ *inbetween2 = forest->nodes()->createNode(18),
+ *child = forest->nodes()->createNode(17);
forest->set_current_base(13);
inbetween2->make_nonlocal(forest->current_rec_);
@@ -595,7 +594,7 @@ class TestForest : public CppUnit::TestCase {
model = Model(2);
model.set_population_number(2);
- model.addSampleSizes(0.0, std::vector<size_t>(2, 1));
+ model.addSampleSizes(0.0, std::vector<size_t>(2, 1));
model.addSymmetricMigration(0.0, 5.0, true, true);
model.finalize();
Forest frst = Forest(&model, rg);
@@ -603,7 +602,7 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( true, frst.checkTree() );
size_t i = 0;
for (size_t j = 0; j < 4; ++j) {
- CPPUNIT_ASSERT( frst.nodes()->at(j)->population() == 0 ||
+ CPPUNIT_ASSERT( frst.nodes()->at(j)->population() == 0 ||
frst.nodes()->at(j)->population() == 1 );
i += frst.nodes()->at(j)->population();
}
@@ -613,10 +612,10 @@ class TestForest : public CppUnit::TestCase {
model.set_population_number(2);
std::vector<size_t> sample_size(2, 1);
sample_size.at(1) = 0;
- model.addSampleSizes(0.0, sample_size);
+ model.addSampleSizes(0.0, sample_size);
sample_size.at(0) = 0;
sample_size.at(1) = 1;
- model.addSampleSizes(1.0, sample_size);
+ model.addSampleSizes(1.0, sample_size);
model.addSymmetricMigration(0.0, 5.0, true, true);
model.finalize();
frst = Forest(&model, rg);
@@ -669,7 +668,7 @@ class TestForest : public CppUnit::TestCase {
event.setToRecombination(forest->active_node(0), 0);
CPPUNIT_ASSERT_NO_THROW( forest->implementRecombination(event, tii) );
- CPPUNIT_ASSERT( forest->nodes()->at(9) == forest->active_node(0) ||
+ CPPUNIT_ASSERT( forest->nodes()->at(9) == forest->active_node(0) ||
forest->nodes()->at(10) == forest->active_node(0) );
CPPUNIT_ASSERT( forest->active_node(0)->local() );
@@ -677,7 +676,7 @@ class TestForest : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( (size_t)1, forest->active_node(0)->countChildren() );
Node* single_branch = forest->nodes()->at(9);
- if( forest->nodes()->at(9) == forest->active_node(0) )
+ if( forest->nodes()->at(9) == forest->active_node(0) )
single_branch = forest->nodes()->at(10);
CPPUNIT_ASSERT( !single_branch->local() );
@@ -687,7 +686,7 @@ class TestForest : public CppUnit::TestCase {
void testImplementCoalescence() {
for (size_t i=0; i<100; ++i) {
- forest->createExampleTree();
+ forest->createExampleTree();
Node* new_root = forest->cut(TreePoint(forest->nodes()->at(1), 1.5, false));
TimeIntervalIterator tii(forest, new_root);
@@ -709,7 +708,7 @@ class TestForest : public CppUnit::TestCase {
if ( !new_root->local() ) {
CPPUNIT_ASSERT( new_root->countChildren() == 2 );
- CPPUNIT_ASSERT( !(new_root->first_child()->local() && new_root->second_child()->local()) );
+ CPPUNIT_ASSERT( !(new_root->first_child()->local() && new_root->second_child()->local()) );
Node* child = new_root->first_child();
if (!new_root->second_child()->local()) child = new_root->second_child();
CPPUNIT_ASSERT( child->last_update() == 1 );
@@ -749,10 +748,11 @@ class TestForest : public CppUnit::TestCase {
// Circe detection
model2->addSingleMigrationEvent(0.5, 2, 0, 1.0);
- CPPUNIT_ASSERT_THROW( forest2->implementFixedTimeEvent(tii),
+ CPPUNIT_ASSERT_THROW( forest2->implementFixedTimeEvent(tii),
std::logic_error );
- delete forest2, model2;
+ delete forest2;
+ delete model2;
}
void testSamplePoint() {
@@ -760,7 +760,7 @@ class TestForest : public CppUnit::TestCase {
forest->createScaledExampleTree();
TreePoint point;
- int n0 = 0, n1 = 0, n2 = 0, n3 = 0,
+ int n0 = 0, n1 = 0, n2 = 0, n3 = 0,
n4 = 0, n5 = 0;
for (int i = 0; i < 240000; ++i) {
point = forest->samplePoint();
@@ -776,27 +776,38 @@ class TestForest : public CppUnit::TestCase {
//std::cout << n0 << " " << n1 << " " << n2 << " "
// << n3 << " " << n4 << " " << n5 << std::endl;
- CPPUNIT_ASSERT( 29500 <= n0 && n0 <= 30500 ); // expected 30000
- CPPUNIT_ASSERT( 29500 <= n1 && n1 <= 30500 ); // expected 30000
- CPPUNIT_ASSERT( 9800 <= n2 && n2 <= 10200 ); // expected 10000
- CPPUNIT_ASSERT( 9800 <= n3 && n3 <= 10200 ); // expected 10000
- CPPUNIT_ASSERT( 89000 <= n4 && n4 <= 91000 ); // expected 90000
- CPPUNIT_ASSERT( 69000 <= n5 && n5 <= 71000 ); // expected 70000
+ CPPUNIT_ASSERT( 29500 <= n0 && n0 <= 30500 ); // expected 30000
+ CPPUNIT_ASSERT( 29500 <= n1 && n1 <= 30500 ); // expected 30000
+ CPPUNIT_ASSERT( 9800 <= n2 && n2 <= 10200 ); // expected 10000
+ CPPUNIT_ASSERT( 9800 <= n3 && n3 <= 10200 ); // expected 10000
+ CPPUNIT_ASSERT( 89000 <= n4 && n4 <= 91000 ); // expected 90000
+ CPPUNIT_ASSERT( 69000 <= n5 && n5 <= 71000 ); // expected 70000
}
void testCopyConstructor() {
forest->createScaledExampleTree();
+ forest->set_next_base(10.0);
CPPUNIT_ASSERT( forest->coalescence_finished_ == true );
- Forest forest2 = Forest(*forest);
- CPPUNIT_ASSERT_EQUAL( forest->nodes()->size(), forest2.nodes()->size() );
- CPPUNIT_ASSERT( forest2.model_ == forest->model_ );
+ Forest forest2(*forest);
+
+ CPPUNIT_ASSERT_EQUAL(forest->nodes()->size(), forest2.nodes()->size() );
+ CPPUNIT_ASSERT(forest2.model_ == forest->model_);
+ CPPUNIT_ASSERT(forest2.local_root() != NULL);
+ CPPUNIT_ASSERT(forest2.local_root() != forest->local_root());
+
+ CPPUNIT_ASSERT_EQUAL(forest->rec_bases_.size(), forest2.rec_bases_.size());
+ CPPUNIT_ASSERT_EQUAL(forest->current_base(), forest2.current_base());
for (auto it = forest2.nodes()->iterator(); it.good(); ++it) {
CPPUNIT_ASSERT( (*it)->label() <= 4 );
CPPUNIT_ASSERT( (*it)->parent() != NULL || (*it)->first_child() != NULL );
}
+ CPPUNIT_ASSERT( forest2.checkLeafsOnLocalTree() );
+ CPPUNIT_ASSERT( forest2.checkTree() );
+
+ forest2.sampleNextGenealogy();
CPPUNIT_ASSERT( forest2.checkTree() );
}
@@ -817,9 +828,9 @@ class TestForest : public CppUnit::TestCase {
void testPrintLocusSumStats() {
ostringstream output;
- forest->printLocusSumStats(output);
-
forest->writable_model()->addSummaryStatistic(std::make_shared<TMRCA>());
+ forest->set_current_base(0);
+ forest->set_next_base(forest->model().loci_length());
forest->calcSegmentSumStats();
forest->printLocusSumStats(output);
CPPUNIT_ASSERT( output.str() != "" );
diff --git a/tests/unittests/test_model.cc b/tests/unittests/test_model.cc
index dca3442..3cc7c33 100644
--- a/tests/unittests/test_model.cc
+++ b/tests/unittests/test_model.cc
@@ -43,41 +43,41 @@ class TestModel : public CppUnit::TestCase {
public:
void testAddChangeTime() {
Model model = Model();
- std::vector<double> v1 = std::vector<double>(1, 1),
- v2 = std::vector<double>(1, 2),
- v3 = std::vector<double>(1, 3);
+ std::vector<double> v1 = std::vector<double>(1, 1),
+ v2 = std::vector<double>(1, 2),
+ v3 = std::vector<double>(1, 3);
// Check basic adding first time;
CPPUNIT_ASSERT( model.addChangeTime(0) == 0 );
CPPUNIT_ASSERT( model.change_times_.size() == 1 );
- CPPUNIT_ASSERT( model.change_times_.at(0) == 0 );
+ CPPUNIT_ASSERT( model.change_times_.at(0) == 0 );
CPPUNIT_ASSERT( model.pop_sizes_list_.size() == 1 );
model.pop_sizes_list_[0] = v1;
// Check adding a time at the end;
CPPUNIT_ASSERT_EQUAL( (size_t)1, model.addChangeTime(3) );
CPPUNIT_ASSERT( model.change_times_.size() == 2 );
- CPPUNIT_ASSERT( model.change_times_.at(1) == 3 );
- CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 );
- CPPUNIT_ASSERT( model.pop_sizes_list_[1].empty() );
+ CPPUNIT_ASSERT( model.change_times_.at(1) == 3 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[1].empty() );
model.pop_sizes_list_[1] = v3;
// Check adding a time in the middle;
CPPUNIT_ASSERT_EQUAL( (size_t)1, model.addChangeTime(2) );
CPPUNIT_ASSERT( model.change_times_.size() == 3 );
- CPPUNIT_ASSERT( model.change_times_.at(0) == 0 );
- CPPUNIT_ASSERT( model.change_times_.at(1) == 2 );
- CPPUNIT_ASSERT( model.change_times_.at(2) == 3 );
- CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 );
- CPPUNIT_ASSERT( model.pop_sizes_list_[1].empty() );
- CPPUNIT_ASSERT( model.pop_sizes_list_[2] == v3 );
+ CPPUNIT_ASSERT( model.change_times_.at(0) == 0 );
+ CPPUNIT_ASSERT( model.change_times_.at(1) == 2 );
+ CPPUNIT_ASSERT( model.change_times_.at(2) == 3 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[1].empty() );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[2] == v3 );
model.pop_sizes_list_[1] = v2;
- CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 );
- CPPUNIT_ASSERT( model.pop_sizes_list_[1] == v2 );
- CPPUNIT_ASSERT( model.pop_sizes_list_[2] == v3 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[0] == v1 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[1] == v2 );
+ CPPUNIT_ASSERT( model.pop_sizes_list_[2] == v3 );
- // Check that we don't add a time twice
+ // Check that we don't add a time twice
CPPUNIT_ASSERT_EQUAL( (size_t)0, model.addChangeTime(0) );
CPPUNIT_ASSERT_EQUAL( (size_t)1, model.addChangeTime(1) );
CPPUNIT_ASSERT_EQUAL( (size_t)2, model.addChangeTime(2) );
@@ -89,56 +89,56 @@ class TestModel : public CppUnit::TestCase {
void testAddSampleSizes() {
Model model = Model();
model.addSampleSizes(0.0, std::vector<size_t>(1,3));
- CPPUNIT_ASSERT_EQUAL( model.sample_size(), (size_t)3 );
- CPPUNIT_ASSERT_EQUAL( model.sample_population(2), (size_t)0 );
- CPPUNIT_ASSERT_EQUAL( model.sample_time(1), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_size(), (size_t)3 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_population(2), (size_t)0 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_time(1), (double)0.0 );
model = Model(0);
model.set_population_number(3);
model.addSampleSizes(0.0, std::vector<size_t>(3,1));
- CPPUNIT_ASSERT_EQUAL( model.sample_size(), (size_t)3 );
- CPPUNIT_ASSERT_EQUAL( model.sample_population(0), (size_t)0 );
- CPPUNIT_ASSERT_EQUAL( model.sample_population(1), (size_t)1 );
- CPPUNIT_ASSERT_EQUAL( model.sample_population(2), (size_t)2 );
- CPPUNIT_ASSERT_EQUAL( model.sample_time(0), (double)0.0 );
- CPPUNIT_ASSERT_EQUAL( model.sample_time(1), (double)0.0 );
- CPPUNIT_ASSERT_EQUAL( model.sample_time(2), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_size(), (size_t)3 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_population(0), (size_t)0 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_population(1), (size_t)1 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_population(2), (size_t)2 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_time(0), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_time(1), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model.sample_time(2), (double)0.0 );
std::vector<size_t> sample_size;
sample_size.push_back(0);
sample_size.push_back(0);
sample_size.push_back(2);
model.addSampleSizes(1.0, sample_size);
- CPPUNIT_ASSERT_EQUAL( (size_t)5, model.sample_size() );
- CPPUNIT_ASSERT_EQUAL( (size_t)2, model.sample_population(3) );
- CPPUNIT_ASSERT_EQUAL( (size_t)2, model.sample_population(4) );
- CPPUNIT_ASSERT_EQUAL( 1.0, model.sample_time(3) );
- CPPUNIT_ASSERT_EQUAL( 1.0, model.sample_time(4) );
+ CPPUNIT_ASSERT_EQUAL( (size_t)5, model.sample_size() );
+ CPPUNIT_ASSERT_EQUAL( (size_t)2, model.sample_population(3) );
+ CPPUNIT_ASSERT_EQUAL( (size_t)2, model.sample_population(4) );
+ CPPUNIT_ASSERT_EQUAL( 1.0, model.sample_time(3) );
+ CPPUNIT_ASSERT_EQUAL( 1.0, model.sample_time(4) );
Model model2 = Model();
std::vector<size_t> sample_sizes;
sample_sizes.push_back(5);
sample_sizes.push_back(2);
model2.addSampleSizes(0.0, sample_sizes);
- CPPUNIT_ASSERT_EQUAL( model2.sample_size(), (size_t)7 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(0), (size_t)0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(4), (size_t)0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(5), (size_t)1 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(6), (size_t)1 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_time(0), (double)0.0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_time(4), (double)0.0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_time(6), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_size(), (size_t)7 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(0), (size_t)0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(4), (size_t)0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(5), (size_t)1 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(6), (size_t)1 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_time(0), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_time(4), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_time(6), (double)0.0 );
sample_sizes.clear();
sample_sizes.push_back(2);
sample_sizes.push_back(1);
model2.addSampleSizes(7.4, sample_sizes);
- CPPUNIT_ASSERT_EQUAL( model2.sample_size(), (size_t)10 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(7), (size_t)0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(8), (size_t)0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_population(9), (size_t)1 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_time(0), (double)0.0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_time(6), (double)0.0 );
- CPPUNIT_ASSERT_EQUAL( model2.sample_time(8), (double)7.4 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_size(), (size_t)10 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(7), (size_t)0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(8), (size_t)0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_population(9), (size_t)1 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_time(0), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_time(6), (double)0.0 );
+ CPPUNIT_ASSERT_EQUAL( model2.sample_time(8), (double)7.4 );
}
void testAddPopulationSizes() {
@@ -153,25 +153,25 @@ class TestModel : public CppUnit::TestCase {
model.addPopulationSizes(1, std::vector<double>(2, 5), true, false);
model.increaseTime();
CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size, model.getCurrentTime() );
- CPPUNIT_ASSERT( model.population_size(0) == 5 );
+ CPPUNIT_ASSERT( areSame(model.population_size(0), 5) );
model.addPopulationSizes(2, 10, true, false);
model.increaseTime();
CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size, model.getCurrentTime() );
- CPPUNIT_ASSERT( model.population_size(0) == 10 );
+ CPPUNIT_ASSERT( areSame(model.population_size(0), 10) );
auto pop_sizes2 = std::vector<double>();
pop_sizes2.push_back(7);
pop_sizes2.push_back(6);
model.addPopulationSizes(0.0, pop_sizes2);
model.resetTime();
- CPPUNIT_ASSERT( model.population_size(0) == 7 );
- CPPUNIT_ASSERT( model.population_size(1) == 6 );
-
+ CPPUNIT_ASSERT( areSame(model.population_size(0), 7) );
+ CPPUNIT_ASSERT( areSame(model.population_size(1), 6) );
+
CPPUNIT_ASSERT_THROW( model.addPopulationSizes(1, std::vector<double>(1, 5)), std::logic_error );
CPPUNIT_ASSERT_THROW( model.addPopulationSizes(1, std::vector<double>(3, 5)), std::logic_error );
}
-
+
void testAddRelativePopulationSizes() {
Model model = Model();
model.set_population_number(2);
@@ -241,7 +241,7 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT_THROW( model.addMigrationRates(1, rates), std::logic_error );
rates.push_back(9);
model.addMigrationRates(1, rates);
-
+
model.resetTime();
model.increaseTime();
CPPUNIT_ASSERT_EQUAL( 2.0, model.migration_rate(0,1) );
@@ -314,13 +314,17 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL(1.0, model.change_position_[2]);
CPPUNIT_ASSERT_EQUAL((size_t)3, model.mutation_rates_.size());
CPPUNIT_ASSERT_EQUAL((size_t)3, model.recombination_rates_.size());
+
+ model.setLocusLength(10.0);
+ CPPUNIT_ASSERT_THROW(model.addChangePosition(-0.5), std::invalid_argument);
+ CPPUNIT_ASSERT_THROW(model.addChangePosition(10.5), std::invalid_argument);
}
void testIncreaseTime() {
Model model = Model(7);
model.addGrowthRates(1.0, std::vector<double>(1, 1.5));
model.addGrowthRates(2.0, std::vector<double>(1, 1));
-
+
CPPUNIT_ASSERT_EQUAL( (size_t)0, model.current_time_idx_ );
CPPUNIT_ASSERT_NO_THROW( model.increaseTime() );
@@ -390,20 +394,21 @@ class TestModel : public CppUnit::TestCase {
void testSetGetMutationRate() {
Model model = Model(5);
+ model.setLocusLength(10);
model.setMutationRate(0.001);
CPPUNIT_ASSERT_EQUAL( 0.001, model.mutation_rate() );
-
+
model.setMutationRate(10, false, false);
CPPUNIT_ASSERT_EQUAL( 10.0, model.mutation_rate() );
model.setMutationRate(10, true, false);
- CPPUNIT_ASSERT_EQUAL( 10.0/model.default_loci_length, model.mutation_rate() );
+ CPPUNIT_ASSERT_EQUAL( 1.0, model.mutation_rate() );
model.setMutationRate(10, false, true);
CPPUNIT_ASSERT_EQUAL( 10.0/(4*model.default_pop_size), model.mutation_rate() );
model.setMutationRate(10, true, true);
- CPPUNIT_ASSERT_EQUAL(10.0/(4*model.default_pop_size*model.default_loci_length),
+ CPPUNIT_ASSERT_EQUAL(10.0/(4*model.default_pop_size*10),
model.mutation_rate() );
// Test setting multiple rates
@@ -462,8 +467,6 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT( areSame(0.00001/(4*model.default_pop_size), model.recombination_rate()) );
CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() );
- CPPUNIT_ASSERT_THROW( model.setRecombinationRate(-0.001, 100), std::invalid_argument );
-
// Test setting multiple rates
model.setRecombinationRate(10, false, false, 0.0);
model.setRecombinationRate(5, false, false, 2.0);
@@ -478,6 +481,12 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL(7.5, model.recombination_rate());
model.increaseSequencePosition();
CPPUNIT_ASSERT_EQUAL(5.0, model.recombination_rate());
+
+ // Check for error on invalid change positions
+ CPPUNIT_ASSERT_THROW( model.setRecombinationRate(-0.001, 100), std::invalid_argument);
+ model.setLocusLength(10);
+ CPPUNIT_ASSERT_THROW(model.setRecombinationRate(7.5, false, false, -1.0), std::invalid_argument);
+ CPPUNIT_ASSERT_THROW(model.setRecombinationRate(7.5, false, false, 11.0), std::invalid_argument);
}
void testCheck() {
@@ -496,11 +505,11 @@ class TestModel : public CppUnit::TestCase {
model.set_population_number(2);
model.addSingleMigrationEvent(1, 0, 1, 1);
CPPUNIT_ASSERT_NO_THROW( model.check() );
- }
+ }
void testPopSizeAfterGrowth() {
// Growth only
- Model model = Model(5);
+ Model model = Model(5);
model.set_population_number(2);
model.addSymmetricMigration(0, 1.0);
model.addPopulationSizes(0, 1000);
@@ -523,9 +532,9 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT( areSame( n_2*std::exp(-2), model.population_size(1) ) );
CPPUNIT_ASSERT( areSame( n_1*std::exp(-2*2), model.population_size(0, 4.5) ) );
CPPUNIT_ASSERT( areSame( n_2*std::exp(-2*2), model.population_size(1, 4.5) ) );
-
+
// Growth with a pop size change in between
- model = Model(5);
+ model = Model(5);
model.set_population_number(2);
model.addSymmetricMigration(0, 1.0);
growth_rates.clear();
@@ -539,7 +548,7 @@ class TestModel : public CppUnit::TestCase {
model.addGrowthRate(3.5, 0, 2.0);
model.finalize();
//std::cout << model << std::endl;
-
+
model.resetTime();
double size0 = model.default_pop_size;
@@ -547,12 +556,12 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) );
CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) );
- model.increaseTime();
+ model.increaseTime();
CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() );
CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) );
CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) );
- model.increaseTime();
+ model.increaseTime();
CPPUNIT_ASSERT_EQUAL( 1.25, model.getCurrentTime() );
CPPUNIT_ASSERT( areSame( size0*std::exp(0.5*0.25), model.population_size(0) ) );
CPPUNIT_ASSERT( areSame( size0*std::exp(-0.5*0.25), model.population_size(1) ) );
@@ -571,7 +580,7 @@ class TestModel : public CppUnit::TestCase {
model.increaseTime();
CPPUNIT_ASSERT_EQUAL( 3.0, model.getCurrentTime() );
- size0 = 500;
+ size0 = 500;
size1 *= exp(-2*0.5);
CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) );
CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) );
@@ -583,7 +592,7 @@ class TestModel : public CppUnit::TestCase {
CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) );
CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) );
- model = Model(5);
+ model = Model(5);
model.set_population_number(2);
model.addSymmetricMigration(0, 1.0);
model.addPopulationSize(0, 1, 500);
@@ -596,16 +605,17 @@ class TestModel : public CppUnit::TestCase {
void testAddSummaryStatistic() {
Model model = Model(5);
- CPPUNIT_ASSERT(model.countSummaryStatistics() == 0);
+ CPPUNIT_ASSERT(model.countSummaryStatistics() == 0);
model.addSummaryStatistic(std::make_shared<TMRCA>());
- CPPUNIT_ASSERT(model.countSummaryStatistics() == 1);
+ CPPUNIT_ASSERT(model.countSummaryStatistics() == 1);
}
void testSetLocusLength() {
Model model = Model(5);
- CPPUNIT_ASSERT_EQUAL(model.default_loci_length, model.loci_length() );
+ model.setLocusLength(10);
+ CPPUNIT_ASSERT_EQUAL((size_t)10, model.loci_length() );
model.setMutationRate(5.0, true, false);
- CPPUNIT_ASSERT_EQUAL(5.0/model.default_loci_length, model.mutation_rate());
+ CPPUNIT_ASSERT_EQUAL(0.5, model.mutation_rate());
model.setLocusLength(2000);
CPPUNIT_ASSERT_EQUAL((size_t)2000, model.loci_length() );
@@ -619,9 +629,9 @@ class TestModel : public CppUnit::TestCase {
void testAddPopToVectorList() {
Model model = Model(5);
std::vector<std::vector<double> > vector_list;
- vector_list.push_back(std::vector<double>(2, 1.0));
- vector_list.push_back(std::vector<double>());
- vector_list.push_back(std::vector<double>(2, 0.5));
+ vector_list.push_back(std::vector<double>(2, 1.0));
+ vector_list.push_back(std::vector<double>());
+ vector_list.push_back(std::vector<double>(2, 0.5));
model.addPopToVectorList(vector_list);
CPPUNIT_ASSERT_EQUAL( (size_t)3, vector_list.at(0).size() );
@@ -638,9 +648,9 @@ class TestModel : public CppUnit::TestCase {
void testAddPopToMatrixList() {
Model model = Model(5);
std::vector<std::vector<double> > vector_list;
- vector_list.push_back(std::vector<double>(2, 1.0));
- vector_list.push_back(std::vector<double>());
- vector_list.push_back(std::vector<double>(2, 0.5));
+ vector_list.push_back(std::vector<double>(2, 1.0));
+ vector_list.push_back(std::vector<double>());
+ vector_list.push_back(std::vector<double>(2, 0.5));
model.set_population_number(3);
model.addPopToMatrixList(vector_list, 2);
diff --git a/tests/unittests/test_node.cc b/tests/unittests/test_node.cc
index 1e6f501..9a80ed4 100644
--- a/tests/unittests/test_node.cc
+++ b/tests/unittests/test_node.cc
@@ -145,12 +145,12 @@ class TestNode : public CppUnit::TestCase {
}
void testLocalNavigation() {
- Node* n1 = new Node(7.5);
+ Node* n1 = forest->nodes()->createNode(7.5);
n1->make_local();
- Node* n2 = new Node(6.5);
+ Node* n2 = forest->nodes()->createNode(6.5);
n2->make_nonlocal(1.0);
- Node *n3 = forest->nodes()->at(4),
+ Node *n3 = forest->nodes()->at(4),
*root = forest->local_root(),
*n4 = forest->nodes()->at(5);
@@ -163,7 +163,7 @@ class TestNode : public CppUnit::TestCase {
CPPUNIT_ASSERT(n1->getLocalParent() == root);
CPPUNIT_ASSERT(n3->getLocalParent() == root);
- CPPUNIT_ASSERT( (root->getLocalChild1() == n3 && root->getLocalChild2() == n4) ||
+ CPPUNIT_ASSERT( (root->getLocalChild1() == n3 && root->getLocalChild2() == n4) ||
(root->getLocalChild1() == n4 && root->getLocalChild2() == n3) );
CPPUNIT_ASSERT(forest->nodes()->at(0)->getLocalChild1() == NULL);
diff --git a/tests/unittests/test_node_container.cc b/tests/unittests/test_node_container.cc
index fc082cd..658d58c 100644
--- a/tests/unittests/test_node_container.cc
+++ b/tests/unittests/test_node_container.cc
@@ -33,10 +33,10 @@ class TestNodeContainer : public CppUnit::TestCase {
public:
void setUp() {
- node1 = new Node(1, 1);
- node2 = new Node(2, 2);
- node3 = new Node(3, 0);
nc = NodeContainer();
+ node1 = nc.createNode(1, 1);
+ node2 = nc.createNode(2, 2);
+ node3 = nc.createNode(3, 0);
nc.add(node1);
nc.add(node2);
nc.add(node3);
@@ -61,28 +61,28 @@ class TestNodeContainer : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( (size_t)3, nc.size() );
// Add new first node
- Node* node = new Node(0);
+ Node* node = nc.createNode(0);
nc.add(node);
CPPUNIT_ASSERT_EQUAL( (size_t)4, nc.size() );
CPPUNIT_ASSERT( nc.get(0) == node );
CPPUNIT_ASSERT( nc.first() == node );
// Add new last node
- node = new Node(10);
+ node = nc.createNode(10);
nc.add(node);
CPPUNIT_ASSERT_EQUAL( (size_t)5, nc.size() );
CPPUNIT_ASSERT( nc.get(4) == node );
CPPUNIT_ASSERT( nc.last() == node );
- // Add something in between
- node = new Node(3);
+ // Add something in between
+ node = nc.createNode(3);
nc.add(node);
CPPUNIT_ASSERT_EQUAL( (size_t)6, nc.size() );
CPPUNIT_ASSERT( nc.get(4) == node );
}
void testRemove() {
- Node *node = new Node(2.5);
+ Node *node = nc.createNode(2.5);
nc.add(node);
nc.remove(node);
CPPUNIT_ASSERT( nc.size() == 3 );
@@ -106,7 +106,7 @@ class TestNodeContainer : public CppUnit::TestCase {
void testMove() {
- Node* node4 = new Node(4);
+ Node* node4 = nc.createNode(4);
nc.add(node4);
// Middle -> Middle (forwards)
@@ -134,7 +134,7 @@ class TestNodeContainer : public CppUnit::TestCase {
CPPUNIT_ASSERT( nc.first() == node2 );
CPPUNIT_ASSERT( nc.last() == node1 );
CPPUNIT_ASSERT( node1->height() == 10 );
-
+
// end -> start
nc.move(node1, 1);
CPPUNIT_ASSERT( nc.get(0) == node1 );
@@ -193,7 +193,7 @@ class TestNodeContainer : public CppUnit::TestCase {
--it; --it;
CPPUNIT_ASSERT_THROW( --it, std::out_of_range );
}
-
+
void testReverseIterator() {
ReverseConstNodeIterator it = nc.reverse_iterator();
@@ -205,7 +205,7 @@ class TestNodeContainer : public CppUnit::TestCase {
CPPUNIT_ASSERT( !it.good() );
CPPUNIT_ASSERT_THROW( ++it, std::out_of_range );
}
-
+
void testNodeIteratorHeight() {
NodeIterator it = nc.iterator();
diff --git a/tests/unittests/test_param.cc b/tests/unittests/test_param.cc
index 4c27f52..5b176df 100644
--- a/tests/unittests/test_param.cc
+++ b/tests/unittests/test_param.cc
@@ -27,6 +27,7 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_TEST( testErrorOnInfintieRecRate );
CPPUNIT_TEST( testScientificNotation );
CPPUNIT_TEST( testApproximation );
+ CPPUNIT_TEST( testTransposeSegSites );
CPPUNIT_TEST_SUITE_END();
@@ -53,40 +54,40 @@ class TestParam : public CppUnit::TestCase {
void testParseBasicCmdStructure() {
Model model;
- CPPUNIT_ASSERT_THROW(Param("scrm").parse(), std::invalid_argument);
+ CPPUNIT_ASSERT_THROW(Param("scrm").parse(), std::invalid_argument);
Param pars = Param("-h");
model = pars.parse();
- CPPUNIT_ASSERT( pars.help() );
+ CPPUNIT_ASSERT( pars.help() );
pars = Param("--help");
model = pars.parse();
- CPPUNIT_ASSERT( pars.help() );
+ CPPUNIT_ASSERT( pars.help() );
pars = Param("2 1 -h");
model = pars.parse();
- CPPUNIT_ASSERT(pars.help());
+ CPPUNIT_ASSERT(pars.help());
- CPPUNIT_ASSERT_THROW( Param("1").parse(), std::invalid_argument );
- CPPUNIT_ASSERT_THROW( Param("1 -t 5").parse(), std::invalid_argument );
- CPPUNIT_ASSERT_THROW( Param("-t 5").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param("1").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param("1 -t 5").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param("-t 5").parse(), std::invalid_argument );
pars = Param("2 1");
model = pars.parse();
- CPPUNIT_ASSERT( !pars.help() );
- CPPUNIT_ASSERT( !pars.version() );
+ CPPUNIT_ASSERT( !pars.help() );
+ CPPUNIT_ASSERT( !pars.version() );
pars = Param("-v");
model = pars.parse();
- CPPUNIT_ASSERT( pars.version() );
+ CPPUNIT_ASSERT( pars.version() );
pars = Param("--version");
model = pars.parse();
- CPPUNIT_ASSERT( pars.version() );
+ CPPUNIT_ASSERT( pars.version() );
pars = Param("2 1 -v");
model = pars.parse();
- CPPUNIT_ASSERT( pars.version() );
+ CPPUNIT_ASSERT( pars.version() );
}
void testParseSeeds() {
@@ -102,8 +103,8 @@ class TestParam : public CppUnit::TestCase {
model = pars.parse();
CPPUNIT_ASSERT_EQUAL(seed, pars.random_seed());
- CPPUNIT_ASSERT_THROW(Param("20 10 -seed 1 2 3 4").parse(), std::invalid_argument);
- CPPUNIT_ASSERT_THROW(Param("20 10 -seed -t 2").parse(), std::invalid_argument);
+ CPPUNIT_ASSERT_THROW(Param("20 10 -seed 1 2 3 4").parse(), std::invalid_argument);
+ CPPUNIT_ASSERT_THROW(Param("20 10 -seed -t 2").parse(), std::invalid_argument);
}
void testParseCommonOptions() {
@@ -117,9 +118,9 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT( model.has_window_seq() );
CPPUNIT_ASSERT_EQUAL( 1000.0, model.window_length_seq() );
- CPPUNIT_ASSERT_THROW( Param("15 10 -I 3 7 8 5").parse(), std::invalid_argument );
- CPPUNIT_ASSERT_THROW( Param("15 10 -tv").parse(), std::invalid_argument );
- CPPUNIT_ASSERT_THROW( Param("15 10 -t 17 1").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param("15 10 -I 3 7 8 5").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param("15 10 -tv").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param("15 10 -t 17 1").parse(), std::invalid_argument );
model = Param("20 10 -t 3.74 -I 3 7 8 5 -T -M 5.0").parse();
CPPUNIT_ASSERT_EQUAL( (size_t)3, model.population_number() );
@@ -156,7 +157,7 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( 12.3 * 4 * model.default_pop_size, model.sample_time(20) );
CPPUNIT_ASSERT_EQUAL( 12.3 * 4 * model.default_pop_size, model.sample_time(22) );
- // -N & -eN
+ // -N & -eN
model = Param("20 10 -t 3.74 -I 3 7 8 5 -N 0.3 -eN 8.2 0.75 -G 1.5 -M 5.0").parse();
model.resetTime();
CPPUNIT_ASSERT_EQUAL( 0.0, model.getCurrentTime() );
@@ -172,14 +173,14 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(1) );
CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(2) );
- // -n & -en
+ // -n & -en
model = Param("20 10 -t 3.74 -I 3 7 8 5 -G 1.5 -n 2 0.3 -eN 1.1 0.75 -en 2 3 0.1 -eG 1.5 2 -M 5.0").parse();
model.finalize();
model.resetTime();
CPPUNIT_ASSERT_EQUAL( 0.0, model.getCurrentTime() );
CPPUNIT_ASSERT( areSame(model.default_pop_size, model.population_size(0)) );
CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size, model.population_size(1)) );
- CPPUNIT_ASSERT( areSame(model.default_pop_size, (size_t)model.population_size(2)) );
+ CPPUNIT_ASSERT( areSame(model.default_pop_size, model.population_size(2)) );
CPPUNIT_ASSERT( areSame(1.5 / 4 / model.default_pop_size, model.growth_rate(0)) );
CPPUNIT_ASSERT( areSame(1.5 / 4 / model.default_pop_size, model.growth_rate(1)) );
CPPUNIT_ASSERT( areSame(1.5 / 4 / model.default_pop_size, model.growth_rate(2)) );
@@ -191,7 +192,7 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT( areSame(2.0 * 4 * model.default_pop_size, model.getCurrentTime()) );
CPPUNIT_ASSERT( 0.75*model.default_pop_size > model.population_size(0) );
CPPUNIT_ASSERT( 0.75*model.default_pop_size > model.population_size(1) );
- CPPUNIT_ASSERT_EQUAL( (size_t)(0.10*model.default_pop_size), (size_t)model.population_size(2) );
+ CPPUNIT_ASSERT_EQUAL( 0.10 * model.default_pop_size, model.population_size(2) );
CPPUNIT_ASSERT_EQUAL( 2.0 / 4 / model.default_pop_size , model.growth_rate(0) );
CPPUNIT_ASSERT_EQUAL( 2.0 / 4 / model.default_pop_size, model.growth_rate(1) );
CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(2) );
@@ -277,13 +278,13 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT_EQUAL( model.default_pop_size, model.population_size(1) );
CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate() );
- CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -eM 0.9 5.0").parse(),
+ CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -eM 0.9 5.0").parse(),
std::invalid_argument );
- CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -eG 0.9 5.0").parse(),
+ CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -eG 0.9 5.0").parse(),
std::invalid_argument );
- CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -eN 0.9 5.0").parse(),
+ CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -eN 0.9 5.0").parse(),
std::invalid_argument );
- CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -es 0.9 3 1").parse(),
+ CPPUNIT_ASSERT_THROW(Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5 -es 0.9 3 1").parse(),
std::invalid_argument );
}
@@ -315,7 +316,7 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT_NO_THROW(model = Param("20 10 -L").parse());
CPPUNIT_ASSERT( model.summary_statistics_.size() == 1 );
- CPPUNIT_ASSERT_THROW(Param("20 10 -oSFS").parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW(Param("20 10 -oSFS").parse(), std::invalid_argument );
CPPUNIT_ASSERT_NO_THROW(model = Param("20 10 -oSFS -t 5").parse());
CPPUNIT_ASSERT( model.summary_statistics_.size() == 2 );
@@ -330,7 +331,7 @@ class TestParam : public CppUnit::TestCase {
void testParseGrowthOptions() {
// -G && -eG
Model model;
- char *argv[] = { "scrm", "20", "10", "-t", "3.74", "-I", "2", "10", "10",
+ char *argv[] = { "scrm", "20", "10", "-t", "3.74", "-I", "2", "10", "10",
"-G", "2", "-eG", "1", "3", "-M", "5.0" };
CPPUNIT_ASSERT_NO_THROW( model = Param(16, argv).parse(); );
model.resetTime();
@@ -342,8 +343,8 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size, model.growth_rate(1)) );
// -g && -eg
- char *argv2[] = {
- "scrm", "20", "10", "-t", "3.74", "-I", "2", "10", "10",
+ char *argv2[] = {
+ "scrm", "20", "10", "-t", "3.74", "-I", "2", "10", "10",
"-g", "2", "0.1", "-eG", "1", "3", "-eg", "2", "1", "2.4", "-M", "5.0"};
CPPUNIT_ASSERT_NO_THROW( model = Param(21, argv2).parse(); );
model.resetTime();
@@ -360,21 +361,21 @@ class TestParam : public CppUnit::TestCase {
void testParseVariableRates() {
// -st
Model model;
- char *argv[] = { "scrm", "20", "10", "-t", "3.74", "-st", "10", "5.1", "-st", "4.5", "3.2"};
+ char *argv[] = { "scrm", "20", "10", "-t", "3.74", "-st", "0.7", "5.1", "-st", "0.5", "3.2"};
CPPUNIT_ASSERT_NO_THROW( model = Param(11, argv).parse(); );
- double scale = 1 / ( 4 * model.default_pop_size * model.default_loci_length );
+ double scale = 1 / ( 4 * model.default_pop_size );
CPPUNIT_ASSERT( areSame(3.74 * scale, model.mutation_rate()) );
CPPUNIT_ASSERT_EQUAL( 0.0, model.getCurrentSequencePosition() );
- CPPUNIT_ASSERT_EQUAL( 4.5, model.getNextSequencePosition() );
+ CPPUNIT_ASSERT_EQUAL( 0.5, model.getNextSequencePosition() );
model.increaseSequencePosition();
CPPUNIT_ASSERT( areSame(3.2 * scale, model.mutation_rate()) );
- CPPUNIT_ASSERT_EQUAL( 4.5, model.getCurrentSequencePosition() );
- CPPUNIT_ASSERT_EQUAL( 10.0, model.getNextSequencePosition() );
+ CPPUNIT_ASSERT_EQUAL( 0.5, model.getCurrentSequencePosition() );
+ CPPUNIT_ASSERT_EQUAL( 0.7, model.getNextSequencePosition() );
model.increaseSequencePosition();
CPPUNIT_ASSERT( areSame(5.1 * scale, model.mutation_rate()) );
- CPPUNIT_ASSERT_EQUAL( 10.0, model.getCurrentSequencePosition() );
+ CPPUNIT_ASSERT_EQUAL( 0.7, model.getCurrentSequencePosition() );
// -sr
char *argv2[] = { "scrm", "20", "10", "-r", "3.74", "100", "-sr", "10", "5.1", "-sr", "4.5", "3.2"};
@@ -421,10 +422,10 @@ class TestParam : public CppUnit::TestCase {
void testParseUnsupportedMsArguments() {
char *argv[] = { "scrm", "20", "10", "-t", "3.74", "-c", "10", "5.1", "-st", "4.5", "3.2"};
- CPPUNIT_ASSERT_THROW( Param(11, argv).parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param(11, argv).parse(), std::invalid_argument );
char *argv2[] = { "scrm", "20", "10", "-s", "5"};
- CPPUNIT_ASSERT_THROW( Param(5, argv2).parse(), std::invalid_argument );
+ CPPUNIT_ASSERT_THROW( Param(5, argv2).parse(), std::invalid_argument );
}
@@ -510,6 +511,20 @@ class TestParam : public CppUnit::TestCase {
CPPUNIT_ASSERT( !model.has_window_seq() );
CPPUNIT_ASSERT_EQUAL( (size_t)10, model.window_length_rec() );
}
+
+ void testTransposeSegSites() {
+ Model model = Param("2 2 -r 10 100 -t 5 -transpose-segsites").parse();
+ SegSites* ss = dynamic_cast<SegSites*>(model.getSummaryStatistic(0));
+ CPPUNIT_ASSERT(ss->get_transpose());
+
+ model = Param("2 2 -t 5 --transpose-segsites").parse();
+ ss = dynamic_cast<SegSites*>(model.getSummaryStatistic(0));
+ CPPUNIT_ASSERT(ss->get_transpose());
+
+ model = Param("2 2 -r 10 100 -t 5").parse();
+ ss = dynamic_cast<SegSites*>(model.getSummaryStatistic(0));
+ CPPUNIT_ASSERT(!ss->get_transpose());
+ }
};
//Uncomment this to activate the test
diff --git a/tests/unittests/test_random_generator.cc b/tests/unittests/test_random_generator.cc
index 1d65ab3..4527fd3 100644
--- a/tests/unittests/test_random_generator.cc
+++ b/tests/unittests/test_random_generator.cc
@@ -42,7 +42,7 @@ class TestRandomGenerator : public CppUnit::TestCase {
for (size_t i = 0; i < 1000; ++i) {
expo = rg->sampleExpoLimit(1, 2);
CPPUNIT_ASSERT( expo == -1 || expo > 0 );
- CPPUNIT_ASSERT( expo <= 2 );
+ CPPUNIT_ASSERT( expo <= 2 );
}
}
@@ -66,12 +66,12 @@ class TestRandomGenerator : public CppUnit::TestCase {
//std::cout << expo << std::endl;
CPPUNIT_ASSERT( 0.199 <= expo && expo <= 0.201 );
}
-
+
void testSampleExpoExpoLimit() {
size_t n = 10000;
double expo = 0, sample = 0;
size_t sample_number = 0;
-
+
// c = 0;
for (size_t i = 0; i < n; ++i) {
sample = rg->sampleExpoExpoLimit(2,0,1);
@@ -81,7 +81,7 @@ class TestRandomGenerator : public CppUnit::TestCase {
}
}
expo /= sample_number;
- // Expected: 0.34
+ // Expected: 0.34
CPPUNIT_ASSERT( 0.32 < expo && expo < 0.36 );
expo = 0;
@@ -94,7 +94,7 @@ class TestRandomGenerator : public CppUnit::TestCase {
}
}
expo /= sample_number;
- // Expected: 0.46
+ // Expected: 0.46
CPPUNIT_ASSERT( 0.44 < expo && expo < 0.48 );
expo = 0;
@@ -107,7 +107,7 @@ class TestRandomGenerator : public CppUnit::TestCase {
}
}
expo /= sample_number;
- // Expected: 0.50
+ // Expected: 0.50
CPPUNIT_ASSERT( 0.48 < expo && expo < 0.52 );
}
@@ -118,12 +118,12 @@ class TestRandomGenerator : public CppUnit::TestCase {
for (size_t i = 0; i < n; ++i) {
sample = rg->sampleInt(5);
- ++result[sample];
+ ++result[sample];
}
for (size_t i = 0; i < 5; ++i) {
//std::cout << i << " : " << result[i] << std::endl;
- CPPUNIT_ASSERT( 9750 < result[i] && result[i] < 10250 );
+ CPPUNIT_ASSERT( 9750 < result[i] && result[i] < 10250 );
}
}
diff --git a/tests/unittests/test_summary_statistics.cc b/tests/unittests/test_summary_statistics.cc
index bb9547e..da3d171 100644
--- a/tests/unittests/test_summary_statistics.cc
+++ b/tests/unittests/test_summary_statistics.cc
@@ -50,6 +50,7 @@ class TestSummaryStatistics : public CppUnit::TestCase {
void testTMRCA() {
forest->createScaledExampleTree();
+ forest->set_next_base(10);
TMRCA tmrca = TMRCA();
tmrca.calculate(*forest);
CPPUNIT_ASSERT_EQUAL( 10.0, tmrca.tmrca().at(0) );
@@ -58,7 +59,7 @@ class TestSummaryStatistics : public CppUnit::TestCase {
tmrca.calculate(*forest);
CPPUNIT_ASSERT_EQUAL( 10.0, tmrca.tmrca().at(1) );
CPPUNIT_ASSERT_EQUAL( 24.0, tmrca.tree_length().at(1) );
-
+
tmrca.clear();
CPPUNIT_ASSERT_EQUAL( (size_t)0, tmrca.tmrca().size() );
CPPUNIT_ASSERT_EQUAL( (size_t)0, tmrca.tree_length().size() );
@@ -164,13 +165,13 @@ class TestSummaryStatistics : public CppUnit::TestCase {
// Check distribution
seg_sites.calculate(*forest);
- CPPUNIT_ASSERT( seg_sites.countMutations() > 0 );
+ CPPUNIT_ASSERT( seg_sites.countMutations() > 0 );
int freqs[4] = { 0 };
int types[2] = { 0 };
int sum = 0;
for (auto it = seg_sites.haplotypes_.begin(); it != seg_sites.haplotypes_.end(); ++it) {
for (size_t i = 0; i < 4; ++i) {
- if ((*it)[i]) ++freqs[i];
+ if ((*it)[i]) ++freqs[i];
sum += (*it)[i];
}
CPPUNIT_ASSERT( sum == 1 || sum == 2 );
@@ -200,7 +201,7 @@ class TestSummaryStatistics : public CppUnit::TestCase {
forest->writable_model()->setMutationRate(0.0001);
forest->set_current_base(0.0);
forest->set_next_base(10.0);
-
+
std::shared_ptr<SegSites> seg_sites = std::make_shared<SegSites>();
seg_sites->positions_.push_back(0.5);
seg_sites->positions_.push_back(0.7);
@@ -248,16 +249,16 @@ class TestSummaryStatistics : public CppUnit::TestCase {
OrientedForest of(4);
size_t pos = 2*forest->sample_size()-2;
double sf = forest->model().scaling_factor();
- of.generateTreeData(forest->local_root(), pos, 0, sf);
+ of.generateTreeData(forest->local_root(), pos, 0, sf);
CPPUNIT_ASSERT( of.heights_.at(0) == 0.0 );
CPPUNIT_ASSERT( of.heights_.at(1) == 0.0 );
CPPUNIT_ASSERT( of.heights_.at(2) == 0.0 );
CPPUNIT_ASSERT( of.heights_.at(3) == 0.0 );
CPPUNIT_ASSERT( of.parents_.at(4) == 7 );
CPPUNIT_ASSERT( of.parents_.at(5) == 7 );
- CPPUNIT_ASSERT( of.heights_.at(6) == 10.0 * sf && of.parents_.at(6) == 0 );
+ CPPUNIT_ASSERT( areSame(of.heights_.at(6), 10.0 * sf) && of.parents_.at(6) == 0.0 );
- CPPUNIT_ASSERT( of.heights_.at(4) == 1.0 * sf || of.heights_.at(4) == 3.0 * sf );
+ CPPUNIT_ASSERT( areSame(of.heights_.at(4), 1.0 * sf) || areSame(of.heights_.at(4), 3.0 * sf) );
if ( of.heights_.at(4) == 1.0 ) {
CPPUNIT_ASSERT( of.parents_.at(0) == 5 );
CPPUNIT_ASSERT( of.parents_.at(1) == 5 );
@@ -276,7 +277,7 @@ class TestSummaryStatistics : public CppUnit::TestCase {
forest->createScaledExampleTree();
forest->set_current_base(0.0);
forest->set_next_base(10.0);
-
+
OrientedForest of(4);
std::ostringstream output;
of.calculate(*forest);
@@ -300,14 +301,14 @@ class TestSummaryStatistics : public CppUnit::TestCase {
forest->createScaledExampleTree();
forest->set_current_base(0.0);
forest->set_next_base(10.0);
-
+
NewickTree of(4, forest->model().has_recombination());
std::ostringstream output;
of.calculate(*forest);
of.printSegmentOutput(output);
//std::cout << output.str() << std::endl;
CPPUNIT_ASSERT( output.str().compare("((1:1,2:1):9,(3:3,4:3):7);\n") == 0 );
-
+
output.str("");
output.clear();
forest->writable_model()->setRecombinationRate(0.0001);
@@ -316,6 +317,30 @@ class TestSummaryStatistics : public CppUnit::TestCase {
of.printSegmentOutput(output);
//std::cout << output.str() << std::endl;
CPPUNIT_ASSERT( output.str().compare("[10]((1:1,2:1):9,(3:3,4:3):7);\n") == 0 );
+
+ /** Scientific notation is disabled when using absolute scaling **/
+ output.str("");
+ output.clear();
+ forest->createScaledExampleTree();
+ forest->set_current_base(0.0);
+ forest->set_next_base(10000000.0);
+ of.calculate(*forest);
+ of.printSegmentOutput(output);
+ //std::cout << output.str() << std::endl;
+ CPPUNIT_ASSERT( output.str().compare("[10000000]((1:1,2:1):9,(3:3,4:3):7);\n") == 0 );
+
+ /** Relative sequence scaling **/
+ forest->writable_model()->setSequenceScaling(relative);
+ forest->writable_model()->setLocusLength(10000);
+ forest->createScaledExampleTree();
+ forest->set_current_base(0.0);
+ forest->set_next_base(100);
+ output.str("");
+ output.clear();
+ of.calculate(*forest);
+ of.printSegmentOutput(output);
+ //std::cout << output.str() << std::endl;
+ CPPUNIT_ASSERT( output.str().compare("[0.01]((1:1,2:1):9,(3:3,4:3):7);\n") == 0 );
}
};
diff --git a/tests/unittests/test_time_interval.cc b/tests/unittests/test_time_interval.cc
index 016ec19..a488c9b 100644
--- a/tests/unittests/test_time_interval.cc
+++ b/tests/unittests/test_time_interval.cc
@@ -23,27 +23,22 @@ class TestTimeInterval : public CppUnit::TestCase {
CPPUNIT_TEST_SUITE_END();
private:
- Forest *forest, *forest2;
- Model *model, *model2;
+ Forest *forest;
+ Model *model;
MersenneTwister *rg;
public:
void setUp() {
rg = new MersenneTwister(5);
model = new Model(5);
- model2 = new Model(5);
- model2->set_population_number(3);
forest = new Forest(model, rg);
forest->createExampleTree();
-
- forest2 = new Forest(model2, rg);
- forest2->createExampleTree();
}
void tearDown() {
- delete forest, forest2;
- delete model, model2;
+ delete forest;
+ delete model;
delete rg;
}
@@ -52,27 +47,27 @@ class TestTimeInterval : public CppUnit::TestCase {
CPPUNIT_ASSERT( (*it).start_height() == 0 );
CPPUNIT_ASSERT( (*it).end_height() == 1 );
CPPUNIT_ASSERT_EQUAL( (size_t)4, it.contemporaries()->size(0) );
-
+
TimeIntervalIterator it2(forest, forest->nodes()->at(4));
CPPUNIT_ASSERT( (*it2).start_height() == 1 );
CPPUNIT_ASSERT( (*it2).end_height() == 3 );
CPPUNIT_ASSERT_EQUAL( (size_t)3, it2.contemporaries()->size(0) );
-
+
TimeIntervalIterator it3(forest, forest->nodes()->at(5));
CPPUNIT_ASSERT( (*it3).start_height() == 3 );
CPPUNIT_ASSERT( (*it3).end_height() == 4 );
CPPUNIT_ASSERT( it3.contemporaries()->size(0) == 2 );
-
+
TimeIntervalIterator it4(forest, forest->nodes()->at(6));
CPPUNIT_ASSERT( (*it4).start_height() == 4 );
CPPUNIT_ASSERT( (*it4).end_height() == 6 );
CPPUNIT_ASSERT( it4.contemporaries()->size(0) == 3 );
-
+
TimeIntervalIterator it5(forest, forest->nodes()->at(7));
CPPUNIT_ASSERT( (*it5).start_height() == 6 );
CPPUNIT_ASSERT( (*it5).end_height() == 10 );
CPPUNIT_ASSERT( it5.contemporaries()->size(0) == 2 );
-
+
TimeIntervalIterator it6(forest, forest->nodes()->at(8));
CPPUNIT_ASSERT( (*it6).start_height() == 10 );
CPPUNIT_ASSERT( (*it6).end_height() == DBL_MAX );
@@ -85,43 +80,43 @@ class TestTimeInterval : public CppUnit::TestCase {
CPPUNIT_ASSERT( (*it).end_height() == 1 );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 4);
CPPUNIT_ASSERT( it.good() );
-
+
++it;
CPPUNIT_ASSERT( (*it).start_height() == 1 );
CPPUNIT_ASSERT( (*it).end_height() == 3 );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 3);
CPPUNIT_ASSERT( it.good() );
-
+
++it;
CPPUNIT_ASSERT( (*it).start_height() == 3 );
CPPUNIT_ASSERT( (*it).end_height() == 4 );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 2 );
CPPUNIT_ASSERT( it.good() );
-
+
++it;
CPPUNIT_ASSERT( (*it).start_height() == 4 );
CPPUNIT_ASSERT( (*it).end_height() == 6 );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 3 );
CPPUNIT_ASSERT( it.good() );
-
+
++it;
CPPUNIT_ASSERT( (*it).start_height() == 6 );
CPPUNIT_ASSERT( (*it).end_height() == 10 );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 2 );
CPPUNIT_ASSERT( it.good() );
-
+
++it;
CPPUNIT_ASSERT( (*it).start_height() == 10 );
CPPUNIT_ASSERT( (*it).end_height() == DBL_MAX );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 0 );
CPPUNIT_ASSERT( it.good() );
-
+
++it;
CPPUNIT_ASSERT( !it.good() );
int i = 0;
for (TimeIntervalIterator events(forest, forest->nodes()->at(0));
- events.good(); ++events) {
+ events.good(); ++events) {
++i;
}
CPPUNIT_ASSERT( i == 6 );
@@ -163,7 +158,7 @@ class TestTimeInterval : public CppUnit::TestCase {
CPPUNIT_ASSERT( (*it).start_height() == 1 );
CPPUNIT_ASSERT( (*it).end_height() == 1.5 );
CPPUNIT_ASSERT( it.contemporaries()->size(0) == 3 );
-
+
CPPUNIT_ASSERT_NO_THROW( it.next() );
CPPUNIT_ASSERT( (*it).start_height() == 1.5 );
CPPUNIT_ASSERT( (*it).end_height() == 3 );
@@ -223,10 +218,10 @@ class TestTimeInterval : public CppUnit::TestCase {
while ( (*it).end_height() < 1000 ) ++it;
forest->nodes()->add(forest->nodes()->createNode(1000));
it.recalculateInterval();
- CPPUNIT_ASSERT( (*it).end_height() == 1000 );
+ CPPUNIT_ASSERT( (*it).end_height() == 1000 );
CPPUNIT_ASSERT_NO_THROW( ++it );
- CPPUNIT_ASSERT( (*it).start_height() == 1000 );
- CPPUNIT_ASSERT( (*it).end_height() == DBL_MAX );
+ CPPUNIT_ASSERT( (*it).start_height() == 1000 );
+ CPPUNIT_ASSERT( (*it).end_height() == DBL_MAX );
}
void testFindContemporaries() {
@@ -234,29 +229,29 @@ class TestTimeInterval : public CppUnit::TestCase {
tii.contemporaries()->clear();
tii.searchContemporariesBottomUp(forest->nodes()->at(4));
- CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
tii.contemporaries()->buffer(forest->nodes()->at(4)->height());
tii.searchContemporariesBottomUp(forest->nodes()->at(5));
- CPPUNIT_ASSERT_EQUAL((size_t)1, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)1, tii.contemporaries()->size(0));
tii.searchContemporariesBottomUp(forest->nodes()->at(5), true);
- CPPUNIT_ASSERT_EQUAL((size_t)1, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)1, tii.contemporaries()->size(0));
tii.contemporaries()->buffer(forest->nodes()->at(4)->height());
tii.searchContemporariesBottomUp(forest->nodes()->at(6));
- CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
tii.searchContemporariesBottomUp(forest->nodes()->at(6), true);
- CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
tii.searchContemporariesBottomUp(forest->nodes()->at(7));
- CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
tii.searchContemporariesBottomUp(forest->nodes()->at(7), true);
- CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)2, tii.contemporaries()->size(0));
tii.searchContemporariesBottomUp(forest->nodes()->at(8));
- CPPUNIT_ASSERT_EQUAL((size_t)0, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)0, tii.contemporaries()->size(0));
tii.searchContemporariesBottomUp(forest->nodes()->at(8), true);
- CPPUNIT_ASSERT_EQUAL((size_t)0, tii.contemporaries()->size(0));
+ CPPUNIT_ASSERT_EQUAL((size_t)0, tii.contemporaries()->size(0));
}
};
diff --git a/tools/compare_versions.sh b/tools/compare_versions.sh
index 7a717b4..a6954bf 100755
--- a/tools/compare_versions.sh
+++ b/tools/compare_versions.sh
@@ -59,10 +59,12 @@ echo "Command | Equal results | time $1 | time $2 | memory $1 | memory $2"
echo "------- | ------------- | ----------- | -------- | ------------- | ---------"
test_scrm 100 100 -t 5 -T
test_scrm 20 100 -r 10 100 -t 3 -T -L
-test_scrm 4 10000 -r 10 100 -t 5
+test_scrm 4 1000 -r 10 100 -t 5 --transpose-segsites
test_scrm 20 2 -r 500 100 -L
-test_scrm 2000 1 -r 10 100 -l 50 -L
+test_scrm 2000 1 -r 10 100 -l 50 -L -T
test_scrm 20 1 -r 4000 10000000 -l 300000 -T
+test_scrm 500 1 -r 1000 10000000 -l 0 -T
+test_scrm 100 1 -r 400 10000000 -l 0 -O
test_scrm 10 100 -r 10 100 -I 3 3 3 4 0.5 -eN 0.1 0.05 -eN 0.2 0.5 -L
test_scrm 50 100 -r 20 200 -G 1.5 -l 100 -L
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/scrm.git
More information about the debian-med-commit
mailing list