Get concentration/factors from most granular PKPD time to WH scope and print.

Author: b-s-code

model/Host/WithinHost/CommonWithinHost.cpp

@@ -238,6 +238,10 @@ void CommonWithinHost::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i
     // Inner vector has one entry "per drug" (really, per element of LSTMModel.m_drugs).
     // Each inner vector element is a drug id (integer), and a drug-specific drug factor.
     std::vector<std::vector<std::pair<size_t, double>>> identifiedComponentDrugFactors;
+
+    // One outer vector entry for each day for each infection.
+    std::vector<std::vector<std::tuple<std::string, double, double>>> pkpdTimeToDrugConcentrationMaps;
+    std::vector<std::vector<std::tuple<std::string, double, double>>> pkpdTimeToTotalFactorMaps;
 
     for( SimTime now = sim::ts0(), end = sim::ts0() + sim::oneTS(); now < end; now = now + sim::oneDay() ){
         // every day, medicate drugs, update each infection, then decay drugs
@@ -250,10 +254,16 @@ void CommonWithinHost::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i
             bool expires = ((*inf)->bloodStage() ? treatmentBlood : treatmentLiver);
 
             if( !expires ){     /* no expiry due to simple treatment model; do update */
-                // TODO : receive structure containing outputs of individual calculateDrugFactor calls here.
+                // Using vector of pairs instead of a map, because otherwise values with same time key will get overwritten.
+                std::vector<std::tuple<std::string, double, double>> pkpdTimeToDrugConcentrationMap;
+                std::vector<std::tuple<std::string, double, double>> pkpdTimeToTotalFactorMap;
+
                 // The values within this structure need to somehow be associated with the specific drug they are from.
                 // In the scenario in question, the particular override calculateDrugFactor being called is LSTMDrugOneComp::calculateDrugFactor.
-                const std::pair<double, std::vector<std::pair<size_t, double>>> output = pkpdModel.getDrugFactor(rng, *inf, body_mass);
+                const std::pair<double, std::vector<std::pair<size_t, double>>> output = pkpdModel.getDrugFactor(rng, *inf, body_mass,
+                    pkpdTimeToDrugConcentrationMap,
+                    pkpdTimeToTotalFactorMap
+                );
                 const double drugFactor = output.first;
                 // Note this drug factor can in some cases be derived from more than one drug type.
                 drugFactors.push_back(drugFactor);
@@ -289,9 +299,15 @@ void CommonWithinHost::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i
     // Outer vector's elements correspond to individual days/infections.
     // Inner vector's elements correspond to individual drugs.
     // Format could possibly be improved w.r.t. days/infections.
-    std::vector<std::vector<std::pair<std::string, double>>> namedComponentDrugFactors;
+    // This drug factor, as opposed to the other drug factors which get printed out,
+    // represents the product of all drugs' drug factors for the given human-day-infection.
+    std::vector<std::vector<std::pair<std::string, double>>> namedComponentDrugFactorProducts;
     for (auto outer_vec : identifiedComponentDrugFactors)
     {
+        if (!outer_vec.empty())
+        {
+            int x = 999; // Dummy instructions for breakpoint.
+        }
         std::vector<std::pair<std::string, double>> namedDrugFactors;
         for (auto inner_vec : outer_vec)
         {
@@ -300,7 +316,7 @@ void CommonWithinHost::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i
             const double drugFactor = inner_vec.second;
             namedDrugFactors.push_back(std::pair<std::string, double>{drugName, drugFactor});
         }
-        namedComponentDrugFactors.push_back(namedDrugFactors);
+        namedComponentDrugFactorProducts.push_back(namedDrugFactors);
     }
 
     json special_info = json::object();
@@ -310,7 +326,9 @@ void CommonWithinHost::update(Host::Human &human, LocalRng& rng, int &nNewInfs_i
     special_info["parasite_density"] = parasite_density;
     special_info["getDrugFactor"] = drugFactors;
     special_info["drugConcentrations"] = drugConcentrationMap;
-    special_info["namedComponentDrugFactors"] = namedComponentDrugFactors;
+    special_info["namedComponentDrugFactorProducts"] = namedComponentDrugFactorProducts;
+    special_info["pkpdTimeToDrugConcentrationMaps"] = pkpdTimeToDrugConcentrationMaps;
+    special_info["pkpdTimeToTotalFactorMaps"] = pkpdTimeToTotalFactorMaps;
     std::cout << special_info.dump() << std::endl;
 
     // As in AJTMH p22, cumulative_h (X_h + 1) doesn't include infections added

model/PkPd/Drug/LSTMDrug.h

@@ -78,8 +78,9 @@ class LSTMDrug {
      * @param body_mass Weight of patient in kg
      */
     virtual double calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
-    std::vector<std::pair<double, double>>& pkpdTimeToDrugConcentrationMap,
-    std::vector<std::pair<double, double>>& pkpdTimeToTotalFactorMap
+    const std::string& drugName,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToDrugConcentrationMap,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToTotalFactorMap
     ) const =0;
 
     /** Updates concentration variable and clears day's doses.

model/PkPd/Drug/LSTMDrugOneComp.cpp

@@ -50,8 +50,9 @@ double LSTMDrugOneComp::getConcentration(size_t index) const {
 // TODO: in high transmission, is this going to get called more often than updateConcentration?
 // When does it make sense to try to optimise (avoid doing decay calcuations here)?
 double LSTMDrugOneComp::calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
-    std::vector<std::pair<double, double>>& pkpdTimeToDrugConcentrationMap,
-    std::vector<std::pair<double, double>>& pkpdTimeToTotalFactorMap
+    const std::string& drugName,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToDrugConcentrationMap,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToTotalFactorMap
 ) const {
     if( concentration == 0.0 && doses.size() == 0 ) return 1.0; // nothing to do
 
@@ -73,27 +74,27 @@ double LSTMDrugOneComp::calculateDrugFactor(LocalRng& rng, WithinHost::CommonInf
         // we iteratate through doses in time order (since doses are sorted)
         if( time_conc.first < 1.0 /*i.e. today*/ ){
             if( time < time_conc.first ){
-                pkpdTimeToDrugConcentrationMap.push_back(std::pair<double, double>{time, concentration_today});
+                pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, concentration_today});
                 totalFactor *= drugPD.calcFactor( Kn, neg_elim_rate, &concentration_today, time_conc.first - time );
-                pkpdTimeToTotalFactorMap.push_back(std::pair<double, double>{time, totalFactor});
-                pkpdTimeToDrugConcentrationMap.push_back(std::pair<double, double>{time, concentration_today});
+                pkpdTimeToTotalFactorMap.push_back(std::tuple<std::string, double, double>{drugName, time, totalFactor});
+                pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, concentration_today});
                 time = time_conc.first;
             }else{ assert( time == time_conc.first ); }
             // add dose (instantaneous absorption):
             concentration_today += time_conc.second / (vol_dist * body_mass);
