deleted file mode 100644

@@ -1,255 +0,0 @@

-//     Copyright 2013-2021 Ramon Diaz-Uriarte

-

-//     This program is free software: you can redistribute it and/or modify

-//     it under the terms of the GNU General Public License as published by

-//     the Free Software Foundation, either version 3 of the License, or

-//     (at your option) any later version.

-

-//     This program is distributed in the hope that it will be useful,

-//     but WITHOUT ANY WARRANTY; without even the implied warranty of

-//     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

-//     GNU General Public License for more details.

-

-//     You should have received a copy of the GNU General Public License

-//     along with this program.  If not, see <http://www.gnu.org/licenses/>.

-

-

-

-#ifndef _BNB_COMMON_H_

-#define _BNB_COMMON_H_

-

-#include <Rcpp.h>

-#include "common_classes.h"

-#include "debug_common.h"

-#include "new_restrict.h" // for the TypeModel enum

-

-

-inline void W_f_st(spParamsP& spP){

-  spP.W = spP.death + spP.birth + spP.mutation;

-}

-

-inline void R_f_st(spParamsP& spP) {

-  spP.R = sqrt( pow( spP.birth - spP.death, 2) +

-		( 2.0 * (spP.birth + spP.death) +

-		  spP.mutation) * spP.mutation );

-}

-

-

-inline double pE_f_st(double& pM, const spParamsP& spP){

-  double pE = (spP.death * (1.0 - pM ) )/(spP.W - spP.death - spP.birth * pM );

-  if( !std::isfinite(pE) ) {

-    DP2(spP.death);  DP2(spP.birth); DP2(pM); DP2(spP.W);

-    DP2(spP.mutation);

-    std::string error_message = R"(pE.f: pE not finite.

-      This is expected to happen when mutationPropGrowth = TRUE

-      and you have have an initMutant with death >> birth,

-      as that inevitably leads to net birth rate of 0

-      and mutation rate of 0)";

-    throw std::range_error(error_message);

-  }

-  return pE;

-}

-

-inline double pB_f_st(const double& pE,

-		      const spParamsP& spP) {

-  return (spP.birth * pE)/spP.death;

-}

-

-

-void print_mapTimes(std::multimap<double, int>& mapTimes);

-  

-void print_initMutant(const std::vector < std::vector<int> >& initMutant);

-

-void print_Genotype(const Genotype& ge);

-

-void mapTimes_updateP(std::multimap<double, int>& mapTimes,

-		      std::vector<spParamsP>& popParams,

-		      const int index,

-		      const double time);

-

-

-void getMinNextMutationTime4(int& nextMutant, double& minNextMutationTime,

-			     const std::multimap<double, int>& mapTimes);

-

-

-void fill_SStats(Rcpp::NumericMatrix& perSampleStats,

-		 const std::vector<double>& sampleTotPopSize,

-		 const std::vector<double>& sampleLargestPopSize,

-		 const std::vector<double>& sampleLargestPopProp,

-		 const std::vector<int>& sampleMaxNDr,

-		 const std::vector<int>& sampleNDrLargestPop);

-

-

-

-void print_spP(const spParamsP& spP);

-

-double pM_f_st(const double& t, const spParamsP& spP);

-

-double ti_nextTime_tmax_2_st(const spParamsP& spP,

-			     const double& currentTime,

-			     const double& tSample,

-			     int& ti_dbl_min,

-			     int& ti_e3);

-

-double Algo2_st(const spParamsP& spP,

-		const double& ti,

-		const int& mutationPropGrowth);

-

-double Algo3_st(const spParamsP& spP, const double& t);

-

-void precissionLoss();

-

-void init_tmpP(spParamsP& tmpParam);

-

-double returnMFE_new(double& en1,

-		     const TypeModel typeModel);

-

-

-void computeMcFarlandError_new(double& en1,

-			       double& en1sc,

-			       double& totPopSize_previous,

-			       double& DA_previous,

-			       const TypeModel typeModel,

-			       const double& totPopSize,

-			       const double& K);

