Merge pull request #15 from rscherrer/develop

Develop

README.md

@@ -86,6 +86,7 @@ General simulation parameters:
 * `tburnin` (0) is the duration of the burn-in period, in generations
 * `tend` (10) is the number of generations to simulate (after burn-in)
 * `tsave` (10) is the frequency at which to record the data
+* `tcomplete` (1000000) is the time at which to force complete reproductive isolation between the two ecotypes (can mimic e.g. genomic incompatibilities between the two species, or the evolution of very good species recognition abilities)
 * `talkative` (1) is either 0 or 1 and sets whether the simulation should print status information to the prompt
 * `choosewhattosave` (0) is either 0 or 1 and sets whether the variables to save are specified in a separate file, the order file (see below). If 0 all of the output variables are saved every `tsave` generations except for whole genomes
 * `datsave` (1) sets whether to save the recorded variables to files

src/MetaPop.cpp

@@ -81,6 +81,7 @@ void MetaPop::disperse(const Param &p)
 
     }
     */
+
     else {
 
         // Use bernoulli if common
@@ -266,6 +267,8 @@ void MetaPop::reproduce(const Param &p, const GenArch &arch)
 
             std::vector<double> probs = fit;
 
+            const size_t ecof = population[f].getEcotype();
+
             // While the season is not over and mating hasn't occured
             while (t) {
 
@@ -274,6 +277,12 @@ void MetaPop::reproduce(const Param &p, const GenArch &arch)
                 // Sample a male and evaluate
                 const size_t idm = getmale(rnd::rng);
                 const size_t m = males[hab][idm];
+
+                // Skip if male and female are incompatible species
+                const size_t ecom = population[m].getEcotype();
+                if (iscomplete && ecom != ecof)
+                    continue;
+
                 const double xm = population[m].getTraitValue(0u);
                 auto ismating = rnd::bernoulli(population[f].mate(xm, p));
 
@@ -330,6 +339,11 @@ void MetaPop::exitburnin()
     isburnin = false;
 }
 
+void MetaPop::complete()
+{
+    iscomplete = true;
+}
+
 // Getters
 //--------
 

src/MetaPop.h

@@ -23,6 +23,7 @@ class MetaPop
     MetaPop(const Param &pars, const GenArch &arch) :
         population(populate(pars, arch)),
         isburnin(pars.tburnin > 0),
+        iscomplete(false),
         resources(utl::zeros(2u, 2u)),
         sexcounts(utl::uzeros(2u, 2u))
     {
@@ -33,6 +34,7 @@ class MetaPop
 
     void cycle(const Param&, const GenArch&);
     void exitburnin();
+    void complete(); // complete speciation (no more reproduction allowed)
     bool isextinct() const;
 
     void disperse(const Param&);
@@ -68,6 +70,7 @@ class MetaPop
 
     std::vector<Individual> population;
     bool isburnin;
+    bool iscomplete;
 
     std::vector<std::vector<double> > resources; // per habitat per resource
     std::vector<std::vector<size_t> > sexcounts; // per habitat per sex

src/Param.cpp

@@ -39,6 +39,7 @@ Param::Param() :
     tburnin(0),
     tend(10),
     tsave(10),
+    tcomplete(1000000),
     talkative(true),
     datsave(true),
     burninsave(false),
@@ -138,6 +139,7 @@ void Param::import(std::ifstream &file)
         else if (input == "tburnin") file >> tburnin;
         else if (input == "tend") file >> tend;
         else if (input == "tsave") file >> tsave;
+        else if (input == "tcomplete") file >> tcomplete;
         else if (input == "talkative") file >> talkative;
         else if (input == "datsave") file >> datsave;
         else if (input == "burninsave") file >> burninsave;
@@ -332,6 +334,7 @@ void Param::write(std::ofstream &file) const
     file << "tburnin " << tburnin << '\n';
     file << "tend " << tend << '\n';
     file << "tsave " << tsave << '\n';
+    file << "tcomplete " << tcomplete << '\n';
     file << "talkative " << talkative << '\n';
     file << "datsave " << datsave << '\n';
     file << "burninsave " << burninsave << '\n';

src/Param.h

@@ -70,6 +70,7 @@ struct Param {
     int tburnin;
     int tend;
     int tsave;
+    int tcomplete;
     bool talkative;
     bool datsave;
     bool burninsave;

src/Simul.cpp

@@ -53,8 +53,12 @@ int simulate(const std::vector<std::string> &args)
         // Loop through time
         for (int t = -pars.tburnin; t <= pars.tend; ++t) {
 
+            // Exit the burnin if needed
             if (t == 0) metapop.exitburnin();
 
+            // Force-complete speciation if needed
+            if (t == pars.tcomplete) metapop.complete();
+
             if (pars.talkative) std::cout << t << '\n';
 
             // Life cycle of the metapopulation

tests/MetaPopTests.cpp

@@ -1,4 +1,5 @@
 #include "src/MetaPop.h"
+#include "src/Collector.h"
 #include "src/Utilities.h"
 #include <boost/test/unit_test.hpp>