add support for Gamma CV=0 and LogNormal sigma=0

Author: acavelan

model/Transmission/PerHost.cpp

@@ -72,13 +72,15 @@ PerHost::PerHost () :
 
 void PerHost::initialise (LocalRng& rng, double availabilityFactor) {
     relativeAvailabilityHet = availabilityFactor;
+    anophEntoAvailabilityRaw.resize(PerHostAnophParams::numSpecies());
     anophEntoAvailability.resize(PerHostAnophParams::numSpecies());
     anophProbMosqBiting.resize(PerHostAnophParams::numSpecies());
     anophProbMosqResting.resize(PerHostAnophParams::numSpecies());
 
     for(size_t i = 0; i < PerHostAnophParams::numSpecies(); ++i) {
         const PerHostAnophParams& base = PerHostAnophParams::get(i);
-        anophEntoAvailability[i] = base.entoAvailabilityFactor * std::min(base.entoAvailability->sample(rng), 25.0) * availabilityFactor;
+        anophEntoAvailabilityRaw[i] = std::min(base.entoAvailability->sample(rng), 25.0);
+        anophEntoAvailability[i] = base.entoAvailabilityFactor * anophEntoAvailabilityRaw[i] * availabilityFactor;
         anophProbMosqBiting[i] = base.probMosqBiting.sample(rng);
         auto pRest1 = base.probMosqFindRestSite.sample(rng);
         auto pRest2 = base.probMosqSurvivalResting.sample(rng);

model/Transmission/PerHost.h

@@ -337,6 +337,7 @@ class PerHost
     /// Checkpointing
     template<class S>
     void operator& (S& stream) {
+        anophEntoAvailabilityRaw & stream;
         anophEntoAvailability & stream;
         anophProbMosqBiting & stream;
         anophProbMosqResting & stream;
@@ -354,7 +355,10 @@ class PerHost
     // entoAvailability param stored in HostMosquitoInteraction.
     double relativeAvailabilityHet;
 
-    /** Species availability rate of human to mosquitoes, including hetergeneity factor
+    /** Species availability rate of human to mosquitoes */
+    vector<double> anophEntoAvailabilityRaw;
+
+    /** Species availability rate of human to mosquitoes, including heterogeneity factor
      * and base rate, but excluding age and intervention factors. */
     vector<double> anophEntoAvailability;
 

model/interventions/Deployments.hpp

@@ -253,50 +253,50 @@ class TimedHumanDeployment : public TimedDeployment, protected HumanDeploymentBa
 
                 if(copula)
                 {
-                    if(Transmission::PerHostAnophParams::numSpecies() > 1)
-                        throw std::runtime_error("Gaussian copula: (deployment with availability correlation) only supports one mosquito species");
-
-                    // Use avail that is already weighted?? or Need raw value?
-                    double scaling_factor = (Transmission::PerHostAnophParams::get(0).entoAvailabilityFactor * human.perHostTransmission.relativeAvailabilityHet);
-                    double ux = Transmission::PerHostAnophParams::get(0).entoAvailability->cdf(availability / scaling_factor);
-                    double xx = gsl_cdf_ugaussian_Pinv(ux);                      // Unit interval to Normal
-                
-                    // Apply the Gaussian copula transformation for correlation
-                    ComponentId cid = subPop;
-
-                    /** human.rng.gauss(0.0, 1.0) is calculated once per component and per human. Re-deployment of
-                     * the same component on the same human must use the same gaussian sample. Therefore we store
-                     * this value in the human the first time it is calculated. */
-                    double g = 0.0;
-                    auto it = human.perHostTransmission.copulaGaussianSamples.find(cid);
-                    if (it != human.perHostTransmission.copulaGaussianSamples.end())
-                        g = it->second;
-                    else
-                    {
-                        g = human.rng.gauss(0.0, 1.0);
-                        human.perHostTransmission.copulaGaussianSamples[cid] = g;
-                    }
-
-                    double yy = coverageCorr * xx + g * sqrt(1 - coverageCorr * coverageCorr);
-
-                    // Transform back to unit interval
-                    double uy = gsl_cdf_ugaussian_P(yy);
-        
                     // Beta distribution for intervention (Beta distributed)
                     double beta_mean = coverage;  // Assuming this value is given
                     double beta_var = coverageVar;       // Given as well
                     double alpha = ((1.0 - beta_mean) / beta_var - 1.0 / beta_mean) * (beta_mean * beta_mean);
                     double beta = alpha * (1.0 / beta_mean - 1.0);
 
-                    if (alpha < 0 || beta < 0)
-                        throw std::runtime_error("Gaussian copula: resulting alpha and beta parameters must be positive");
+                    if(Transmission::PerHostAnophParams::numSpecies() > 1)
+                        throw std::runtime_error("Gaussian copula: (deployment with availability correlation) only supports one mosquito species");
 
-                    
+                    double ux = Transmission::PerHostAnophParams::get(0).entoAvailability->cdf(human.perHostTransmission.anophEntoAvailabilityRaw[0]);
+                    if(ux < 1)
+                    {
+                        double xx = gsl_cdf_ugaussian_Pinv(ux);                      // Unit interval to Normal
+
+                        // Apply the Gaussian copula transformation for correlation
+                        ComponentId cid = subPop;
+
+                        /** human.rng.gauss(0.0, 1.0) is calculated once per component and per human. Re-deployment of
+                         * the same component on the same human must use the same gaussian sample. Therefore we store
+                         * this value in the human the first time it is calculated. */
+                        double g = 0.0;
+                        auto it = human.perHostTransmission.copulaGaussianSamples.find(cid);
+                        if (it != human.perHostTransmission.copulaGaussianSamples.end())
+                            g = it->second;
+                        else
+                        {
+                            g = human.rng.gauss(0.0, 1.0);
+                            human.perHostTransmission.copulaGaussianSamples[cid] = g;
+                        }
+                        double yy = coverageCorr * xx + g * sqrt(1 - coverageCorr * coverageCorr);
 
-                    // Get the final probability from Beta distribution
-                    probability = gsl_cdf_beta_Pinv(uy, alpha, beta);
+                        // Transform back to unit interval
+                        double uy = gsl_cdf_ugaussian_P(yy);
+
+                        if (alpha < 0 || beta < 0)
+                            throw std::runtime_error("Gaussian copula: resulting alpha and beta parameters must be positive");
+
+                        // Get the final probability from Beta distribution
+                        probability = gsl_cdf_beta_Pinv(uy, alpha, beta);
+                    }
+                    else // Gamma with CV=0
+                        probability = human.rng.beta(alpha, beta);
 
-                    //cout << availability << " " << ux << " " << xx << " " << yy << " " << uy << " " << coverageCorr << " " << g << endl;;
+                    // cout << availability << " " << probability << " " << alpha << " " << beta << endl;
 
                     file << availability << " " << probability << endl;
                 }

model/util/sampler.cpp

@@ -171,10 +171,11 @@ double LognormalSampler::sample(LocalRng& rng) const{
 }
 
 double LognormalSampler::cdf(double x) const {
-    // Check if sigma is zero, and handle the degenerate log-normal case
     if (sigma == 0.0)
-        throw std::runtime_error("LognormalSampler::cdf() does not support sigma = 0 (default const distirbution)");
-
+    {
+        if(x >= mu) return 1.0;
+        else return 0.0;
+    }
     return gsl_cdf_lognormal_P(x, mu, sigma);
 }
 
@@ -266,12 +267,17 @@ double GammaSampler::mean() const {
 }
 
 double GammaSampler::sample(LocalRng& rng) const{
-    if(isnan(theta))
+    if(isnan(theta)) // CV=0
         return mu;
     return rng.gamma(k, theta);
 }
 
 double GammaSampler::cdf(double x) const {
+    if(isnan(theta)) // CV=0
+    {
+        if(x >= mu) return 1.0;
+        else return 0.0;
+    }
     return gsl_cdf_gamma_P(x, k, theta);
 }