-

-void updateRatesMcFarland(std::vector<spParamsP>& popParams,

-			  double& adjust_fitness_MF,

-			  const double& K,

-			  const double& totPopSize);

-

-void updateRatesMcFarlandLog(std::vector<spParamsP>& popParams,

-			     double& adjust_fitness_MF,

-			     const double& K,

-			     const double& totPopSize);

-

-void updateRatesFDFMcFarlandLog(std::vector<spParamsP>& popParams,

-				const std::vector<Genotype>& Genotypes,

-				const fitnessEffectsAll& fitnessEffects,

-				double& adjust_fitness_MF,

-				const double& K,

-				const double& totPopSize,

-				const double& currentTime);

-

-void updateRatesMcFarlandLog_D(std::vector<spParamsP>& popParams,

-			       double& adjust_fitness_MF,

-			       const double& K,

-			       const double& totPopSize);

-

-void updateRatesFDFMcFarlandLog_D(std::vector<spParamsP>& popParams,

-				  const std::vector<Genotype>& Genotypes,

-				  const fitnessEffectsAll& fitnessEffects,

-				  double& adjust_fitness_MF,

-				  const double& K,

-				  const double& totPopSize,

-				  const double& currentTime);

-

-

-void updateRatesFDFExp(std::vector<spParamsP>& popParams,

-		       const std::vector<Genotype>& Genotypes,

-		       const fitnessEffectsAll& fitnessEffects,

-		       const double& currentTime);

-

-void updateRatesFDFBozic(std::vector<spParamsP>& popParams,

-			 const std::vector<Genotype>& Genotypes,

-			 const fitnessEffectsAll& fitnessEffects,

-			 const double& currentTime);

-

-

-void updateRatesArbitrary(std::vector<spParamsP>& popParams,

-		       const std::vector<Genotype>& Genotypes,

-		       const fitnessEffectsAll& fitnessEffects,

-		       const double& currentTime);

-

-

-void updateRatesConstant(std::vector<spParamsP>& popParams,

-		       const std::vector<Genotype>& Genotypes,

-		       const fitnessEffectsAll& fitnessEffects,

-               const double& cte,

-           	   const double& sampleEvery,

-		       const double& currentTime);

-

-void updateBirthDeathRates(std::vector<spParamsP>& popParams,

-			   const std::vector<Genotype>& Genotypes,

-			   const fitnessEffectsAll& fitnessEffects,

-			   double& adjust_fitness_MF,

-			   const double& K,

-			   const double& totPopSize,

-			   const double& cteSize,

-			   const double& sampleEvery,

-			   const double& currentTime,

-			   const TypeModel typeModel);

-

-void detect_ti_duplicates(const std::multimap<double, int>& m,

-			  const double ti,

-			  const int spcies);

-

-

-

-void message1(const int verbosity, const std::string message,

-	      const int iteration, const double currentTime,

-	      const unsigned int numSpecies,

-	      const double totalPopulationSize,

-	      const double timeNextPopSample,

-	      const double minNextMutationTime);

-

-void messageNewSpecies(const int verbosity,

-		       const int iteration, 

-		       const unsigned int numSpecies,

-		       const int nextMutant);

-

-void vvmessageNewSpecies(const int verbosity,

-			 const unsigned int sp,

-			 const Genotype& newGenotype,

-			 const Genotype& parentGenotype,

-			 const spParamsP& tmpParam,

-			 const spParamsP& parentParam);

-

-void messageSampling(const int verbosity,

-		     const double tSample,

-		     const double finalTime,

-		     std::vector<spParamsP>& popParams);

-

-void messagePostSampling(const int verbosity,

-			 std::vector<spParamsP>& popParams,

-			 const double totPopSize);

-

-double setDummyMutationRate(std::vector<double> mu, const int verbosity);

-

-