-            pkpdTimeToDrugConcentrationMap.push_back(std::pair<double, double>{time, concentration_today});
+            pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, concentration_today});
         }else/*i.e. tomorrow or later*/{
             break;
         }
     }
     if( time < 1.0 ){
         totalFactor *= drugPD.calcFactor( Kn, neg_elim_rate, &concentration_today, 1.0 - time );
-        pkpdTimeToTotalFactorMap.push_back(std::pair<double, double>{time, totalFactor});
-        pkpdTimeToDrugConcentrationMap.push_back(std::pair<double, double>{time, concentration_today});
+        pkpdTimeToTotalFactorMap.push_back(std::tuple<std::string, double, double>{drugName, time, totalFactor});
+        pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, concentration_today});
 
     }
 
-    pkpdTimeToTotalFactorMap.push_back(std::pair<double, double>{time, totalFactor});
+    pkpdTimeToTotalFactorMap.push_back(std::tuple<std::string, double, double>{drugName, time, totalFactor});
     return totalFactor; // Drug effect per day per drug per parasite
 }
 

model/PkPd/Drug/LSTMDrugOneComp.h

@@ -51,8 +51,9 @@ class LSTMDrugOneComp : public LSTMDrug {
     virtual double getConcentration(size_t index) const;
 
     virtual double calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
-    std::vector<std::pair<double, double>>& pkpdTimeToDrugConcentrationMap,
-    std::vector<std::pair<double, double>>& pkpdTimeToTotalFactorMap
+    const std::string& drugName,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToDrugConcentrationMap,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToTotalFactorMap
     ) const;
     virtual void updateConcentration (double body_mass);
 

model/PkPd/LSTMModel.cpp

@@ -133,21 +133,22 @@ double LSTMModel::getDrugConc (size_t drug_index) const{
     return c;
 }
 
-std::pair<double, std::vector<std::pair<size_t, double>>> LSTMModel::getDrugFactor (LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass) const{
+std::pair<double, std::vector<std::pair<size_t, double>>> LSTMModel::getDrugFactor (LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToDrugConcentrationMap,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToTotalFactorMap
+) const{
     double factor = 1.0; //no effect
 
     // Using vector of pairs instead of map in case multiple elements in m_drugs have same drug index.
     std::vector<std::pair<size_t, double>> drugIndexToDrugFactor;
 
     for( auto drug = m_drugs.begin(), end = m_drugs.end(); drug != end; ++drug ){
-        // Can be used by caller to determine drug name.
+        // Get name of current drug.
         const size_t drugIndex = (*drug)->getIndex();
-
-        // Using vector of pairs instead of a map, because otherwise values with same time key will get overwritten.
-        std::vector<std::pair<double, double>> pkpdTimeToDrugConcentrationMap;
-        std::vector<std::pair<double, double>> pkpdTimeToTotalFactorMap;
+        const std::string& drugName = OM::PkPd::LSTMDrugType::getDrugAbbrev(drugIndex);
 
         double drugFactor = (*drug)->calculateDrugFactor(rng, inf, body_mass,
+            drugName,
             pkpdTimeToDrugConcentrationMap,
             pkpdTimeToTotalFactorMap
         );

model/PkPd/LSTMModel.h

@@ -124,7 +124,10 @@ class LSTMModel {
      * Each time step, on each infection, the parasite density is multiplied by
      * the return value of this infection. The WithinHostModels are responsible
      * for clearing infections once the parasite density is negligible. */
-    std::pair<double, std::vector<std::pair<size_t, double>>> getDrugFactor (LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass) const;
+    std::pair<double, std::vector<std::pair<size_t, double>>> getDrugFactor (LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
+        std::vector<std::tuple<std::string, double, double>>& pkpdTimeToDrugConcentrationMap,
+        std::vector<std::tuple<std::string, double, double>>& pkpdTimeToTotalFactorMap
+    ) const;
 
     /** After any resident infections have been reduced by getDrugFactor(),
      * this function is called to update drug levels to their effective level