Add visibility of drugs in pkpd time to LSTMDrugThreeComp and remove exit guards.

Author: b-s-code

model/PkPd/Drug/LSTMDrugConversion.cpp

@@ -135,8 +135,6 @@ gsl_integration_workspace *gsl_intgr_conv_wksp = gsl_integration_workspace_alloc
 //NOTE: we "should" free, but mem-leaks at end of program aren't really important
 // gsl_integration_workspace_free (gsl_intgr_conv_wksp);
 double LSTMDrugConversion::calculateFactor(const Params_convFactor& p, double duration) const{
-    std::cout << "INSIDE LSTMDrugConversion::calculateFactor" << std::endl;
-    exit(0);
     gsl_function F;
     F.function = &func_convFactor;
     // gsl_function doesn't accept const; we re-apply const later

model/PkPd/Drug/LSTMDrugThreeComp.cpp

@@ -130,8 +130,6 @@ gsl_integration_workspace *gsl_intgr_wksp = gsl_integration_workspace_alloc (GSL
 // gsl_integration_workspace_free (gsl_intgr_wksp);
 
 double LSTMDrugThreeComp::calculateFactor(const Params_fC& p, double duration) const{
-    std::cout << "INSIDE LSTMDrugThreeComp::calculateFactor" << std::endl;
-    exit(0);
     gsl_function F;
     F.function = &func_fC;
     // gsl_function doesn't accept const; we re-apply const later
@@ -160,7 +158,11 @@ double LSTMDrugThreeComp::calculateFactor(const Params_fC& p, double duration) c
 
 // 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 LSTMDrugThreeComp::calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass) const {
+double LSTMDrugThreeComp::calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
+    const std::string& drugName,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToDrugConcentrationMap,
+    std::vector<std::tuple<std::string, double, double>>& pkpdTimeToTotalFactorMap
+) const {
     if( conc() == 0.0 && doses.size() == 0 ) return 1.0; // nothing to do
     updateCached(body_mass);
 
@@ -180,26 +182,61 @@ double LSTMDrugThreeComp::calculateDrugFactor(LocalRng& rng, WithinHost::CommonI
         if( time_conc.first < 1.0 /*i.e. today*/ ){
             if( time < time_conc.first ){
                 double duration = time_conc.first - time;
+                
+                // Add concentration and factor tracking before duration
+                double current_conc = p.cA + p.cB + p.cC - p.cABC;
+                pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, current_conc});
+                
                 totalFactor *= calculateFactor(p, duration);
+                pkpdTimeToTotalFactorMap.push_back(std::tuple<std::string, double, double>{drugName, time, totalFactor});
+                
+                // Update concentrations after duration
                 p.cA *= exp(p.na * duration);
                 p.cB *= exp(p.nb * duration);
                 p.cC *= exp(p.ng * duration);
                 p.cABC *= exp(p.nka * duration);
+                
+                // Add concentration tracking after duration
+                current_conc = p.cA + p.cB + p.cC - p.cABC;
+                pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, current_conc});
+                
                 time = time_conc.first;
             }else{ assert( time == time_conc.first ); }
             // add dose:
             p.cA += A * time_conc.second;
             p.cB += B * time_conc.second;
             p.cC += C * time_conc.second;
             p.cABC += (A + B + C) * time_conc.second;
+            
+            // Add concentration tracking after dose
+            double current_conc = p.cA + p.cB + p.cC - p.cABC;
+            pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, current_conc});
         }else /*i.e. tomorrow or later*/{
             // ignore
         }
     }
     if( time < 1.0 ){
-        totalFactor *= calculateFactor(p, 1.0 - time);
+        double duration = 1.0 - time;
+        
+        // // Add concentration and factor tracking before final duration
+        // double current_conc = p.cA + p.cB + p.cC - p.cABC;
+        // pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, current_conc});
+        
+        totalFactor *= calculateFactor(p, duration);
+        pkpdTimeToTotalFactorMap.push_back(std::tuple<std::string, double, double>{drugName, time, totalFactor});
+        
+        // // Add concentration tracking after final duration
+        // p.cA *= exp(p.na * duration);
+        // p.cB *= exp(p.nb * duration);
+        // p.cC *= exp(p.ng * duration);
+        // p.cABC *= exp(p.nka * duration);
+        const double current_conc = p.cA + p.cB + p.cC - p.cABC;
+        pkpdTimeToDrugConcentrationMap.push_back(std::tuple<std::string, double, double>{drugName, time, current_conc});
     }
 
+    // Add final total factor entry
+    pkpdTimeToTotalFactorMap.push_back(std::tuple<std::string, double, double>{drugName, time, totalFactor});
+    
     return totalFactor;
 }
 

model/PkPd/Drug/LSTMDrugThreeComp.h

@@ -53,7 +53,11 @@ class LSTMDrugThreeComp : public LSTMDrug {
     virtual size_t getIndex() const;
     virtual double getConcentration(size_t index) const;
 
-    virtual double calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass) const;
+    virtual double calculateDrugFactor(LocalRng& rng, WithinHost::CommonInfection *inf, double body_mass,
+    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);
 
 protected:

model/PkPd/Drug/LSTMDrugType.cpp

@@ -155,9 +155,7 @@ unique_ptr<LSTMDrug> LSTMDrugType::createInstance(LocalRng& rng, size_t index) {
         throw std::runtime_error("Drug conversion not supported");
     }else if( typeData.k12 ){
         // k21 is set when k12 is set; k13 and k31 may be set
-        // return unique_ptr<LSTMDrug>(new LSTMDrugThreeComp( typeData, rng ));
-        // High velocity hack.
-        throw std::runtime_error("Multi-compartment drug models not supported");
+        return unique_ptr<LSTMDrug>(new LSTMDrugThreeComp( typeData, rng ));
     }else{
         // none of k12/k21/k13/k31 should be set in this case
         return unique_ptr<LSTMDrug>(new LSTMDrugOneComp( typeData, rng ));