-void initPops(

-	     unsigned int& numSpecies,

-	     double& totPopSize,

-	     int& outNS_i,

-	     double& lastStoredSample,

-	     std::vector<Genotype>& Genotypes,

-	     std::vector<spParamsP>& popParams,

-	     std::vector<Genotype>& genot_out,

-	     std::vector<double>& popSizes_out,

-	     std::vector<int>& index_out,

-	     std::vector<double>& time_out,

-	     std::vector<double>& sampleTotPopSize,

-	     std::vector<double>& sampleLargestPopSize,

-	     std::vector<int>& sampleMaxNDr,

-	     std::vector<int>& sampleNDrLargestPop,

-	     POM& pom,

-	     std::mt19937& ran_gen,

-	     const std::vector<std::vector<int> >& initMutant,

-	     const std::vector<double>& initSize,

-	     const fitnessEffectsAll& fitnessEffects,

-	     const std::vector<double>& mu,

-	     const fitnessEffectsAll& muEF,

-	     const std::vector<int>& full2mutator,

-	     const std::map<int, std::string>& intName,

-	     const fitness_as_genes& genesInFitness,	     

-	     const double& dummyMutationRate,

-	     const double& K,

-	     const double& death,

-	     const double& currentTime,	     

-	     const double& keepEvery,

-	     const int& mutationPropGrowth,

-	     const TypeModel typeModel,

-	     const int& verbosity

-	      );

-

-#endif

@@ -18,7 +18,7 @@

 #ifndef _BNB_COMMON_H_

 #define _BNB_COMMON_H_

-#include<Rcpp.h>

+#include <Rcpp.h>

 #include "common_classes.h"

 #include "debug_common.h"

 #include "new_restrict.h" // for the TypeModel enum

@@ -156,12 +156,27 @@ void updateRatesFDFBozic(std::vector<spParamsP>& popParams,

 			 const double& currentTime);

+void updateRatesArbitrary(std::vector<spParamsP>& popParams,

+		       const std::vector<Genotype>& Genotypes,

+		       const fitnessEffectsAll& fitnessEffects,

+		       const double& currentTime);

+

+

+void updateRatesConstant(std::vector<spParamsP>& popParams,

+		       const std::vector<Genotype>& Genotypes,

+		       const fitnessEffectsAll& fitnessEffects,

+               const double& cte,

+           	   const double& sampleEvery,

+		       const double& currentTime);

+

 void updateBirthDeathRates(std::vector<spParamsP>& popParams,

 			   const std::vector<Genotype>& Genotypes,

 			   const fitnessEffectsAll& fitnessEffects,

 			   double& adjust_fitness_MF,

 			   const double& K,

 			   const double& totPopSize,

+			   const double& cteSize,

+			   const double& sampleEvery,

 			   const double& currentTime,

 			   const TypeModel typeModel);

@@ -232,7 +247,7 @@ void initPops(

 	     const double& death,

 	     const double& currentTime,	     

 	     const double& keepEvery,

-	     const int& mutationPropGrowth,	     

+	     const int& mutationPropGrowth,

 	     const TypeModel typeModel,

 	     const int& verbosity

 	      );

@@ -101,9 +101,6 @@ void precissionLoss();

 void init_tmpP(spParamsP& tmpParam);

-double returnMFE_new(double& en1,

-		     const std::string& typeFitness);

-

 double returnMFE_new(double& en1,

 		     const TypeModel typeModel);

@@ -116,14 +113,6 @@ void computeMcFarlandError_new(double& en1,

 			       const double& totPopSize,

 			       const double& K);

-void computeMcFarlandError_new(double& en1,

-			       double& en1sc,

-			       double& totPopSize_previous,

-			       double& DA_previous,

-			       const std::string& typeFitness,

-			       const double& totPopSize,

-			       const double& K);

