Adding gen int functions for sir and seir conn models

.vscode/launch.json

@@ -11,7 +11,13 @@
             "environment": [],
             "externalConsole": false,
             "MIMode": "gdb",
-            "preLaunchTask": "C/C++: g++ build active file"
+            "preLaunchTask": "C/C++: g++ build active file for debug",
+            "setupCommands": [
+                {
+                    "description": "Enable pretty-printing for gdb",
+                    "text": "-enable-pretty-printing",
+                    "ignoreFailures": true
+                }
         },
     ]
 }

.vscode/tasks.json

@@ -23,6 +23,31 @@
                 "isDefault": true
             },
             "detail": "Task generated by Debugger."
+        },
+        {
+            "type": "cppbuild",
+            "label": "C/C++: g++ build active file for debug",
+            "command": "/usr/bin/g++",
+            "args": [
+                "-fdiagnostics-color=always",
+                "-O0",
+                "-g",
+                "${file}",
+                "-fopenmp",
+                "-o",
+                "${fileDirname}/${fileBasenameNoExtension}"
+            ],
+            "options": {
+                "cwd": "${fileDirname}"
+            },
+            "problemMatcher": [
+                "$gcc"
+            ],
+            "group": {
+                "kind": "build",
+                "isDefault": true
+            },
+            "detail": "Task generated by Debugger."
         }
     ],
     "version": "2.0.0"

examples/09-sir-connected/09-sir-connected.cpp

