[med-svn] [scrm] 01/03: Imported Upstream version 1.7.1

Kevin Murray daube-guest at moszumanska.debian.org
Thu Mar 24 00:37:04 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 cc1373121043f0a7c3acda324febbfc329e3407e
Author: Kevin Murray <spam at kdmurray.id.au>
Date:   Thu Mar 24 11:12:07 2016 +1100

    Imported Upstream version 1.7.1
---
 NEWS.md                                  |  14 +++
 configure.ac                             |   2 +-
 src/forest-debug.cc                      |  10 +--
 src/forest.cc                            |  19 ++--
 src/forest.h                             |   4 +-
 src/model.cc                             |  61 +++++++------
 src/model.h                              |  13 +--
 src/summary_statistics/newick_tree.h     |  19 ++--
 src/summary_statistics/oriented_forest.h |   3 +
 src/summary_statistics/seg_sites.cc      |   2 +-
 tests/test_binaries.sh                   |   1 +
 tests/unittests/test_forest.cc           |  19 +++-
 tests/unittests/test_model.cc            |  44 +++++-----
 tests/unittests/test_param.cc            | 146 +++++++++++++++----------------
 14 files changed, 194 insertions(+), 163 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index f655922..73781d2 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,6 +1,20 @@
 scrm Version History
 ========================
 