-

 void updateRatesMcFarland(std::vector<spParamsP>& popParams,

 			  double& adjust_fitness_MF,

 			  const double& K,

@@ -23,27 +23,6 @@

 #include "debug_common.h"

 #include "new_restrict.h" // for the TypeModel enum

-// // Simple custom exception for exceptions that lead to re-runs.

-// class rerunExcept: public std::runtime_error {

-// public:

-//   rerunExcept(const std::string &s) :

-//     std::runtime_error(s) {}

-// };

-

-

-// struct spParamsP {

-//   double popSize;

-//   double birth;

-//   double death;

-//   double W;

-//   double R;

-//   double mutation;

-//   double timeLastUpdate;

-//   std::multimap<double, int>::iterator pv;

-//   double absfitness; //convenient for Beerenwinkel

-//   int numMutablePos; //for mutator if need update of mutation

-// };

-

 inline void W_f_st(spParamsP& spP){

   spP.W = spP.death + spP.birth + spP.mutation;

@@ -72,14 +51,21 @@ inline double pE_f_st(double& pM, const spParamsP& spP){

 }

 inline double pB_f_st(const double& pE,

-			     const spParamsP& spP) {

+		      const spParamsP& spP) {

   return (spP.birth * pE)/spP.death;

 }

+

+void print_mapTimes(std::multimap<double, int>& mapTimes);

+  

+void print_initMutant(const std::vector < std::vector<int> >& initMutant);

+

+void print_Genotype(const Genotype& ge);

+

 void mapTimes_updateP(std::multimap<double, int>& mapTimes,

-			     std::vector<spParamsP>& popParams,

-			     const int index,

-			     const double time);

+		      std::vector<spParamsP>& popParams,

+		      const int index,

+		      const double time);

 void getMinNextMutationTime4(int& nextMutant, double& minNextMutationTime,

@@ -87,15 +73,13 @@ void getMinNextMutationTime4(int& nextMutant, double& minNextMutationTime,

 void fill_SStats(Rcpp::NumericMatrix& perSampleStats,

-			       const std::vector<double>& sampleTotPopSize,

-			       const std::vector<double>& sampleLargestPopSize,

-			       const std::vector<double>& sampleLargestPopProp,

-			       const std::vector<int>& sampleMaxNDr,

+		 const std::vector<double>& sampleTotPopSize,

+		 const std::vector<double>& sampleLargestPopSize,

+		 const std::vector<double>& sampleLargestPopProp,

+		 const std::vector<int>& sampleMaxNDr,

 		 const std::vector<int>& sampleNDrLargestPop);

-void print_mapTimes(std::multimap<double, int>& mapTimes);

-  

-void print_initMutant(const std::vector < std::vector<int> >& initMutant);

+

 void print_spP(const spParamsP& spP);

@@ -117,31 +101,12 @@ void precissionLoss();

 void init_tmpP(spParamsP& tmpParam);

-// double returnMFE(double& e1,

-// 		 const std::string& typeFitness);

-

-// double returnMFE(double& e1,

-// 		 const TypeModel typeModel);

-

 double returnMFE_new(double& en1,

 		     const std::string& typeFitness);

 double returnMFE_new(double& en1,

 		     const TypeModel typeModel);

-// void computeMcFarlandError(double& e1,

-// 			   double& n_0,

-// 			   double& tps_0,

-// 			   const std::string& typeFitness,

-// 			   const double& totPopSize,

-// 			   const double& K);

-

-// void computeMcFarlandError(double& e1,

-// 			   double& n_0,

-// 			   double& tps_0,

-// 			   const TypeModel typeModel,

-// 			   const double& totPopSize,

-// 			   const double& K);

 void computeMcFarlandError_new(double& en1,

 			       double& en1sc,

@@ -160,65 +125,127 @@ void computeMcFarlandError_new(double& en1,

 			       const double& K);

 void updateRatesMcFarland(std::vector<spParamsP>& popParams,

-				 double& adjust_fitness_MF,

-				 const double& K,

+			  double& adjust_fitness_MF,

+			  const double& K,

 			  const double& totPopSize);

 void updateRatesMcFarlandLog(std::vector<spParamsP>& popParams,

-				    double& adjust_fitness_MF,

-				    const double& K,

+			     double& adjust_fitness_MF,

+			     const double& K,

 			     const double& totPopSize);

 void updateRatesFDFMcFarlandLog(std::vector<spParamsP>& popParams,

-  const std::vector<Genotype>& Genotypes,

-  const fitnessEffectsAll& fitnessEffects,

-  double& adjust_fitness_MF,

-  const double& K,

-  const double& totPopSize,

-  const double& currentTime);

+				const std::vector<Genotype>& Genotypes,

+				const fitnessEffectsAll& fitnessEffects,

+				double& adjust_fitness_MF,

+				const double& K,

+				const double& totPopSize,

+				const double& currentTime);

 void updateRatesMcFarlandLog_D(std::vector<spParamsP>& popParams,

-				    double& adjust_fitness_MF,

-				    const double& K,

-			     const double& totPopSize);

+			       double& adjust_fitness_MF,

+			       const double& K,

+			       const double& totPopSize);

 void updateRatesFDFMcFarlandLog_D(std::vector<spParamsP>& popParams,

-  const std::vector<Genotype>& Genotypes,

-  const fitnessEffectsAll& fitnessEffects,

-  double& adjust_fitness_MF,

-  const double& K,

-  const double& totPopSize,

-  const double& currentTime);

+				  const std::vector<Genotype>& Genotypes,

+				  const fitnessEffectsAll& fitnessEffects,

+				  double& adjust_fitness_MF,

+				  const double& K,

+				  const double& totPopSize,

+				  const double& currentTime);

 void updateRatesFDFExp(std::vector<spParamsP>& popParams,

-  const std::vector<Genotype>& Genotypes,

-  const fitnessEffectsAll& fitnessEffects,

-  const double& currentTime);

+		       const std::vector<Genotype>& Genotypes,

+		       const fitnessEffectsAll& fitnessEffects,

+		       const double& currentTime);

 void updateRatesFDFBozic(std::vector<spParamsP>& popParams,

-  const std::vector<Genotype>& Genotypes,

-  const fitnessEffectsAll& fitnessEffects,

-  const double& currentTime);

