Fixed convergence issues in PEB model

embecosm · Nov 11, 2022 · 4e55393 · 4e55393
1 parent e746224
commit 4e55393
Show file tree

Hide file tree

Showing 38 changed files with 1,174 additions and 1,721 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -9,6 +9,8 @@
 
 cmake_minimum_required(VERSION 3.16)
 
+project(dcm_demo LANGUAGES CXX)
+
 SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}-std=c++14 -O3 -flto -march=native\
  -D_GLIBCXX_PARALLEL -DNDEBUG")
 SET(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}-O3")
@@ -17,7 +19,7 @@ SET(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS}-O3")
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
-project(dcm_demo LANGUAGES CXX)
+find_package(OpenMP)
 
 SET(EIGEN)
 add_subdirectory(${CMAKE_SOURCE_DIR}/lib/googletest)
@@ -79,8 +81,6 @@ set(SOURCES_SERIALIZATION_TESTS
     tests/serialization_test.cc
 )
 
-find_package(OpenMP)
-
 add_executable(dcm_covid ${SOURCES_COVID})
 if(OpenMP_FOUND)
     target_link_libraries(dcm_covid PUBLIC OpenMP::OpenMP_CXX)

diff --git a/README.md b/README.md
@@ -11,6 +11,10 @@ This codebase is documented with Doxygen and Sphinx. Online versions of
 documentation can be found [here](https://embecosm.github.io/dcEmb_docs/),
 and also come attached to the repository. 
 
+To build sphinx documentation:
+1) sphinx-build -b html man/source man/build
+2) doxygen doc/Doxyfile.in
+
 # Getting Started
 See the online manual [here](https://embecosm.github.io/dcEmb_docs/).
 

diff --git a/include/peb_model.hh b/include/peb_model.hh
@@ -137,7 +137,9 @@ class peb_model : public dynamic_model {
                               this->GCM.at(i).conditional_parameter_covariances(
                                   this->random_effects, this->random_effects) *
                               singular_vec;
-      c_p_c = (c_p_c.inverse() + (p_p_c.inverse() / 16)).inverse();
+
+      c_p_c = utility::inverse_tol(utility::inverse_tol(c_p_c) +
+                                   (utility::inverse_tol(p_p_c) / 16));
       perDCM_prior_p_e.push_back(p_p_e);
       perDCM_prior_p_c.push_back(p_p_c);
       perDCM_posterior_p_e.push_back(c_p_e);
@@ -153,7 +155,7 @@ class peb_model : public dynamic_model {
 
     Eigen::MatrixXd prior_p_q =
         singular_vec.transpose() * prior_p_c * singular_vec;
-    prior_p_q = prior_p_q.inverse();
+    prior_p_q = utility::inverse_tol(prior_p_q);
 
     Eigen::VectorXd kron_b_e =
         Eigen::MatrixXd::Identity(num_between_effects, num_between_effects)
@@ -169,13 +171,13 @@ class peb_model : public dynamic_model {
         kroneckerProduct(kron_b_e, singular_vec.transpose() * tmp_prior_b_e);
     Eigen::MatrixXd prior_b_c = kroneckerProduct(
         kron_b_c, (singular_vec.transpose() * tmp_prior_b_c * singular_vec));
-    Eigen::MatrixXd prior_b_q = prior_b_c.inverse();
+    Eigen::MatrixXd prior_b_q = utility::inverse_tol(prior_b_c);
 
     Eigen::VectorXd prior_g_e = Eigen::MatrixXd(1, 1);
     prior_g_e << 0;
     Eigen::MatrixXd prior_g_c = Eigen::MatrixXd(1, 1);
     prior_g_c << 1.0 / 16;
-    Eigen::MatrixXd prior_g_q = prior_g_c.inverse();
+    Eigen::MatrixXd prior_g_q = utility::inverse_tol(prior_g_c);
     Eigen::MatrixXd prior_bg_q =
         Eigen::MatrixXd::Zero(prior_b_q.rows() + prior_g_q.rows(),
                               prior_b_q.cols() + prior_g_q.cols());
@@ -205,15 +207,15 @@ class peb_model : public dynamic_model {
     Eigen::MatrixXd current_dFdgg;
     double t = -4;
     double dF = 0;
-    double Fc;
+    double Fc = 0;
     for (int i = 0; i < this->max_invert_it; i++) {
       if (num_precision_comps > 0) {
         prior_r_p_q = prior_p_q * exp(-8);
         for (int j = 0; j < num_precision_comps; j++) {
           prior_r_p_q += exp(g_estimate(j)) * prior_p_q;
         }
       }
-      prior_r_p_c = prior_r_p_q.inverse();
+      prior_r_p_c = utility::inverse_tol(prior_r_p_q);
 
       double free_energy_tmp = 0;
 
@@ -298,8 +300,7 @@ class peb_model : public dynamic_model {
 
       double Fb = b_estimate.transpose() * prior_b_q * b_estimate;
       double Fg = g_estimate.transpose() * prior_g_q * g_estimate;
-      utility::print_matrix("prior_bg_q", prior_bg_q);
-      Fc = Fb / 2 + Fg / 2 + utility::logdet(prior_bg_q * posterior_p_c) / 2;
+      Fc = Fb / 2 + Fg / 2 - utility::logdet(prior_bg_q * posterior_p_c) / 2;
 
       free_energy_tmp = free_energy_tmp - Fc;
 
@@ -330,7 +331,18 @@ class peb_model : public dynamic_model {
       }
 
       Eigen::VectorXd dp = utility::dx(dFdpp, dFdp, t);
-      utility::print_matrix("dp", dp);
+
+      if (dp.norm() >= 8) {
+        dFdpp = Eigen::MatrixXd::Zero(dFdbb.rows() + dFdgg.rows(),
+                                      dFdbb.cols() + dFdgg.cols());
+        dFdpp(Eigen::seq(0, dFdbb.rows() - 1),
+              Eigen::seq(0, dFdbb.cols() - 1)) = dFdbb;
+        dFdpp(Eigen::seq(dFdbb.rows(), dFdbb.rows() + dFdgg.rows() - 1),
+              Eigen::seq(dFdbb.cols(), dFdbb.cols() + dFdgg.cols() - 1)) =
+            dFdgg;
+        dp = utility::dx(dFdpp, dFdp, t);
+      }
+
       Eigen::VectorXd db = dp(Eigen::seq(0, b_estimate.size() - 1));
       Eigen::VectorXd dg = dp(Eigen::seq(
           b_estimate.size(), b_estimate.size() + g_estimate.size() - 1));
@@ -351,6 +363,9 @@ class peb_model : public dynamic_model {
       prior_e_p_c = utility::inverse_tol(prior_e_p_c);
 
       Eigen::VectorXd prior_e_p_e = GCM.at(i).prior_parameter_expectations;
+      //   std::cout << "perDCM_prior_r_p_e[i]" << perDCM_prior_r_p_e[i] <<
+      //   '\n'; std::cout << "GCM.at(i).prior_parameter_expectations"
+      //             << GCM.at(i).prior_parameter_expectations << '\n';
       prior_e_p_e(this->random_effects) = singular_vec * perDCM_prior_r_p_e[i];
 
       Eigen::MatrixXd tmp_prior_p_c =

diff --git a/src/COVID/DEM_COVID.cc b/src/COVID/DEM_COVID.cc
@@ -32,7 +32,7 @@
  * Run the COVID example
  */
 int run_COVID_test() {
-  bool OCTAVE_INVERSIONS = 0;
+  bool OCTAVE_INVERSIONS = 1;
   int num_countries = 5;
   std::vector<dynamic_COVID_model> GCM;
   std::vector<country_data> countries = read_country_data(num_countries);
@@ -125,16 +125,16 @@ int run_COVID_test() {
   PEB.max_invert_it = 64;
   PEB.invert_model();
 
+  std::cout << "PEB.prior_parameter_expectations" << PEB.prior_parameter_expectations << '\n';
+  std::cout << "PEB.posterior_parameter_expectations" << PEB.conditional_parameter_expectations << '\n';
   generate_PEB_values(PEB);
+  std::cout << "PEB.prior_parameter_expectations" << PEB.prior_parameter_expectations << '\n';
+  std::cout << "PEB.posterior_parameter_expectations" << PEB.conditional_parameter_expectations << '\n';
 
   bmr_model<peb_model<dynamic_COVID_model>> BMR;
   BMR.DCM_in = PEB;
   BMR.reduce();
   bma_model<peb_model<dynamic_COVID_model>> BMA = BMR.BMA;
-  utility::print_matrix("BMR.free_energy_vector", BMR.free_energy_vector);
-  utility::print_matrix("BMR.posterior_probability_vector",
-                        BMR.posterior_probability_vector);
-  utility::print_matrix("BMR.model_space", BMR.model_space);
 
   std::vector<dynamic_COVID_model> GCM_empirical = PEB.empirical_GCM;
   for (int i = 0; i < GCM_empirical.size(); i++) {
@@ -341,10 +341,10 @@ std::vector<dynamic_COVID_model> generate_comparison_data(
 void generate_us_posteriors(dynamic_COVID_model& model) {
   model.conditional_parameter_expectations =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/US_conditional_parameter_expectations.csv");
+          "../src/data/US_posterior_p_e.csv");
   model.conditional_parameter_covariances =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/US_conditional_parameter_covariances.csv");
+          "../src/data/US_posterior_p_c.csv");
   model.free_energy = -11263.38592908148;
 }
 
@@ -354,10 +354,10 @@ void generate_us_posteriors(dynamic_COVID_model& model) {
 void generate_brazil_posteriors(dynamic_COVID_model& model) {
   model.conditional_parameter_expectations =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/Brazil_conditional_parameter_expectations.csv");
+          "../src/data/Brazil_posterior_p_e.csv");
   model.conditional_parameter_covariances =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/Brazil_conditional_parameter_covariances.csv");
+          "../src/data/Brazil_posterior_p_c.csv");
   model.free_energy = -9205.455950138792;
 }
 
@@ -367,10 +367,10 @@ void generate_brazil_posteriors(dynamic_COVID_model& model) {
 void generate_india_posteriors(dynamic_COVID_model& model) {
   model.conditional_parameter_expectations =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/India_conditional_parameter_expectations.csv");
+          "../src/data/India_posterior_p_e.csv");
   model.conditional_parameter_covariances =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/India_conditional_parameter_covariances.csv");
+          "../src/data/India_posterior_p_c.csv");
   model.free_energy = -10481.23075593396;
 }
 
@@ -380,10 +380,10 @@ void generate_india_posteriors(dynamic_COVID_model& model) {
 void generate_russia_posteriors(dynamic_COVID_model& model) {
   model.conditional_parameter_expectations =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/Russia_conditional_parameter_expectations.csv");
+          "../src/data/Russia_posterior_p_e.csv");
   model.conditional_parameter_covariances =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/Russia_conditional_parameter_covariances.csv");
+          "../src/data/Russia_posterior_p_c.csv");
   model.free_energy = -8396.914141674371;
 }
 
@@ -393,10 +393,10 @@ void generate_russia_posteriors(dynamic_COVID_model& model) {
 void generate_mexico_posteriors(dynamic_COVID_model& model) {
   model.conditional_parameter_expectations =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/Mexico_conditional_parameter_expectations.csv");
+          "../src/data/Mexico_posterior_p_e.csv");
   model.conditional_parameter_covariances =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/Mexico_conditional_parameter_covariances.csv");
+          "../src/data/Mexico_posterior_p_c.csv");
   model.free_energy = -7379.296952775237;
 }
 
@@ -405,13 +405,13 @@ void generate_mexico_posteriors(dynamic_COVID_model& model) {
  */
 void generate_PEB_values(peb_model<dynamic_COVID_model>& model) {
   model.prior_parameter_expectations = utility::read_matrix<Eigen::MatrixXd>(
-      "../src/data/peb_prior_parameter_expectations.csv");
+      "../src/data/peb_prior_p_e.csv");
   model.prior_parameter_covariances = utility::read_matrix<Eigen::MatrixXd>(
-      "../src/data/peb_prior_parameter_covariances.csv");
+      "../src/data/peb_prior_p_c.csv");
   model.conditional_parameter_expectations =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/peb_conditional_parameter_expectations.csv");
+          "../src/data/peb_posterior_p_e.csv");
   model.conditional_parameter_covariances =
       utility::read_matrix<Eigen::MatrixXd>(
-          "../src/data/peb_conditional_parameter_covariances.csv");
+          "../src/data/peb_posterior_p_c.csv");
 }
diff --git a/src/COVID/import_COVID.cc b/src/COVID/import_COVID.cc
@@ -33,8 +33,8 @@ std::vector<country_data> read_country_data(int num_countries) {
   std::vector<country_data> countries;
   std::ifstream cases_file;
   std::ifstream deaths_file;
-  cases_file.open("../src/data/time_series_covid19_confirmed_global_new.csv");
-  deaths_file.open("../src/data/time_series_covid19_deaths_global_new.csv");
+  cases_file.open("../src/data/time_series_covid19_confirmed_global.csv");
+  deaths_file.open("../src/data/time_series_covid19_deaths_global.csv");
 
   std::string cases_line;
   std::string deaths_line;
@@ -52,8 +52,12 @@ std::vector<country_data> read_country_data(int num_countries) {
     splitstr(cases_split_tmp, cases_line, ',');
     splitstr(deaths_split_tmp, deaths_line, ',');
     int first_case = 4;
+    // First case is when cases are > 1 for parity with DCM for COVID example
+    // Should probably be >= 1
+    int cumsum = 0;
     for (int i = 4; i < cases_split_tmp.size(); i++) {
-      if (cases_split_tmp[i].compare("0")) {
+      cumsum += stoi(cases_split_tmp[i]);
+      if (cumsum > 1) {
         first_case = i;
         break;
       }

diff --git a/src/data/Brazil_conditional_parameter_covariances.csv b/src/data/Brazil_conditional_parameter_covariances.csv
diff --git a/src/data/Brazil_conditional_parameter_expectations.csv b/src/data/Brazil_conditional_parameter_expectations.csv