Merge pull request #14 from rscherrer/develop

save data during burnin

README.md

@@ -89,6 +89,7 @@ General simulation parameters:
 * `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
+* `burninsave` (0) sets whether to save data during the burn-in phase too (time points belonging to the burn-in are negative)
 * `gensave` (0) is either 0 or 1 and sets whether whole genomes should be saved every `tfreeze` generations
 * `archsave` (0) is either 0 or 1 and sets whether the genetic architecture should be saved into a file (see below)
 * `archload` (0) sets whether the genetic architecture of the simulation should be loaded from a file instead of generated anew

src/Param.cpp

@@ -41,6 +41,7 @@ Param::Param() :
     tsave(10),
     talkative(true),
     datsave(true),
+    burninsave(false),
     choosewhattosave(false),
     gensave(false),
     archsave(false),
@@ -139,6 +140,7 @@ void Param::import(std::ifstream &file)
         else if (input == "tsave") file >> tsave;
         else if (input == "talkative") file >> talkative;
         else if (input == "datsave") file >> datsave;
+        else if (input == "burninsave") file >> burninsave;
         else if (input == "choosewhattosave") file >> choosewhattosave;
         else if (input == "gensave") file >> gensave;
         else if (input == "archsave") file >> archsave;
@@ -332,6 +334,7 @@ void Param::write(std::ofstream &file) const
     file << "tsave " << tsave << '\n';
     file << "talkative " << talkative << '\n';
     file << "datsave " << datsave << '\n';
+    file << "burninsave " << burninsave << '\n';
     file << "choosewhattosave " << choosewhattosave << '\n';
     file << "gensave " << gensave << '\n';
     file << "archsave " << archsave << '\n';

src/Param.h

@@ -72,6 +72,7 @@ struct Param {
     int tsave;
     bool talkative;
     bool datsave;
+    bool burninsave;
     bool choosewhattosave;
     bool gensave;
     bool archsave;

src/Printer.cpp

@@ -38,13 +38,13 @@ void Printer::shutdown()
         files[f]->close();
 }
 
-void Printer::print(const size_t &t, const Collector &c, const MetaPop &m)
+void Printer::print(const int &t, const Collector &c, const MetaPop &m)
 {
 
     for (size_t f = 0u; f < filenames.size(); ++f) {
 
         if (filenames[f] == "time")
-            stf::write(utl::size2dbl(t), files[f]);
+            stf::write(utl::int2dbl(t), files[f]);
         else if (filenames[f] == "population_size")
             stf::write(utl::size2dbl(c.counts[2u][2u]), files[f]);
         else if (filenames[f] == "ecotype_size")

src/Printer.h

@@ -18,7 +18,7 @@ class Printer
     ~Printer();
 
     void open();
-    void print(const size_t&, const Collector&, const MetaPop&);
+    void print(const int&, const Collector&, const MetaPop&);
     void shutdown();
 
 private:

src/Simul.cpp

@@ -1,10 +1,5 @@
 #include "Simul.h"
 
-bool timetosave(const int &t, const int &tsave)
-{
-    return t >= 0 && t % tsave == 0;
-}
-
 int simulate(const std::vector<std::string> &args)
 {
 
@@ -66,15 +61,16 @@ int simulate(const std::vector<std::string> &args)
             metapop.disperse(pars);
             metapop.consume(pars);
 
+            const bool timetosave = t % pars.tsave;
+
             // Analyze the metapopulation if needed
-            if (pars.datsave && timetosave(t, pars.tsave)) {
+            if (pars.datsave && (t >= 0 || pars.burninsave) && timetosave) {
 
                 // Collect stats
                 collector.analyze(metapop, pars, arch);
 
                 // Save them to files
-                const size_t tu = static_cast<size_t>(t);
-                printer.print(tu, collector, metapop);
+                printer.print(t, collector, metapop);
 
                 // Save whole genomes if needed (space-consuming)
                 if (pars.gensave) freezer.freeze(metapop, pars.nloci);

src/Utilities.cpp

@@ -190,6 +190,12 @@ void utl::correct(double &x, const double &x0, const double &d)
     x = x < x0 + d ? x0 : x;
 }
 
+// Convert integer to double
+double utl::int2dbl(const int &x)
+{
+    return static_cast<double>(x);
+}
+
 // Convert unsigned integer to double
 double utl::size2dbl(const size_t &x)
 {

src/Utilities.h

@@ -29,6 +29,7 @@ namespace utl
     std::vector<std::vector<double> > dividemat(
      const std::vector<std::vector<double> >&,
       const std::vector<std::vector<size_t> >&);
+    double int2dbl(const int&);
     double size2dbl(const size_t&);
     size_t dbl2size(const double&);
     double round(const double&, const size_t&);

tests/OutputTests.cpp

@@ -76,6 +76,44 @@ BOOST_AUTO_TEST_CASE(OutputFilesAreCorrectlyWritten)
 
 }
 
+// Test that burnin is saved correctly (negative time values)
+BOOST_AUTO_TEST_CASE(SaveBurnin)
+{
+
+    Param pars;
+    pars.tburnin = 10;
+    pars.tend = 10;
+
+    GenArch arch = GenArch(pars);
+    MetaPop metapop = MetaPop(pars, arch);
+    Collector collector = Collector(arch);
+    Printer printer = Printer();
+    printer.open();
+
+    // This should save also during the burnin
+    for (int t = -pars.tburnin; t <= pars.tend; ++t) {
+
+        if (t == 0) metapop.exitburnin();
+
+        // Collect stats
+        collector.analyze(metapop, pars, arch);
+
+        // Save them to files
+        printer.print(t, collector, metapop);
+
+    }
+
+    printer.shutdown();
+
+    // Read output files
+    std::vector<double> time = tst::readfile("time.dat");
+
+    // Check that the first time point saved is negative (burnin)
+    BOOST_CHECK(time[0u] < 0.0);
+
+}
+
+
 BOOST_AUTO_TEST_CASE(SaveOneGenome)
 {