+			 const std::vector<Genotype>& Genotypes,

+			 const fitnessEffectsAll& fitnessEffects,

+			 const double& currentTime);

-void updateRatesMcFarland0(std::vector<spParamsP>& popParams,

-				  double& adjust_fitness_MF,

-				  const double& K,

-				  const double& totPopSize,

-				  const int& mutationPropGrowth,

-			   const double& mu);

-

-void updateRatesBeeren(std::vector<spParamsP>& popParams,

-			      double& adjust_fitness_B,

-			      const double& initSize,

-			      const double& currentTime,

-			      const double& alpha,

-			      const double& totPopSize,

-			      const int& mutationPropGrowth,

-		       const double& mu);

+

+void updateBirthDeathRates(std::vector<spParamsP>& popParams,

+			   const std::vector<Genotype>& Genotypes,

+			   const fitnessEffectsAll& fitnessEffects,

+			   double& adjust_fitness_MF,

+			   const double& K,

+			   const double& totPopSize,

+			   const double& currentTime,

+			   const TypeModel typeModel);

 void detect_ti_duplicates(const std::multimap<double, int>& m,

 			  const double ti,

 			  const int spcies);

+

+

+void message1(const int verbosity, const std::string message,

+	      const int iteration, const double currentTime,

+	      const unsigned int numSpecies,

+	      const double totalPopulationSize,

+	      const double timeNextPopSample,

+	      const double minNextMutationTime);

+

+void messageNewSpecies(const int verbosity,

+		       const int iteration, 

+		       const unsigned int numSpecies,

+		       const int nextMutant);

+