+scrm 1.7.1
+------------------------
+Released: 2016-03-23
+
+### Bug Fixes
++ Fixes that bug that could lead to false results when the "-eI" option
+  was used to include non-contemporary samples. In that case, scrm
+  could skip the time interval between the MRCA of the contemporary samples
+  and the time of the "-eI" event (if any), and neglect all events 
+  during that time. We recommend to repeat all simulations which
+  used "-eI" with the updated version (#102).
+
+
+
 scrm 1.7.0
 ------------------------
 Released: 2016-02-07
diff --git a/configure.ac b/configure.ac
index 73009df..f348720 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1,4 +1,4 @@
-AC_INIT([scrm], [1.7.0],[https://github.com/paulstaab/scrm/issues])
+AC_INIT([scrm], [1.7.1],[https://github.com/paulstaab/scrm/issues])
 AM_INIT_AUTOMAKE([subdir-objects -Wall -Werror foreign])
 
 # Use -O3 as default optimization level
diff --git a/src/forest-debug.cc b/src/forest-debug.cc
index ab83be2..0bf25d9 100644
--- a/src/forest-debug.cc
+++ b/src/forest-debug.cc
@@ -91,11 +91,11 @@ void Forest::createExampleTree() {
 void Forest::createScaledExampleTree() {
   this->createExampleTree();
 
-  this->nodes()->at(4)->set_height(1 * 4 * model().default_pop_size); 
-  this->nodes()->at(5)->set_height(3 * 4 * model().default_pop_size); 
-  this->nodes()->at(6)->set_height(4 * 4 * model().default_pop_size); 
-  this->nodes()->at(7)->set_height(6 * 4 * model().default_pop_size); 
-  this->nodes()->at(8)->set_height(10 * 4 * model().default_pop_size); 
+  this->nodes()->at(4)->set_height(1 * 4 * model().default_pop_size()); 
+  this->nodes()->at(5)->set_height(3 * 4 * model().default_pop_size()); 
+  this->nodes()->at(6)->set_height(4 * 4 * model().default_pop_size()); 
+  this->nodes()->at(7)->set_height(6 * 4 * model().default_pop_size()); 
+  this->nodes()->at(8)->set_height(10 * 4 * model().default_pop_size()); 
 
   updateAbove(nodes()->at(4));
   updateAbove(nodes()->at(5));
diff --git a/src/forest.cc b/src/forest.cc
index efaaef1..3f02869 100644
--- a/src/forest.cc
+++ b/src/forest.cc
@@ -474,20 +474,21 @@ void Forest::sampleCoalescences(Node *start_node) {
   assert( start_node->is_root() );
   // We can have one or active local nodes: If the coalescing node passes the
   // local root, it also starts a coalescence.
-  set_active_node(0, start_node);
-  set_active_node(1, this->local_root());
+  if (start_node->height() > this->local_root()->height()) {
+    set_active_node(0, this->local_root());
+    set_active_node(1, start_node);
+  } else {
+    set_active_node(0, start_node);
+    set_active_node(1, this->local_root());
+  }
+  assert ( active_node(0)->height() <= active_node(1)->height() );
 
   // Initialize Temporary Variables
-  tmp_event_ = Event(start_node->height());
+  tmp_event_ = Event(active_node(0)->height());
   coalescence_finished_ = false;
 
-  // This assertion needs an exception for building the initial tree
-  assert ( current_rec() == 1 ||
-           active_node(1)->in_sample() ||
-           start_node->height() <= active_node(1)->height() );
-
   // Only prune every second round
-  for (TimeIntervalIterator ti(this, start_node); ti.good(); ++ti) {
+  for (TimeIntervalIterator ti(this, active_node(0)); ti.good(); ++ti) {
 
     dout << "* * Time interval: " << (*ti).start_height() << " - "
          << (*ti).end_height() << " (Last event at " << tmp_event_.time() << ")" << std::endl;
diff --git a/src/forest.h b/src/forest.h
index db2a867..0cd57ca 100644
--- a/src/forest.h
+++ b/src/forest.h
@@ -165,12 +165,12 @@ class Forest
   NodeContainer *nodes() { return &(this->nodes_); }
 
   double getTMRCA(const bool &scaled = false) const {
-    if (scaled) return local_root()->height() / (4 * this->model_->default_pop_size);
+    if (scaled) return local_root()->height() / (4 * this->model_->default_pop_size());
     else return local_root()->height();
   }
 
   double getLocalTreeLength(const bool &scaled = false) const {
-    if (scaled) return local_root()->length_below() / (4 * this->model_->default_pop_size);
+    if (scaled) return local_root()->length_below() / (4 * this->model_->default_pop_size());
     else return local_root()->length_below();
   }
 
diff --git a/src/model.cc b/src/model.cc
index aa3127c..bd08d31 100644
--- a/src/model.cc
+++ b/src/model.cc
@@ -22,14 +22,9 @@
 #include "model.h"
 
 
-Model::Model() {
-  default_pop_size = 10000;
-  default_growth_rate = 0.0;
-  default_mig_rate = 0.0;
-  scaling_factor_ = 1.0 / (4 * default_pop_size);
-
-  has_migration_ = false;
-  has_recombination_ = false;
+Model::Model() : 
+  has_migration_(false),
+  has_recombination_(false) {
 
   this->set_loci_number(1);
   this->setLocusLength(1);
@@ -51,24 +46,32 @@ Model::Model() {
 }
 
 
-Model::Model(size_t sample_size) : Model() {
+Model::Model(size_t sample_size) : 
+  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->setMutationRate(0.0);
+  this->setRecombinationRate(0.0);
+
+  this->window_length_seq_ = 0;
+  this->set_window_length_rec(500);
+
+  this->setSequenceScaling(ms);
+
   this->addSampleSizes(0.0, std::vector<size_t>(1, sample_size));
   this->setLocusLength(1000);
   this->resetTime();
+  this->resetSequencePosition();
 }
 
 
-/*
-void Model::reset() {
-  pop_sizes_list_.clear();
-  growth_rates_list_.clear();
-  mig_rates_list_.clear();
-  total_mig_rates_list_.clear();
-  single_mig_probs_list_.clear();
-  summary_statistics_.clear();
-}
-*/
-
 /**
  * Function to add a new change time to the model.
  *
@@ -83,7 +86,7 @@ void Model::reset() {
  * @returns The position the time has now in the vector
  */
 size_t Model::addChangeTime(double time, const bool &scaled) {
-  if (scaled) time *= 4 * default_pop_size;
+  if (scaled) time *= 4 * default_pop_size();
 
   size_t position = 0;
   if ( change_times_.size() == 0 ) {
@@ -159,7 +162,7 @@ size_t Model::addChangePosition(const double position) {
 
 
 void Model::addSampleSizes(double time, const std::vector<size_t> &samples_sizes, const bool &scaled) {
-  if (scaled) time *= 4 * default_pop_size;
+  if (scaled) time *= 4 * default_pop_size();
 
   for (size_t pop = 0; pop < samples_sizes.size(); ++pop) {
     for (size_t i = 0; i < samples_sizes.at(pop); ++i) {
@@ -196,7 +199,7 @@ void Model::addPopulationSizes(double time, const std::vector<double> &pop_sizes
   for (double pop_size : pop_sizes) {
     if (!std::isnan(pop_size)) {
       // Scale to absolute values if necessary
-      if (relative) { pop_size *= this->default_pop_size; }
+      if (relative) { pop_size *= this->default_pop_size(); }
 
       // Save inverse double value
       if (pop_size <= 0.0) throw std::invalid_argument("population size <= 0");
@@ -247,7 +250,7 @@ void Model::addPopulationSize(const double time, const size_t pop, double popula
                               const bool &time_scaled, const bool &relative) {
   checkPopulation(pop);
   size_t position = addChangeTime(time, time_scaled);
-  if (relative) population_size *= default_pop_size;
+  if (relative) population_size *= default_pop_size();
 
   if (population_size <= 0.0) throw std::invalid_argument("population size <= 0");
   if (pop_sizes_list_.at(position).empty()) addPopulationSizes(time, nan("value to replace"), time_scaled);
@@ -450,7 +453,7 @@ 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 << "N0 is assumed to be " << model.default_pop_size() << std::endl;
 
   model.resetSequencePosition();
   for (size_t idx = 0; idx < model.countChangePositions(); ++idx) {
@@ -548,12 +551,12 @@ 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) {
       if (std::isnan(pop_sizes_list_.at(0).at(pop)))
-        addPopulationSize(0, pop, default_pop_size);
+        addPopulationSize(0, pop, default_pop_size());
     }
   }
 
@@ -676,7 +679,7 @@ void Model::setRecombinationRate(double rate,
     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");
@@ -700,7 +703,7 @@ void Model::setRecombinationRate(double rate,
  */
 void Model::setMutationRate(double rate, const bool &per_locus, const bool &scaled,
                             const double seq_position) {
-  if (scaled) rate /= 4.0 * default_pop_size;
+  if (scaled) rate /= 4.0 * default_pop_size();
 
   size_t idx = addChangePosition(seq_position);
   if (per_locus) {
diff --git a/src/model.h b/src/model.h
index b358732..0a38cc4 100644
--- a/src/model.h
+++ b/src/model.h
@@ -69,11 +69,14 @@ class Model
 
    
    // Default values;
-   double default_pop_size;
-   double default_growth_rate;
-   double default_mig_rate;
+   constexpr static double default_pop_size_ = 10000.0;
+   constexpr static double default_growth_rate = 0.0;
+   constexpr static double default_mig_rate = 0.0;
+   constexpr static double scaling_factor_ = 1.0 / (4 * default_pop_size_);
 
    // Getters & Setters
+   double default_pop_size() const { return Model::default_pop_size_; };
+
 
    /**
     * @brief Returns the scaling factor for times and many parameters
@@ -157,7 +160,7 @@ class Model
    double inv_double_pop_size(const size_t pop = 0, const double time = -1) const { 
      double pop_size;
      if (current_pop_sizes_ == NULL) 
-       pop_size = 1/(2*default_pop_size);
+       pop_size = 1/(2*default_pop_size());
      else 
        pop_size = current_pop_sizes_->at(pop);   // population size basal to this segment
 
@@ -457,8 +460,6 @@ class Model
                           double default_value = nan("value to replace"));
   void addPopToVectorList(std::vector<std::vector<double> > &vector_list);
 
-   double scaling_factor_; // 1 / (4N0);
-
    // Stores information about samples. Each index represents a sample.
    std::vector<size_t> sample_populations_;
    std::vector<double> sample_times_;
diff --git a/src/summary_statistics/newick_tree.h b/src/summary_statistics/newick_tree.h
index 4d58248..ab32fbb 100644
--- a/src/summary_statistics/newick_tree.h
+++ b/src/summary_statistics/newick_tree.h
@@ -41,27 +41,22 @@ struct NewickBuffer {
   std::string tree;      ///< The subtree itself.
 };
 
+
 class NewickTree : public SummaryStatistic
 {
  public:
-  NewickTree() : NewickTree(6, true) { }
-  NewickTree(size_t precision) : NewickTree(precision, true) { }
-  NewickTree(size_t precision, bool has_recombination) { 
-    precision_ = precision; 
-    has_rec_ = has_recombination;
-  }
-
-  ~NewickTree() {}
+  NewickTree() : precision_(6), has_rec_(true) {}
+  NewickTree(size_t precision) : precision_(precision), has_rec_(true) {}
+  NewickTree(size_t precision, bool has_recombination) : 
+      precision_(precision), 
+      has_rec_(has_recombination) {}
 
   //Virtual methods
   void calculate(const Forest &forest);
   void printSegmentOutput(std::ostream &output) const;
 
   NewickTree* clone() const { return new NewickTree(precision_, has_rec_); };
-
-  void clear() {
-    buffer_.clear();
-  }
+  void clear() { buffer_.clear(); }
 
  private:
   std::string generateTree(Node const* node, const Forest &forest, const bool use_buffer);
diff --git a/src/summary_statistics/oriented_forest.h b/src/summary_statistics/oriented_forest.h
index 67f17d8..c8b7bbf 100644
--- a/src/summary_statistics/oriented_forest.h
+++ b/src/summary_statistics/oriented_forest.h
@@ -43,6 +43,9 @@ class OrientedForest : public SummaryStatistic
   void printSegmentOutput(std::ostream &output) const;
   void clear() { }
 
+  std::vector<int> parents() const { return parents_; }
+  std::vector<double> heights() const { return heights_; }
+
   OrientedForest* clone() const {
     return new OrientedForest(this->parents_.size());
   }
diff --git a/src/summary_statistics/seg_sites.cc b/src/summary_statistics/seg_sites.cc
index 1025b40..656ecff 100644
--- a/src/summary_statistics/seg_sites.cc
+++ b/src/summary_statistics/seg_sites.cc
@@ -32,7 +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));
+    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);
diff --git a/tests/test_binaries.sh b/tests/test_binaries.sh
index 30c5960..807304b 100755
--- a/tests/test_binaries.sh
+++ b/tests/test_binaries.sh
@@ -88,5 +88,6 @@ 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 1 -es 1.0 1 0.5 -ej 1.0 2 1 || exit 1
+ test_scrm 4 1 -t 1.0 -I 2 2 0 -ej 1.0 1 2 -eI 3.5 0 2
 echo ""
 
diff --git a/tests/unittests/test_forest.cc b/tests/unittests/test_forest.cc
index 65e9f24..a3833bc 100644
--- a/tests/unittests/test_forest.cc
+++ b/tests/unittests/test_forest.cc
@@ -1,4 +1,3 @@
-#include <cppunit/TestCase.h>
 #include <cppunit/extensions/HelperMacros.h>
 
 #include "../../src/forest.h"
@@ -19,6 +18,7 @@ class TestForest : public CppUnit::TestCase {
   CPPUNIT_TEST( testSamplePoint );
   CPPUNIT_TEST( testCalcRecombinationRate );
   CPPUNIT_TEST( testCalcRate );
+  CPPUNIT_TEST( testCalcRateWithArachicSamples );
   CPPUNIT_TEST( testNodeIsOld );
   CPPUNIT_TEST( testPrune );
   CPPUNIT_TEST( testSelectFirstTime );
@@ -177,6 +177,23 @@ class TestForest : public CppUnit::TestCase {
     delete node2;
   }
 
+  void testCalcRateWithArachicSamples() {
+    Node* node1 = forest->nodes()->createNode(20, 5);
+    forest->nodes()->add(node1);
+    TimeIntervalIterator tii(forest, node1);
+    double pop_size = 2*forest->model().population_size(0);
+
+    forest->set_active_node(0, node1);
+    forest->set_active_node(1, forest->local_root());
+    forest->states_[0] = 1;
+    forest->states_[1] = 1;
+  
+    CPPUNIT_ASSERT_NO_THROW( forest->calcRates(*tii) );
+    CPPUNIT_ASSERT( areSame(1.0/pop_size, forest->rates_[0]) );
+    CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[1] );
+    CPPUNIT_ASSERT_EQUAL( 0.0, forest->rates_[2] );
+  }
+
   void testGetNodeState() {
     CPPUNIT_ASSERT( forest->getNodeState(forest->getNodes()->get(8), 11) == 1 );
     CPPUNIT_ASSERT( forest->getNodeState(forest->getNodes()->get(6), 11) == 2 );
diff --git a/tests/unittests/test_model.cc b/tests/unittests/test_model.cc
index 3cc7c33..f95a120 100644
--- a/tests/unittests/test_model.cc
+++ b/tests/unittests/test_model.cc
@@ -152,12 +152,12 @@ 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_EQUAL( 1.0 * 4 * model.default_pop_size(), model.getCurrentTime() );
     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_EQUAL( 2.0 * 4 * model.default_pop_size(), model.getCurrentTime() );
     CPPUNIT_ASSERT( areSame(model.population_size(0), 10) );
 
     auto pop_sizes2 = std::vector<double>();
@@ -180,17 +180,17 @@ class TestModel : public CppUnit::TestCase {
     model.increaseTime();
 
     CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() );
-    CPPUNIT_ASSERT_EQUAL( 0.5 * model.default_pop_size, model.population_size(0) );
+    CPPUNIT_ASSERT_EQUAL( 0.5 * model.default_pop_size(), model.population_size(0) );
 
     model.addPopulationSizes(1, std::vector<double>(2, .4), true, true);
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size, model.getCurrentTime() );
-    CPPUNIT_ASSERT_EQUAL( 0.4 * model.default_pop_size, model.population_size(0) );
+    CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size(), model.getCurrentTime() );
+    CPPUNIT_ASSERT_EQUAL( 0.4 * model.default_pop_size(), model.population_size(0) );
 
     model.addPopulationSizes(2, 10, true, true);
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size, model.getCurrentTime() );
-    CPPUNIT_ASSERT( areSame(10.0 * model.default_pop_size, model.population_size(0)) );
+    CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size(), model.getCurrentTime() );
+    CPPUNIT_ASSERT( areSame(10.0 * model.default_pop_size(), model.population_size(0)) );
 
     CPPUNIT_ASSERT_THROW( model.addPopulationSizes(3, 0, true, true), std::invalid_argument);
   }
@@ -221,13 +221,13 @@ class TestModel : public CppUnit::TestCase {
     model.finalize();
     model.resetTime();
 
-    CPPUNIT_ASSERT( areSame( model.default_pop_size,  model.population_size() ) );
-    CPPUNIT_ASSERT( areSame( model.default_pop_size,  model.population_size(0) ) );
-    CPPUNIT_ASSERT( areSame( model.default_pop_size,  model.population_size(1) ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size(),  model.population_size() ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size(),  model.population_size(0) ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size(),  model.population_size(1) ) );
     model.increaseTime();
     CPPUNIT_ASSERT_EQUAL( 1.0, model.getCurrentTime() );
-    CPPUNIT_ASSERT( areSame( model.default_pop_size,  model.population_size(0) ) );
-    CPPUNIT_ASSERT( areSame( model.default_pop_size * std::exp( -1.5 ),  model.population_size(0, 2.0) ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size(),  model.population_size(0) ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size()* std::exp( -1.5 ),  model.population_size(0, 2.0) ) );
   }
 
   void testAddMigRates() {
@@ -359,11 +359,9 @@ class TestModel : public CppUnit::TestCase {
 
     CPPUNIT_ASSERT_EQUAL( (size_t)7, model.sample_size() );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(0) );
-    CPPUNIT_ASSERT_EQUAL( model.default_pop_size, model.population_size(0) );
 
     CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); // 1.0
     CPPUNIT_ASSERT_EQUAL( 1.5, model.growth_rate(0) );
-    CPPUNIT_ASSERT_EQUAL( model.default_pop_size, model.population_size(0) );
 
     CPPUNIT_ASSERT_NO_THROW( model.increaseTime() ); // 2.0
     CPPUNIT_ASSERT_EQUAL( 1.5, model.growth_rate(0) );
@@ -405,10 +403,10 @@ class TestModel : public CppUnit::TestCase {
     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() );
+    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*10),
+    CPPUNIT_ASSERT_EQUAL(10.0/(4*model.default_pop_size()*10),
                          model.mutation_rate() );
 
     // Test setting multiple rates
@@ -460,11 +458,11 @@ class TestModel : public CppUnit::TestCase {
     CPPUNIT_ASSERT( model.has_recombination() );
 
     model.setRecombinationRate(0.001, false, true);
-    CPPUNIT_ASSERT_EQUAL( 0.001/(4*model.default_pop_size), model.recombination_rate() );
+    CPPUNIT_ASSERT_EQUAL( 0.001/(4*model.default_pop_size()), model.recombination_rate() );
     CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() );
 
     model.setRecombinationRate(0.001, true, true);
-    CPPUNIT_ASSERT( areSame(0.00001/(4*model.default_pop_size), model.recombination_rate()) );
+    CPPUNIT_ASSERT( areSame(0.00001/(4*model.default_pop_size()), model.recombination_rate()) );
     CPPUNIT_ASSERT_EQUAL( (size_t)101, model.loci_length() );
 
     // Test setting multiple rates
@@ -551,8 +549,8 @@ class TestModel : public CppUnit::TestCase {
 
 
     model.resetTime();
-    double size0 = model.default_pop_size;
-    double size1 = model.default_pop_size;
+    double size0 = model.default_pop_size();
+    double size1 = model.default_pop_size();
     CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) );
     CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) );
 
@@ -569,12 +567,12 @@ class TestModel : public CppUnit::TestCase {
     model.increaseTime();
     CPPUNIT_ASSERT_EQUAL( 1.5, model.getCurrentTime() );
     CPPUNIT_ASSERT( areSame( 500, model.population_size(0) ) );
-    CPPUNIT_ASSERT( areSame( model.default_pop_size*std::exp(-0.5*0.5), model.population_size(1) ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size()*std::exp(-0.5*0.5), model.population_size(1) ) );
 
     model.increaseTime();
     CPPUNIT_ASSERT_EQUAL( 2.5, model.getCurrentTime() );
     size0 = 500*std::exp(0.5*1.0);
-    size1 = model.default_pop_size*std::exp(-0.5*1.5);
+    size1 = model.default_pop_size()*std::exp(-0.5*1.5);
     CPPUNIT_ASSERT( areSame( size0, model.population_size(0) ) );
     CPPUNIT_ASSERT( areSame( size1, model.population_size(1) ) );
 
@@ -599,7 +597,7 @@ class TestModel : public CppUnit::TestCase {
     model.finalize();
     //std::cout << model << std::endl;
     model.resetTime();
-    CPPUNIT_ASSERT( areSame( model.default_pop_size, model.population_size(0) ) );
+    CPPUNIT_ASSERT( areSame( model.default_pop_size(), model.population_size(0) ) );
     CPPUNIT_ASSERT( areSame( 500.0, model.population_size(1) ) );
   }
 
diff --git a/tests/unittests/test_param.cc b/tests/unittests/test_param.cc
index 5b176df..b2a6682 100644
--- a/tests/unittests/test_param.cc
+++ b/tests/unittests/test_param.cc
@@ -111,8 +111,8 @@ class TestParam : public CppUnit::TestCase {
     Model model = Param("4 7 -t 40.04 -r 1.24 1001 -l 1000").parse();
     CPPUNIT_ASSERT_EQUAL( (size_t)4, model.sample_size() );
     CPPUNIT_ASSERT_EQUAL( (size_t)7, model.loci_number() );
-    CPPUNIT_ASSERT( areSame(40.04/(4*model.default_pop_size*1001), model.mutation_rate()) );
-    CPPUNIT_ASSERT( areSame(1.24/(4*model.default_pop_size*1000), model.recombination_rate()) );
+    CPPUNIT_ASSERT( areSame(40.04/(4*model.default_pop_size()*1001), model.mutation_rate()) );
+    CPPUNIT_ASSERT( areSame(1.24/(4*model.default_pop_size()*1000), model.recombination_rate()) );
     CPPUNIT_ASSERT_EQUAL( (size_t)1001, model.loci_length() );
     CPPUNIT_ASSERT( model.has_approximation() );
     CPPUNIT_ASSERT( model.has_window_seq() );
@@ -141,9 +141,9 @@ class TestParam : public CppUnit::TestCase {
     CPPUNIT_ASSERT_EQUAL( model.sample_time(17), (double)0.0 );
     //std::cout << model << std::endl;
     model.finalize();
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(0, 1)) );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(0, 2)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(0, 2)) );
 
     model = Param("2 1 -t 3.74 -I 2 1 1 5.0").parse();
     CPPUNIT_ASSERT_EQUAL( (size_t)2, model.population_number() );
@@ -154,21 +154,21 @@ class TestParam : public CppUnit::TestCase {
     model = Param("23 10 -t 3.74 -I 3 7 8 5 -eI 12.3 2 0 1 -M 5.0").parse();
     CPPUNIT_ASSERT_EQUAL( model.sample_population(20), (size_t)0 );
     CPPUNIT_ASSERT_EQUAL( model.sample_population(22), (size_t)2 );
-    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) );
+    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
     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() );
-    CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size, model.population_size(0)) );
-    CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size, model.population_size(1)) );
-    CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size, model.population_size(2)) );
+    CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size(), model.population_size(0)) );
+    CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size(), model.population_size(1)) );
+    CPPUNIT_ASSERT( areSame(0.3*model.default_pop_size(), model.population_size(2)) );
     model.increaseTime();
-    CPPUNIT_ASSERT( areSame(8.2 * 4 * model.default_pop_size, model.getCurrentTime()) );
-    CPPUNIT_ASSERT( areSame(0.75*model.default_pop_size, model.population_size(0)) );
-    CPPUNIT_ASSERT( areSame(0.75*model.default_pop_size, model.population_size(1)) );
-    CPPUNIT_ASSERT( areSame(0.75*model.default_pop_size, model.population_size(2)) );
+    CPPUNIT_ASSERT( areSame(8.2 * 4 * model.default_pop_size(), model.getCurrentTime()) );
+    CPPUNIT_ASSERT( areSame(0.75*model.default_pop_size(), model.population_size(0)) );
+    CPPUNIT_ASSERT( areSame(0.75*model.default_pop_size(), model.population_size(1)) );
+    CPPUNIT_ASSERT( areSame(0.75*model.default_pop_size(), model.population_size(2)) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(0) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(1) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(2) );
@@ -178,23 +178,23 @@ class TestParam : public CppUnit::TestCase {
     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, 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)) );
+    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(), 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)) );
     model.increaseTime();
-    CPPUNIT_ASSERT( areSame(1.1 * 4 * model.default_pop_size, model.getCurrentTime()) );
+    CPPUNIT_ASSERT( areSame(1.1 * 4 * model.default_pop_size(), model.getCurrentTime()) );
     model.increaseTime();
-    CPPUNIT_ASSERT( areSame(1.5 * 4 * model.default_pop_size, model.getCurrentTime()) );
+    CPPUNIT_ASSERT( areSame(1.5 * 4 * model.default_pop_size(), model.getCurrentTime()) );
     model.increaseTime();
-    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( 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( 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( 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) );
   }
 
@@ -202,52 +202,52 @@ class TestParam : public CppUnit::TestCase {
     // -ma
     Model model = Param("20 10 -I 2 10 10 -ma x 5 7 x").parse();
     model.resetTime();
-    CPPUNIT_ASSERT( areSame(5.0/(4*model.default_pop_size), model.migration_rate(0, 1)) );
-    CPPUNIT_ASSERT( areSame(7.0/(4*model.default_pop_size), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(5.0/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(7.0/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
 
     // -ema
     model = Param("20 10 -I 2 10 10 -ema 1.6 x 5 7 x").parse();
     model.resetTime();
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size, model.getCurrentTime() );
-    CPPUNIT_ASSERT( areSame(5.0/(4*model.default_pop_size), model.migration_rate(0, 1)) );
-    CPPUNIT_ASSERT( areSame(7.0/(4*model.default_pop_size), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size(), model.getCurrentTime() );
+    CPPUNIT_ASSERT( areSame(5.0/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(7.0/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
 
     // -M
     model = Param("30 10 -I 3 10 10 10 0 -M 5").parse();
     model.resetTime();
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(0, 1)) );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(0, 2)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(0, 2)) );
 
     // -eM
     model = Param("30 10 -I 3 10 10 10 0 -eM 1.6 5").parse();
     model.resetTime();
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size,  model.getCurrentTime() );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(0, 1)) );
-    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size), model.migration_rate(0, 2)) );
+    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size(),  model.getCurrentTime() );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(2.5/(4*model.default_pop_size()), model.migration_rate(0, 2)) );
 
     model = Param("20 1 -I 2 13 7 -m 1 2 1.5").parse();
     //std::cout << model << std::endl;
     CPPUNIT_ASSERT( areSame(0.0, model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
 
     model = Param("20 1 -I 2 13 7 -m 1 2 1.5 -m 2 1 0.5").parse();
     //std::cout << model << std::endl;
-    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size), model.migration_rate(0, 1)) );
-    CPPUNIT_ASSERT( areSame(0.5/(4*model.default_pop_size), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(0.5/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
 
     model = Param("20 1 -I 2 13 7 -em 1.0 2 1 0.5 -em 2.0 2 1 0.6").parse();
-    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size), model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
     model.increaseTime();
-    CPPUNIT_ASSERT( areSame(0.5/(4*model.default_pop_size), model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(0.5/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
     model.increaseTime();
-    CPPUNIT_ASSERT( areSame(0.6/(4*model.default_pop_size), model.migration_rate(1, 0)) );
-    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size), model.migration_rate(0, 1)) );
+    CPPUNIT_ASSERT( areSame(0.6/(4*model.default_pop_size()), model.migration_rate(1, 0)) );
+    CPPUNIT_ASSERT( areSame(0.0/(4*model.default_pop_size()), model.migration_rate(0, 1)) );
   }
 
   void testParseMergeOptions() {
@@ -256,26 +256,24 @@ class TestParam : public CppUnit::TestCase {
     CPPUNIT_ASSERT_NO_THROW(model = Param("20 10 -I 2 10 10 1.5 -es 1.6 2 0.5").parse());
     model.resetTime();
     CPPUNIT_ASSERT_EQUAL( (size_t)3, model.population_number() );
-    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size), model.migration_rate(0,1)) );
-    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size), model.migration_rate(1,0)) );
+    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size()), model.migration_rate(0,1)) );
+    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size()), model.migration_rate(1,0)) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(0,2) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(2,0) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(0,2) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(2,1) );
-    CPPUNIT_ASSERT_EQUAL( model.default_pop_size, model.population_size(2) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(2) );
 
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size, model.getCurrentTime() );
+    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size(), model.getCurrentTime() );
     CPPUNIT_ASSERT_EQUAL( 0.5, model.single_mig_pop(1, 2) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.single_mig_pop(2, 1) );
-    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size), model.migration_rate(0,1)) );
-    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size), model.migration_rate(1,0)) );
+    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size()), model.migration_rate(0,1)) );
+    CPPUNIT_ASSERT( areSame(1.5/(4*model.default_pop_size()), model.migration_rate(1,0)) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(0,2) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(2,0) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(1,2) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(2,1) );
-    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(),
@@ -294,12 +292,12 @@ class TestParam : public CppUnit::TestCase {
     CPPUNIT_ASSERT_NO_THROW(model = Param("20 10 -I 2 10 10 -ej 1.6 2 1 -M 1.3").parse());
     model.resetTime();
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size, model.getCurrentTime() );
+    CPPUNIT_ASSERT_EQUAL( 1.6 * 4 * model.default_pop_size(), model.getCurrentTime() );
     CPPUNIT_ASSERT_EQUAL( 1.0, model.single_mig_pop(1, 0) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.single_mig_pop(0, 1) );
 
     CPPUNIT_ASSERT_EQUAL( 0.0, model.migration_rate(0, 1) );
-    CPPUNIT_ASSERT_EQUAL( 1.3/(4*model.default_pop_size), model.migration_rate(1, 0) );
+    CPPUNIT_ASSERT_EQUAL( 1.3/(4*model.default_pop_size()), model.migration_rate(1, 0) );
   }
 
   void testParseOutputOptions() {
@@ -335,12 +333,12 @@ class TestParam : public CppUnit::TestCase {
       "-G", "2", "-eG", "1", "3", "-M", "5.0" };
     CPPUNIT_ASSERT_NO_THROW( model = Param(16, argv).parse(); );
     model.resetTime();
-    CPPUNIT_ASSERT( areSame(2.0 / 4 / model.default_pop_size, model.growth_rate(0)) );
-    CPPUNIT_ASSERT( areSame(2.0 / 4 / model.default_pop_size, model.growth_rate(1)) );
+    CPPUNIT_ASSERT( areSame(2.0 / 4 / model.default_pop_size(), model.growth_rate(0)) );
+    CPPUNIT_ASSERT( areSame(2.0 / 4 / model.default_pop_size(), model.growth_rate(1)) );
     model.increaseTime();
-    CPPUNIT_ASSERT( areSame(1.0 * 4 * model.default_pop_size, model.getCurrentTime()) );
-    CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size, model.growth_rate(0)) );
-    CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size, model.growth_rate(1)) );
+    CPPUNIT_ASSERT( areSame(1.0 * 4 * model.default_pop_size(), model.getCurrentTime()) );
+    CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size(), model.growth_rate(0)) );
+    CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size(), model.growth_rate(1)) );
 
     // -g && -eg
     char *argv2[] = {
@@ -348,14 +346,14 @@ class TestParam : public CppUnit::TestCase {
       "-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();
-    CPPUNIT_ASSERT_EQUAL( model.default_growth_rate, model.growth_rate(0) );
-    CPPUNIT_ASSERT_EQUAL( 0.1 / 4 / model.default_pop_size, model.growth_rate(1) );
+    CPPUNIT_ASSERT_EQUAL( 0.0, model.growth_rate(0) );
+    CPPUNIT_ASSERT_EQUAL( 0.1 / 4 / model.default_pop_size(), model.growth_rate(1) );
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size, model.getCurrentTime() );
+    CPPUNIT_ASSERT_EQUAL( 1.0 * 4 * model.default_pop_size(), model.getCurrentTime() );
     model.increaseTime();
-    CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size, model.getCurrentTime() );
-    CPPUNIT_ASSERT( areSame(2.4 / 4 / model.default_pop_size, model.growth_rate(0)) );
-    CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size, model.growth_rate(1)) );
+    CPPUNIT_ASSERT_EQUAL( 2.0 * 4 * model.default_pop_size(), model.getCurrentTime() );
+    CPPUNIT_ASSERT( areSame(2.4 / 4 / model.default_pop_size(), model.growth_rate(0)) );
+    CPPUNIT_ASSERT( areSame(3.0 / 4 / model.default_pop_size(), model.growth_rate(1)) );
   }
 
   void testParseVariableRates() {
@@ -363,7 +361,7 @@ class TestParam : public CppUnit::TestCase {
     Model model;
     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 );
+    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( 0.5, model.getNextSequencePosition() );
@@ -381,7 +379,7 @@ class TestParam : public CppUnit::TestCase {
     char *argv2[] = { "scrm", "20", "10", "-r", "3.74", "100", "-sr", "10", "5.1", "-sr", "4.5", "3.2"};
     CPPUNIT_ASSERT_NO_THROW( model = Param(12, argv2).parse(); );
     model.resetSequencePosition();
-    scale = 1 / ( 4 * model.default_pop_size * 99 );
+    scale = 1 / ( 4 * model.default_pop_size() * 99 );
     CPPUNIT_ASSERT( areSame(3.74 * scale, model.recombination_rate()) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.getCurrentSequencePosition() );
     CPPUNIT_ASSERT_EQUAL( 4.5, model.getNextSequencePosition() );
@@ -400,8 +398,8 @@ class TestParam : public CppUnit::TestCase {
     CPPUNIT_ASSERT_NO_THROW( model = Param(12, argv3).parse(); );
     model.resetSequencePosition();
     model.resetTime();
-    scale = 1 / ( 4 * model.default_pop_size * 99 );
-    double scale2 = 1 / ( 4 * model.default_pop_size * 100 );
+    scale = 1 / ( 4 * model.default_pop_size() * 99 );
+    double scale2 = 1 / ( 4 * model.default_pop_size() * 100 );
 
     CPPUNIT_ASSERT( areSame(3.74 * scale, model.recombination_rate()) );
     CPPUNIT_ASSERT_EQUAL( 0.0, model.mutation_rate() );
@@ -467,7 +465,7 @@ class TestParam : public CppUnit::TestCase {
     Param pars;
     pars = Param("4 7 -r 1e3 1001");
     CPPUNIT_ASSERT_NO_THROW(model = pars.parse());
-    CPPUNIT_ASSERT_EQUAL(1.0/(4*model.default_pop_size), model.recombination_rate());
+    CPPUNIT_ASSERT_EQUAL(1.0/(4*model.default_pop_size()), model.recombination_rate());
 
     pars = Param("4 7 -r 0.5 2e3");
     CPPUNIT_ASSERT_NO_THROW(model = pars.parse());

-- 
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