//     Copyright 2013, 2014, 2015 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
 
 // // 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;
 }
 
 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 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(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,
 				 const double& K,
 			  const double& totPopSize);
 
 void updateRatesMcFarlandLog(std::vector<spParamsP>& popParams,
 				    double& adjust_fitness_MF,
 				    const double& K,
 			     const double& totPopSize);
 
 
 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 detect_ti_duplicates(const std::multimap<double, int>& m,
 			  const double ti,
 			  const int spcies);

 #endif