+void vvmessageNewSpecies(const int verbosity,

+			 const unsigned int sp,

+			 const Genotype& newGenotype,

+			 const Genotype& parentGenotype,

+			 const spParamsP& tmpParam,

+			 const spParamsP& parentParam);

+

+void messageSampling(const int verbosity,

+		     const double tSample,

+		     const double finalTime,

+		     std::vector<spParamsP>& popParams);

+

+void messagePostSampling(const int verbosity,

+			 std::vector<spParamsP>& popParams,

+			 const double totPopSize);

+

+double setDummyMutationRate(std::vector<double> mu, const int verbosity);

+

+

+void initPops(

+	     unsigned int& numSpecies,

+	     double& totPopSize,

+	     int& outNS_i,

+	     double& lastStoredSample,

+	     std::vector<Genotype>& Genotypes,

+	     std::vector<spParamsP>& popParams,

+	     std::vector<Genotype>& genot_out,

+	     std::vector<double>& popSizes_out,

+	     std::vector<int>& index_out,

+	     std::vector<double>& time_out,

+	     std::vector<double>& sampleTotPopSize,

+	     std::vector<double>& sampleLargestPopSize,

+	     std::vector<int>& sampleMaxNDr,

+	     std::vector<int>& sampleNDrLargestPop,

+	     POM& pom,

+	     std::mt19937& ran_gen,

+	     const std::vector<std::vector<int> >& initMutant,

+	     const std::vector<double>& initSize,

+	     const fitnessEffectsAll& fitnessEffects,

+	     const std::vector<double>& mu,

+	     const fitnessEffectsAll& muEF,

+	     const std::vector<int>& full2mutator,

+	     const std::map<int, std::string>& intName,

+	     const fitness_as_genes& genesInFitness,	     

+	     const double& dummyMutationRate,

+	     const double& K,

+	     const double& death,

+	     const double& currentTime,	     

+	     const double& keepEvery,

+	     const int& mutationPropGrowth,	     

+	     const TypeModel typeModel,

+	     const int& verbosity

+	      );

+

 #endif

@@ -1,4 +1,4 @@

-//     Copyright 2013, 2014, 2015 Ramon Diaz-Uriarte

+//     Copyright 2013-2021 Ramon Diaz-Uriarte

 //     This program is free software: you can redistribute it and/or modify

 //     it under the terms of the GNU General Public License as published by

@@ -21,7 +21,7 @@

 #include<Rcpp.h>

 #include "common_classes.h"

 #include "debug_common.h"

-// #include "new_restrict.h" // for the TypeModel enum

+#include "new_restrict.h" // for the TypeModel enum

 // // Simple custom exception for exceptions that lead to re-runs.

 // class rerunExcept: public std::runtime_error {

@@ -37,7 +37,7 @@

 //   double death;

 //   double W;

 //   double R;

-//   double mutation; 

+//   double mutation;

 //   double timeLastUpdate;

 //   std::multimap<double, int>::iterator pv;

 //   double absfitness; //convenient for Beerenwinkel