@@ -54,19 +54,20 @@ int main(int argc, char* argv[]) {
     std::vector< std::vector< int > > results(nexperiments);
     std::vector< std::vector< int > > date(nexperiments);
     std::vector< std::string > labels;
-    epiworld_fast_uint nreplica = 0u;
+    int nreplica = -1;
 
     auto record =
         [&results,&date,&nreplica,&labels](size_t s, epiworld::Model<> * m)
         {
 
+            #pragma omp critical
+            nreplica++;
+
             if (nreplica == 0)
                 m->get_db().get_hist_total(&date[nreplica], &labels, &results[nreplica]);
             else
                 m->get_db().get_hist_total(&date[nreplica], nullptr, &results[nreplica]);
 
-            nreplica++;
-
             return;
 
         };
@@ -78,7 +79,7 @@ int main(int argc, char* argv[]) {
         result["seed"].as<int>(),
         record,       // Function to call after each experiment
         true,         // Whether to reset the population
-        true         // Whether to print a progress bar
+        true          // Whether to print a progress bar
         #ifdef _OPENMP
         ,threads 
         #endif
@@ -106,6 +107,48 @@ int main(int argc, char* argv[]) {
         "total_hist.txt", "transmission.txt", "transition.txt", "", ""
         );
 
+
+    // We can compute the expected generation time
+    auto gen_time = model.generation_time_expected();
+
+    // Compute observed generation interval
+    std::vector< int > agent_id;
+    std::vector< int > virus_id;
+    std::vector< int > time;
+    std::vector< int > gentim;
+    model.get_db().generation_time(
+        agent_id,
+        virus_id,
+        time,
+        gentim
+    );
+
+    // Averaging out at the time level
+    std::vector< double > gen_time_observed;
+    std::vector< double > gen_time_observed_count;
+
+    for (size_t i = 0; i < gentim.size(); i++)
+    {
+        if (gentim[i] >= 0)
+        {
+            if (time[i] >= gen_time_observed.size())
+            {
+                gen_time_observed.resize(time[i] + 1, 0.0);
+                gen_time_observed_count.resize(time[i] + 1, 0.0);
+            }
+
+            gen_time_observed[time[i]] += gentim[i];
+            gen_time_observed_count[time[i]] += 1.0;
+        }
+    }
+
+    for (size_t i = 0; i < gen_time_observed.size(); i++)
+    {
+        if (gen_time_observed_count[i] > 0)
+            gen_time_observed[i] /= gen_time_observed_count[i];
+    }
+
+
     return 0;
 
 }

examples/13-genint/main.cpp

@@ -6,115 +6,6 @@
 
 using namespace epiworld;
 
-// Implementing the factorial function
-inline double log_factorial(int n)
-{
-    if (n == 0)
-        return 0.0;
-    return std::log(static_cast<double>(n)) + log_factorial(n-1);
-}
-
-inline double dpois(int k, double lambda, int max_n = 100, bool as_log = false)
-{
-
-    if (max_n < k)
-        throw std::runtime_error("max_n must be greater than k");
-
-    double res = k * std::log(lambda) - lambda - log_factorial(
-        std::min(k, max_n)
-        );
-
-    return as_log ? res : std::exp(res);
-}
-
-/**
- * @brief Compute the probability of the generation interval
- * 
- * @details
- * If `p_0_approx` is negative, it will be computed using the Poisson
- * distribution.
- * 
- * @param g Generation interval
- * @param S Population size
- * @param p_c Probability of contact
- * @param p_i Probability of infection
- * @param p_r Probability of recovery
- * @param p_0_approx Approximation of the probability of not being infected
- * @param max_n Maximum number of contacts
- * 
- * @return The probability of the generation interval
- * 
- */
-double dgenint(
-    int g,
-    double S,
-    double p_c,
-    double p_i,
-    double p_r,
-    double & p_0_approx,
-    int max_n = 500
-    ) {
-
-    if (p_0_approx < 0.0)
-    {
-
-        p_0_approx = 0.0;
-        for (int i = 0; i < max_n; ++i)
-        {
-
-            p_0_approx +=
-                dpois(i, S * p_c, max_n, false) *
-                std::pow(1.0 - p_i, static_cast<double>(i)) ;
-
-        }
-    }
-
-    double g_dbl = static_cast<double>(g);
-
-    return
-        std::pow(1 - p_r, g_dbl) *
-        std::pow(p_0_approx, g_dbl - 1.0) * 
-        (1.0 - p_0_approx);
-
-}
-
-// Mean of the generation interval
-/**
- * @brief Compute the mean of the generation interval
- * @param S Population size
- * @param p_c Probability of contact
- * @param p_i Probability of infection
- * @param p_r Probability of recovery
- * @param max_n Maximum number of contacts
- * 
- * @return The mean of the generation interval
- */
-double gen_int_mean(
-    double S,
-    double p_c,
-    double p_i,
-    double p_r,
-    int max_n = 200
-    ) {
-
-    double mean = 0.0;
-    double p_0_approx = -1.0;
-    for (int i = 0; i < max_n; ++i)
-    {
-        mean += 
-            static_cast<double>(i) *
-            dgenint(
-                static_cast<double>(i), S, p_c, p_i, p_r, p_0_approx, max_n
-                );
-
-        if (i == 0)
-            std::cout << "p_0_approx: " << p_0_approx << std::endl;
-    }
-
-    return mean;
-
-}
-
 int main()
 {
 
@@ -234,7 +125,7 @@ int main()
     std::cout << "Mean2 (from sim) : " << mean2 << std::endl;
     std::cout << "Mean  (SEIRCONN) : " << mean_seirconn << std::endl;
 
-    double theory_mean = gen_int_mean(popsize, p_contact, p_i, p_r);
+    double theory_mean = epiworld::gen_int_mean(popsize, p_contact, p_i, p_r);
     std::cout << "Mean  (theory) +1: " << theory_mean + 1.0 << std::endl;
 
 

include/epiworld/database-bones.hpp

@@ -246,10 +246,13 @@ class DataBase {
         std::string fn_generation_time
         ) const;
 
+    /***
+     * @brief Record a transmission event
+     */
     void record_transmission(int i, int j, int virus, int i_expo_date);
 
-    size_t get_n_viruses() const;
-    size_t get_n_tools() const;
+    size_t get_n_viruses() const; ///< Get the number of viruses
+    size_t get_n_tools() const; ///< Get the number of tools
 
     void set_user_data(std::vector< std::string > names);
     void add_user_data(std::vector< epiworld_double > x);
@@ -288,19 +291,23 @@ class DataBase {
 
     /**
      * Calculates the generating time
-     * @param agent_id,virus_id,time,gentime vectors where to save the values agent_id
+     * @param agent_id,virus_id,time,gentime vectors where to save the values 
+     * 
+     * @details
+     * The generation time is the time between the infection of the source and 
+     * the infection of the target.
     */
    ///@{
     void generation_time(
         std::vector< int > & agent_id,
         std::vector< int > & virus_id,
         std::vector< int > & time,
         std::vector< int > & gentime
-    ) const;
+    ) const; ///< Get the generation time
 
     void generation_time(
         std::string fn
-    ) const;
+    ) const; ///< Write the generation time to a file
     ///@}
 
 };

include/epiworld/epiworld.hpp

@@ -19,7 +19,7 @@
 /* Versioning */
 #define EPIWORLD_VERSION_MAJOR 0
 #define EPIWORLD_VERSION_MINOR 3
-#define EPIWORLD_VERSION_PATCH 2
+#define EPIWORLD_VERSION_PATCH 3
 
 static const int epiworld_version_major = EPIWORLD_VERSION_MAJOR;
 static const int epiworld_version_minor = EPIWORLD_VERSION_MINOR;
@@ -33,7 +33,7 @@ namespace epiworld {
     #include "misc.hpp"
     #include "progress.hpp"
 
-    // #include "math/summary-stats.hpp"
+    #include "math/distributions.hpp"
 
     #include "math/lfmcmc.hpp"