@@ -50,8 +50,8 @@ inline void W_f_st(spParamsP& spP){

 }

 inline void R_f_st(spParamsP& spP) {

-  spP.R = sqrt( pow( spP.birth - spP.death, 2) + 

-		( 2.0 * (spP.birth + spP.death) + 

+  spP.R = sqrt( pow( spP.birth - spP.death, 2) +

+		( 2.0 * (spP.birth + spP.death) +

 		  spP.mutation) * spP.mutation );

 }

@@ -62,7 +62,7 @@ inline double pE_f_st(double& pM, const spParamsP& spP){

     DP2(spP.death);  DP2(spP.birth); DP2(pM); DP2(spP.W);

     DP2(spP.mutation);

     std::string error_message = R"(pE.f: pE not finite.

-      This is expected to happen when mutationPropGrowth = TRUE 

+      This is expected to happen when mutationPropGrowth = TRUE

       and you have have an initMutant with death >> birth,

       as that inevitably leads to net birth rate of 0

       and mutation rate of 0)";

@@ -73,7 +73,7 @@ inline double pE_f_st(double& pM, const spParamsP& spP){

 inline double pB_f_st(const double& pE,

 			     const spParamsP& spP) {

-  return (spP.birth * pE)/spP.death; 

+  return (spP.birth * pE)/spP.death;

 }

 void mapTimes_updateP(std::multimap<double, int>& mapTimes,

@@ -93,6 +93,9 @@ void fill_SStats(Rcpp::NumericMatrix& perSampleStats,

 			       const std::vector<int>& sampleMaxNDr,

 		 const std::vector<int>& sampleNDrLargestPop);

+void print_mapTimes(std::multimap<double, int>& mapTimes);

+  

+void print_initMutant(const std::vector < std::vector<int> >& initMutant);

 void print_spP(const spParamsP& spP);

@@ -166,6 +169,37 @@ void updateRatesMcFarlandLog(std::vector<spParamsP>& popParams,

 				    const double& K,

 			     const double& totPopSize);

+void updateRatesFDFMcFarlandLog(std::vector<spParamsP>& popParams,

+  const std::vector<Genotype>& Genotypes,

+  const fitnessEffectsAll& fitnessEffects,

+  double& adjust_fitness_MF,

+  const double& K,

+  const double& totPopSize,

+  const double& currentTime);

+

+void updateRatesMcFarlandLog_D(std::vector<spParamsP>& popParams,

+				    double& adjust_fitness_MF,

+				    const double& K,

+			     const double& totPopSize);

+

+void updateRatesFDFMcFarlandLog_D(std::vector<spParamsP>& popParams,

+  const std::vector<Genotype>& Genotypes,

+  const fitnessEffectsAll& fitnessEffects,

+  double& adjust_fitness_MF,

+  const double& K,

+  const double& totPopSize,

+  const double& currentTime);

+

+

+void updateRatesFDFExp(std::vector<spParamsP>& popParams,

+  const std::vector<Genotype>& Genotypes,

+  const fitnessEffectsAll& fitnessEffects,

+  const double& currentTime);

+

+void updateRatesFDFBozic(std::vector<spParamsP>& popParams,

+  const std::vector<Genotype>& Genotypes,

+  const fitnessEffectsAll& fitnessEffects,

+  const double& currentTime);

 void updateRatesMcFarland0(std::vector<spParamsP>& popParams,

 				  double& adjust_fitness_MF,

@@ -188,4 +222,3 @@ void detect_ti_duplicates(const std::multimap<double, int>& m,

 			  const int spcies);

 #endif

-

@@ -186,5 +186,6 @@ void updateRatesBeeren(std::vector<spParamsP>& popParams,

 void detect_ti_duplicates(const std::multimap<double, int>& m,

 			  const double ti,

 			  const int spcies);

+

 #endif

@@ -183,6 +183,8 @@ void updateRatesBeeren(std::vector<spParamsP>& popParams,

 			      const int& mutationPropGrowth,

 		       const double& mu);

-

+void detect_ti_duplicates(const std::multimap<double, int>& m,

+			  const double ti,

+			  const int spcies);

 #endif

@@ -114,33 +114,47 @@ void precissionLoss();

 void init_tmpP(spParamsP& tmpParam);

-double returnMFE(double& e1,

-		 // const double& K,

-		 const std::string& typeFitness);

-

-double returnMFE(double& e1,

-		 // const double& K,

-		 const TypeModel typeModel);

-

-void computeMcFarlandError(double& e1,

-			   double& n_0,

-			   double& n_1,

-			   double& tps_0,

-			   double& tps_1,

-			   const std::string& typeFitness,

-			   const double& totPopSize,

-			   const double& K);

-

-void computeMcFarlandError(double& e1,

-			   double& n_0,

-			   double& n_1,

-			   double& tps_0,

-			   double& tps_1,

-			   const TypeModel typeModel,

-			   const double& totPopSize,

-			   const double& K);

-

-

+// double returnMFE(double& e1,

+// 		 const std::string& typeFitness);

+

+// double returnMFE(double& e1,

+// 		 const TypeModel typeModel);

+

+double returnMFE_new(double& en1,

+		     const std::string& typeFitness);

+

+double returnMFE_new(double& en1,

+		     const TypeModel typeModel);

+

+// void computeMcFarlandError(double& e1,

+// 			   double& n_0,

+// 			   double& tps_0,

+// 			   const std::string& typeFitness,

+// 			   const double& totPopSize,

+// 			   const double& K);

+

+// void computeMcFarlandError(double& e1,

+// 			   double& n_0,

+// 			   double& tps_0,

+// 			   const TypeModel typeModel,

+// 			   const double& totPopSize,

+// 			   const double& K);

+

+void computeMcFarlandError_new(double& en1,

+			       double& en1sc,

+			       double& totPopSize_previous,

+			       double& DA_previous,

+			       const TypeModel typeModel,

+			       const double& totPopSize,

+			       const double& K);

+

+void computeMcFarlandError_new(double& en1,

+			       double& en1sc,

+			       double& totPopSize_previous,

+			       double& DA_previous,

+			       const std::string& typeFitness,

+			       const double& totPopSize,

+			       const double& K);

 void updateRatesMcFarland(std::vector<spParamsP>& popParams,

 				 double& adjust_fitness_MF,

@@ -19,29 +19,30 @@

 #define _BNB_COMMON_H_

 #include<Rcpp.h>

-#include"debug_common.h"

-#include "new_restrict.h" // for the TypeModel enum

-

-// Simple custom exception for exceptions that lead to re-runs.

-class rerunExcept: public std::runtime_error {

-public:

-  rerunExcept(const std::string &s) :

-    std::runtime_error(s) {}

-};

-

-

-struct spParamsP {

-  double popSize;

-  double birth;

-  double death;

-  double W;

-  double R;

-  double mutation; 

-  double timeLastUpdate;

-  std::multimap<double, int>::iterator pv;

-  double absfitness; //convenient for Beerenwinkel

-  int numMutablePos; //for mutator if need update of mutation

-};

+#include "common_classes.h"

+#include "debug_common.h"

+// #include "new_restrict.h" // for the TypeModel enum

+

+// // Simple custom exception for exceptions that lead to re-runs.

+// class rerunExcept: public std::runtime_error {

+// public:

+//   rerunExcept(const std::string &s) :

+//     std::runtime_error(s) {}

+// };

+

+

+// struct spParamsP {

+//   double popSize;

+//   double birth;

+//   double death;

+//   double W;

+//   double R;

+//   double mutation; 

+//   double timeLastUpdate;

+//   std::multimap<double, int>::iterator pv;

+//   double absfitness; //convenient for Beerenwinkel

+//   int numMutablePos; //for mutator if need update of mutation

+// };

 inline void W_f_st(spParamsP& spP){

@@ -58,7 +59,14 @@ inline void R_f_st(spParamsP& spP) {

 inline double pE_f_st(double& pM, const spParamsP& spP){

   double pE = (spP.death * (1.0 - pM ) )/(spP.W - spP.death - spP.birth * pM );

   if( !std::isfinite(pE) ) {

-    throw std::range_error("pE.f: pE not finite");

+    DP2(spP.death);  DP2(spP.birth); DP2(pM); DP2(spP.W);

+    DP2(spP.mutation);

+    std::string error_message = R"(pE.f: pE not finite.

+      This is expected to happen when mutationPropGrowth = TRUE 

+      and you have have an initMutant with death >> birth,

+      as that inevitably leads to net birth rate of 0

+      and mutation rate of 0)";

+    throw std::range_error(error_message);

   }

   return pE;

 }

@@ -97,7 +105,8 @@ double ti_nextTime_tmax_2_st(const spParamsP& spP,