//     Copyright 2013, 2014, 2015, 2016 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/>.
 
 

 // #include "randutils.h" //Nope, until we have gcc-4.8 in Win; full C++11

 #include "debug_common.h"

 #include "common_classes.h"
 #include "new_restrict.h"

 #include <Rcpp.h>
 #include <iomanip> 
 #include <algorithm>
 #include <random>

 #include <string>

 
 
 using namespace Rcpp;
 using std::vector;
 using std::back_inserter;
 
 

 std::string concatIntsString(const std::vector<int>& ints,
 			     const std::string sep = ", ") {
   std::string strout;
   std::string comma = "";
   for(auto const &g : ints) {
     strout += (comma + std::to_string(g));
     comma = sep;
   }
   return strout;
 }
 
 

 double prodFitness(const std::vector<double>& s) {

   return accumulate(s.begin(), s.end(), 1.0,
 		    [](double x, double y) {return (x * std::max(0.0, (1 + y)));});
 }
 

 // // This is a better idea? If yes, change code in nr_fitness so that
 // // birth 0 is not extinction.
 // double prodFitnessNegInf(std::vector<double> s) {
 //   double f = 1.0;
 //   for(auto y : s) {
 //     if( y == -std::numeric_limits<double>::infinity()) {
 //       return -std::numeric_limits<double>::infinity();
 //     } else {
 //       f *= std::max(0.0, (1 + y));
 //     }
 //   }
 //   return f;
 // }

 double prodDeathFitness(const std::vector<double>& s) {

   return accumulate(s.begin(), s.end(), 1.0,
 		    [](double x, double y) {return (x * std::max(0.0, (1 - y)));});
 }
 

 double prodMuts(const std::vector<double>& s) {
   // From prodFitness
   // return accumulate(s.begin(), s.end(), 1.0,
   // 		    [](double x, double y) {return (x * y);});
   return accumulate(s.begin(), s.end(), 1.0,
 		    std::multiplies<double>());
 }

 double set_cPDetect(const double n2, const double p2,
 		    const double PDBaseline) {
   return (-log(1.0 - p2)/(n2 - PDBaseline));
 }
 
 // Next two are identical, except for scaling the k. Use the simplest.
 double probDetectSize(const double n, const double cPDetect,
 		      const double PDBaseline) {
   if(n <= PDBaseline) {
     return 0;
   } else {
     return (1 - exp( -cPDetect * (n - PDBaseline)));
   }
 }
 
 // double prob_exit_ratio(const double n, const double k, const double baseline) {
 //   if(n <= baseline) {
 //     return 0;
 //   } else {
 //     return (1 - exp( -c * ( (n - baseline)/baseline)));
 //   }
 // }
 
 
 bool detectedSizeP(const double n, const double cPDetect,
 			const double PDBaseline, std::mt19937& ran_gen) {
   if(cPDetect < 0) {
     // As we OR, return false if this condition does not apply
     return false;
   } else {
     std::uniform_real_distribution<double> runif(0.0, 1.0);
     double prob = probDetectSize(n, cPDetect, PDBaseline);

     if(prob <= 0.0) return false;

     if(runif(ran_gen) <= prob) {
       return true;
     } else {
       return false;
     }
   }
 }

 
 bool operator==(const Genotype& lhs, const Genotype& rhs) {
   return (lhs.orderEff == rhs.orderEff) &&
     (lhs.epistRtEff == rhs.epistRtEff) &&

     (lhs.rest == rhs.rest) &&
     (lhs.flGenes == rhs.flGenes);

 }
 

 // Added for completeness, but not used now
 // bool operator<(const Genotype& lhs, const Genotype& rhs) {
 //   std::vector<int> lh = genotypeSingleVector(lhs);
 //   std::vector<int> rh = genotypeSingleVector(rhs);
 //   if( lh.size() < rh.size() ) return true;
 //   else if ( lh.size() > rh.size() ) return false;
 //   else {
 //     for(size_t i = 0; i != lh.size(); ++i) {
 //       if( lh[i] < rh[i] ) return true;
 //     }
 //     return false;
 //   }
 // }
 

 TypeModel stringToModel(const std::string& mod) {
   if(mod == "exp")
     return TypeModel::exp;
   else if(mod == "bozic1")
     return TypeModel::bozic1;
   else if(mod == "mcfarlandlog")
     return TypeModel::mcfarlandlog;

   // else if(mod == "mcfarland")
   //   return TypeModel::mcfarland;
   // else if(mod == "beerenwinkel")
   //   return TypeModel::beerenwinkel;
   // else if(mod == "mcfarland0")
   //   return TypeModel::mcfarland0;
   // else if(mod == "bozic2")
   //   return TypeModel::bozic2;

   else 
     throw std::out_of_range("Not a valid TypeModel");
 }

 
 
 Dependency stringToDep(const std::string& dep) {
   if(dep == "monotone") // AND, CBN, CMPN
     return Dependency::monotone;
   else if(dep == "semimonotone") // OR, SMN, DMPN
     return Dependency::semimonotone;
   else if(dep == "xmpn") // XOR, XMPN
     return Dependency::xmpn;
   else if(dep == "--") // for root, for example
     return Dependency::single;
   else 
     throw std::out_of_range("Not a valid typeDep");
   // We never create the NA from entry data. NA is reserved for Root.
 }
 

 // std::string depToString(const Dependency dep) {
 //   switch(dep) {
 //   case Dependency::monotone:
 //     return "CMPN or monotone";
 //   case Dependency::semimonotone:
 //     return "DMPN or semimonotone";
 //   case Dependency::xmpn:
 //     return "XMPN (XOR)";
 //   case Dependency::single:
 //     return "--";
 //   default:
 //     throw std::out_of_range("Not a valid dependency");
 //   }
 // }

 
 
 void print_Genotype(const Genotype& ge) {
   Rcpp::Rcout << "\n Printing Genotype";
   Rcpp::Rcout << "\n\t\t order effects genes:";
   for(auto const &oo : ge.orderEff) Rcpp::Rcout << " " << oo;
   Rcpp::Rcout << "\n\t\t epistasis and restriction effects genes:";
   for(auto const &oo : ge.epistRtEff) Rcpp::Rcout << " " << oo;
   Rcpp::Rcout << "\n\t\t non interaction genes :";
   for(auto const &oo : ge.rest) Rcpp::Rcout << " " << oo;

   Rcpp::Rcout << "\n\t\t fitness landscape genes :";
   for(auto const &oo : ge.flGenes) Rcpp::Rcout << " " << oo;

   Rcpp::Rcout << std::endl;
 }
 
 vector<int> genotypeSingleVector(const Genotype& ge) {
   // orderEff in the order they occur. All others are sorted.
   std::vector<int> allgG;
   allgG.insert(allgG.end(), ge.orderEff.begin(), ge.orderEff.end());
   allgG.insert(allgG.end(), ge.epistRtEff.begin(), ge.epistRtEff.end());
   allgG.insert(allgG.end(), ge.rest.begin(), ge.rest.end());

   allgG.insert(allgG.end(), ge.flGenes.begin(), ge.flGenes.end());

   // this should not be unique'd as it aint' sorted

   return allgG;
 }
 
 
 vector<int> allGenesinFitness(const fitnessEffectsAll& F) {
   // Sorted
   std::vector<int> g0;
 
   if(F.Gene_Module_tabl.size()) {
     if( F.Gene_Module_tabl[0].GeneNumID != 0 )

       throw std::logic_error("\n Gene module table's first element must be 0."
 			     " This should have been caught in R.");

     // for(vector<Gene_Module_struct>::size_type i = 1;
     //     i != F.Gene_Module_tabl.size(); i++) {
     for(decltype(F.Gene_Module_tabl.size()) i = 1;
   	i != F.Gene_Module_tabl.size(); i++) {
       g0.push_back(F.Gene_Module_tabl[i].GeneNumID);
     }
   }
   // for(auto const &a : F.Gene_Module_tabl) {
   //    if(a.GeneNumID != 0) g0.push_back(a.GeneNumID);
   // }
   for(auto const &b: F.genesNoInt.NumID) {
     g0.push_back(b);
   }

   // sort(g0.begin(), g0.end());
   for(auto const &b: F.fitnessLandscape.NumID) {
     g0.push_back(b);
   }

   sort(g0.begin(), g0.end());

   // Can we assume the fitness IDs go from 0 to n? Nope: because of
   // muEF. But we assume in several places that there are no repeated
   // elements in the output from this function.
 
   // FIXME we verify there are no repeated elements. That is to strongly
   // check our assumptions are right. Alternatively, return the "uniqued"
   // vector and do not check anything.
   std::vector<int> g0_cp(g0);
   g0.erase( unique( g0.begin(), g0.end() ), g0.end() );
   if(g0.size() != g0_cp.size())

     throw std::logic_error("\n allGenesinFitness: repeated genes. "
 			   " This should have been caught in R.");

   return g0;
 }
 
 vector<int> allGenesinGenotype(const Genotype& ge){

   // Like genotypeSingleVector, but sorted

   std::vector<int> allgG;
   for(auto const &g1 : ge.orderEff)
     allgG.push_back(g1);
   for(auto const &g2 : ge.epistRtEff)
     allgG.push_back(g2);
   for(auto const &g3 : ge.rest)
     allgG.push_back(g3);

   for(auto const &g4 : ge.flGenes)
     allgG.push_back(g4);
   

   sort(allgG.begin(), allgG.end());

   // Remove duplicates see speed comparisons here:
   // http://stackoverflow.com/questions/1041620/whats-the-most-efficient-way-to-erase-duplicates-and-sort-a-vector
   // We assume here there are no duplicates. Yes, a gene can have both
   // fitness effects and order effects and be in the DAG. But it will be
   // in only one of the buckets.
   std::vector<int> g0_cp(allgG);
   allgG.erase( unique( allgG.begin(), allgG.end() ), allgG.end() );
   if(allgG.size() != g0_cp.size())

     throw std::logic_error("\n allGenesinGenotype: repeated genes."
 			   " This should have been caught in R.");

   return allgG;
 }
 
 
 // For users: if something depends on 0, that is it. No further deps.
 // And do not touch the 0 in Gene_Module_table.
 std::vector<Poset_struct> rTable_to_Poset(Rcpp::List rt) { 
 
   // The restriction table, or Poset, has a first element
   // with nothing, so that all references by mutated gene
   // are simply accessing the Poset[mutated gene] without
   // having to remember to add 1, etc.
 
   std::vector<Poset_struct> Poset;
   
   Poset.resize(rt.size() + 1);
   Poset[0].child = "0"; //should this be Root?? I don't think so.
   Poset[0].childNumID = 0;
   Poset[0].typeDep = Dependency::NA;
   Poset[0].s = std::numeric_limits<double>::quiet_NaN();
   Poset[0].sh = std::numeric_limits<double>::quiet_NaN();
   Poset[0].parents.resize(0);
   Poset[0].parentsNumID.resize(0);
 
   Rcpp::List rt_element;
   // std::string tmpname;
   Rcpp::IntegerVector parentsid;
   Rcpp::CharacterVector parents;
 
   for(int i = 1; i != (rt.size() + 1); ++i) {
     rt_element = rt[i - 1];
     Poset[i].child = Rcpp::as<std::string>(rt_element["child"]);
     Poset[i].childNumID = as<int>(rt_element["childNumID"]);
     Poset[i].typeDep = stringToDep(as<std::string>(rt_element["typeDep"]));
     Poset[i].s = as<double>(rt_element["s"]);
     Poset[i].sh = as<double>(rt_element["sh"]);
 
     // if(i != Poset[i].childNumID) {
     // Nope: this assumes we only deal with posets.
     //   // Rcpp::Rcout << "\n childNumID, original = " << as<int>(rt_element["childNumID"]);
     //   // Rcpp::Rcout << "\n childNumID, Poset = " << Poset[i].childNumID;
     //   // Rcpp::Rcout << "\n i = " << i << std::endl;
     //   throw std::logic_error("childNumID != index");
     // }
     // Rcpp::IntegerVector parentsid = as<Rcpp::IntegerVector>(rt_element["parentsNumID"]);
     // Rcpp::CharacterVector parents = as<Rcpp::CharacterVector>(rt_element["parents"]);
 
     parentsid = as<Rcpp::IntegerVector>(rt_element["parentsNumID"]);
     parents = as<Rcpp::CharacterVector>(rt_element["parents"]);
 
     if( parentsid.size() != parents.size() ) {

       throw std::logic_error("parents size != parentsNumID size. Bug in R code.");

     }
 
     for(int j = 0; j != parentsid.size(); ++j) {
       Poset[i].parentsNumID.push_back(parentsid[j]);
       Poset[i].parents.push_back( (Rcpp::as< std::string >(parents[j])) );
       // tmpname = Rcpp::as< std::string >(parents[j]);
       // Poset[i].parents.push_back(tmpname);
     }
 
     // Should not be needed if R always does what is should. Disable later?
     if(! is_sorted(Poset[i].parentsNumID.begin(), Poset[i].parentsNumID.end()) )

       throw std::logic_error("ParentsNumID not sorted. Bug in R code.");

     if(std::isinf(Poset[i].s))

       Rcpp::Rcout << "WARNING: at least one s is infinite" 
 		  << std::endl;

     if(std::isinf(Poset[i].sh) && (Poset[i].sh > 0))

       Rcpp::Rcout << "WARNING: at least one sh is positive infinite" 
 		  << std::endl;
   }
   return Poset;
 }
 
 
 std::vector<Gene_Module_struct> R_GeneModuleToGeneModule(Rcpp::List rGM) {
 
   std::vector<Gene_Module_struct> geneModule;
 
   Rcpp::IntegerVector GeneNumID = rGM["GeneNumID"];
   Rcpp::IntegerVector ModuleNumID = rGM["ModuleNumID"];
   Rcpp::CharacterVector GeneName = rGM["Gene"];
   Rcpp::CharacterVector ModuleName = rGM["Module"];
   // geneModule.resize(GeneNumID.size());
   geneModule.resize(GeneNumID.size()); // remove later
 
   for(size_t i = 0; i != geneModule.size(); ++i) {
     if( static_cast<int>(i) != GeneNumID[i])

       throw std::logic_error(" i != GeneNumID. Bug in R code.");

     // geneModule[i].GeneNumID = GeneNumID[i];
     // geneModule[i].ModuleNumID = ModuleNumID[i];
     // remove these later?
     geneModule[i].GeneNumID = GeneNumID[i];
     geneModule[i].ModuleNumID = ModuleNumID[i];
     geneModule[i].GeneName = GeneName[i];
     geneModule[i].ModuleName = ModuleName[i];
   }
  
   return geneModule;
 }
 
 
 std::vector<int> GeneToModule(const std::vector<int>& Drv,
 			     const 
 			      std::vector<Gene_Module_struct>& Gene_Module_tabl,
 			      const bool sortout, const bool uniqueout) {
   
   std::vector<int>  mutatedModules;
   
   for(auto it = Drv.begin(); it != Drv.end(); ++it) {
     mutatedModules.push_back(Gene_Module_tabl[(*it)].ModuleNumID);
   }
   // sortout and uniqueout returns a single element of each. uniqueout only removes
   // successive duplicates. sortout without unique is just useful for knowing
   // what happens for stats, etc. Neither sortout nor uniqueout for keeping
   // track of order of module events.
   if(sortout) {
     sort( mutatedModules.begin(), mutatedModules.end() );
   }
   if(uniqueout) {
     mutatedModules.erase( unique( mutatedModules.begin(), 
 				  mutatedModules.end() ), 
 			  mutatedModules.end() );
   }
   return mutatedModules;
 }
 
 

 // fitnessLandscape_struct convertFitnessLandscape(Rcpp::List flg,
 // 						Rcpp::DataFrame fl_df) {

 fitnessLandscape_struct convertFitnessLandscape(Rcpp::List flg,

 						Rcpp::List fl_df) {
 

   fitnessLandscape_struct flS;
   flS.names = Rcpp::as<std::vector<std::string> >(flg["Gene"]);
   flS.NumID = Rcpp::as<std::vector<int> >(flg["GeneNumID"]);
 
   std::vector<std::string> genotNames =
     Rcpp::as<std::vector<std::string> >(fl_df["Genotype"]);
   // Rcpp::CharacterVector genotNames = fl_df["Genotype"];
   Rcpp::NumericVector fitness = fl_df["Fitness"];
   
   // Fill up the map genotypes(as string) to fitness

   //for(size_t i = 0; i != fl_df.nrows(); ++i) {
   for(size_t i = 0; i != genotNames.size(); ++i) {

     flS.flmap.insert({genotNames[i], fitness[i]});
   }

   return flS;
 }

 
 genesWithoutInt convertNoInts(Rcpp::List nI) {
 
   genesWithoutInt genesNoInt;
   // FIXME: I think I want to force Gene in long.geneNoInt to be a char.vector
   // to avoid this transformation.
   // Rcpp::CharacterVector names = nI["Gene"];
   // Rcpp::IntegerVector id = nI["GeneNumID"];
   // Rcpp::NumericVector s1 = nI["s"];
   
   // genesNoInt.names = Rcpp::as<std::vector<std::string> >(names);
   genesNoInt.names = Rcpp::as<std::vector<std::string> >(nI["Gene"]);
   genesNoInt.NumID = Rcpp::as<std::vector<int> >(nI["GeneNumID"]);
   genesNoInt.s = Rcpp::as<std::vector<double> >(nI["s"]);
   genesNoInt.shift = genesNoInt.NumID[0]; // FIXME: we assume mutations always
 				      // indexed 1 to something. Not 0 to
 				      // something.
   return genesNoInt;
 }
 
 
 std::vector<epistasis> convertEpiOrderEff(Rcpp::List ep) {
 
   std::vector<epistasis> Epistasis;
   
   Rcpp::List element;
   // For epistasis, the numID must be sorted, but never with order effects.
   // Things come sorted (or not) from R.
   Epistasis.resize(ep.size());
   for(int i = 0; i != ep.size(); ++i) {
     element = ep[i];
     Epistasis[i].NumID = Rcpp::as<std::vector<int> >(element["NumID"]);
     Epistasis[i].names = Rcpp::as<std::vector<std::string> >(element["ids"]);
     Epistasis[i].s = as<double>(element["s"]);
   }
   return Epistasis;
 }
 

 std::vector<int> sortedAllOrder(const std::vector<epistasis>& E) {

   
   std::vector<int> allG;
   for(auto const &ec : E) {
     for(auto const &g : ec.NumID) {
       allG.push_back(g);
     }
   }
   sort(allG.begin(), allG.end());
   allG.erase( unique( allG.begin(), allG.end()),
 		      allG.end());
   return allG;
 }
 

 std::vector<int> sortedAllPoset(const std::vector<Poset_struct>& Poset) {

   // Yes, this could be done inside rTable_to_Poset but this is cleaner
   // and will only add very little time. 
   std::vector<int> allG;
   for(auto const &p : Poset) {
     allG.push_back(p.childNumID);
   }
   sort(allG.begin(), allG.end());
   allG.erase( unique( allG.begin(), allG.end()),
 		      allG.end());
   return allG; 
 }
 
 fitnessEffectsAll convertFitnessEffects(Rcpp::List rFE) {
   // Yes, some of the things below are data.frames in R, but for
   // us that is used just as a list.
 
   fitnessEffectsAll fe;
   
   Rcpp::List rrt = rFE["long.rt"];
   Rcpp::List re = rFE["long.epistasis"];
   Rcpp::List ro = rFE["long.orderEffects"];
   Rcpp::List rgi = rFE["long.geneNoInt"];
   Rcpp::List rgm = rFE["geneModule"];

   bool rone = as<bool>(rFE["gMOneToOne"]);

   Rcpp::IntegerVector drv = rFE["drv"];
 

   Rcpp::List flg = rFE["fitnessLandscape_gene_id"];

   // clang does not like this
   // Rcpp::DataFrame fl_df = rFE["fitnessLandscape_df"];
   Rcpp::List fl_df = rFE["fitnessLandscape_df"];
   

   
  
 
   // In the future, if we want noInt and fitnessLandscape, all
   // we need is use the fitness landscape with an index smaller than those
   // of noInt. So we can use noInt with shift being those in fitnessLandscape.
   // BEWARE: will need to modify also createNewGenotype.
   

   // if(fl_df.nrows()) {
   if(fl_df.size()) {

     fe.fitnessLandscape = convertFitnessLandscape(flg, fl_df);
   }
   

   if(rrt.size()) {
     fe.Poset = rTable_to_Poset(rrt);
   } 
   if(re.size()) {
     fe.Epistasis = convertEpiOrderEff(re);
   } 
   if(ro.size()) {
     fe.orderE = convertEpiOrderEff(ro);
   } 
   if(rgi.size()) {
     fe.genesNoInt = convertNoInts(rgi);
   } else {
     fe.genesNoInt.shift = -9L;
   }

   // If this is null, use the nullFitnessEffects function; never
   // end up here.

   
   //   if( (rrt.size() + re.size() + ro.size() + rgi.size() + fl_df.nrows()) == 0) {
   if( (rrt.size() + re.size() + ro.size() + rgi.size() + fl_df.size()) == 0) {

       throw std::logic_error("\n Nothing inside this fitnessEffects; why are you here?"
 			     "  Bug in R code.");

   }
   

   // At least for now, if fitness landscape nothing else allowed

   // if(fl_df.nrows() && ((rrt.size() + re.size() + ro.size() + rgi.size()) > 0)) {
   if(fl_df.size() && ((rrt.size() + re.size() + ro.size() + rgi.size()) > 0)) {

     throw std::logic_error("\n Fitness landscape specification."
 			   " There should be no other terms. "
 			   " Bug in R code");
   }
 
   // This is silly
   // if(fl_df.nrows() && (rgm.size() > 4) ) {
   //   throw std::logic_error("\n Fitness landscape specification."
   // 			   " Cannot use modules. "
   // 			   " Bug in R code");
   // }
   

   fe.Gene_Module_tabl = R_GeneModuleToGeneModule(rgm);
   fe.allOrderG = sortedAllOrder(fe.orderE);
   fe.allPosetG = sortedAllPoset(fe.Poset);
   fe.gMOneToOne = rone;
   fe.allGenes = allGenesinFitness(fe);

   fe.genomeSize =  fe.Gene_Module_tabl.size() - 1 + fe.genesNoInt.s.size() +
     fe.fitnessLandscape.NumID.size();

   fe.drv = as<std::vector<int> > (drv);
   sort(fe.drv.begin(), fe.drv.end()); //should not be needed, but just in case
   // cannot trust R gives it sorted
   // check_disable_later
   if(fe.genomeSize != static_cast<int>(fe.allGenes.size())) {

     throw std::logic_error("\n genomeSize != allGenes.size(). Bug in R code.");

   }

   // At least for now
   if(fe.fitnessLandscape.NumID.size() > 0) {
     if(fe.genomeSize != static_cast<int>(fe.fitnessLandscape.NumID.size())) {
       throw std::logic_error("\n genomeSize != genes in fitness landscape."
 			     "Bug in R code.");
     }
   }
   

   return fe;
 }
 

 // Before making allGenesinGenotype return a unique vector: we do a
 // set_difference below. If we look at the help
 // (http://en.cppreference.com/w/cpp/algorithm/set_difference) if we had
 // more repetitions of an element in allGenes than in sortedparent we
 // could have a problem. But if you look at function "allgenesinFitness",
 // which is the one used to give the allgenes vector, you will see that
 // that one returns only one entry per gene, as it parses the geneModule
 // structure. So even if allGenesinGenotype returns multiple entries
 // (which they don't), there will be no bugs as the maximum number of
 // entries in the output of setdiff will be 0 or 1 as m is 1. But
 // allGenesinGenotype cannot return more than one element as can be seen
 // in createNewGenotype: an element only ends in the order component if it
 // is not in the epistasis component.  So there are no repeated elements
 // in allGenes or in sortedparent below. Also, beware that does not break
 // correct fitness evaluation of fitnessEffects where the same gene is in
 // epistasis and order, as can be seen in the tests and because of how we
 // evaluate fitness, where genes in a genotype in the orderEffects bucket
 // are placed also in the epist for fitness eval. See evalGenotypeFitness
 // and createNewGenotype.

 void obtainMutations(const Genotype& parent,
 		     const fitnessEffectsAll& fe,

 		     int& numMutablePosParent, 

 		     std::vector<int>& newMutations,

 		     //randutils::mt19937_rng& ran_gen

 		     std::mt19937& ran_gen,
 		     std::vector<double> mu) {

   //Ugly: we return the mutations AND the numMutablePosParent This is
   // almost ready to accept multiple mutations. And it returns a vector,
   // newMutations.
   std::vector<int> sortedparent = allGenesinGenotype(parent);
   std::vector<int> nonmutated;
   set_difference(fe.allGenes.begin(), fe.allGenes.end(),
 		 sortedparent.begin(), sortedparent.end(),
 		 back_inserter(nonmutated));

   // numMutablePos is used not only for mutation but also to decide about
   // the dummy or null mutation case.
   numMutablePosParent = nonmutated.size();
   if(nonmutated.size() < 1)

     throw std::out_of_range("Trying to obtain a mutation when nonmutated.size is 0."
 			    " Bug in R code; let us know.");

   if(mu.size() == 1) { // common mutation rate
     // FIXME:chromothr would not use this, or this is the limit case with a
   // single mutant
     std::uniform_int_distribution<int> rpos(0, nonmutated.size() - 1);
     newMutations.push_back(nonmutated[rpos(ran_gen)]);
   } else { // per-gene mutation rate.
     // Remember that mutations always indexed from 1, not from 0.
     // We take an element from a discrete distribution, with probabilities
     // proportional to the rates.
     // FIXME:varmutrate give a warning if the only mu is for mu = 0.
     std::vector<double> mu_nm;
     for(auto const &nm : nonmutated) mu_nm.push_back(mu[nm - 1]);
     std::discrete_distribution<int> rpos(mu_nm.begin(), mu_nm.end());
     newMutations.push_back(nonmutated[rpos(ran_gen)]);
   }

   // randutils
   // // Yes, the next will work, but pick is simpler!
   // // size_t rpos = ran_gen.uniform(static_cast<size_t>(0), nonmutated.size() - 1);
   // //  newMutations.push_back(nonmutated[rpos]);
   // int posmutated = ran_gen.pick(nonmutated);
   // newMutations.push_back(posmutated);

 }
 
 
 // std::vector<int> genesInOrderModules(const fitnessEffectsAll& fe) {
 //   vector<int> genes;
 //   if(fe.gMOneToOne)
 //     genes = fe.alOrderG;
 //   else {
 //     for(auto const &m : fe.allOrderG) {
 //       genes.push_back()
 //     }
 
 //   }
 //   return genes;  
 // }
 
 
 fitness_as_genes fitnessAsGenes(const fitnessEffectsAll& fe) {

   // Give the fitnessEffects in terms of genes, not modules.

   
   // Extract the noInt. Then those in order effects by creating a multimap
   // to go from map to genes. Then all remaining genes are those only in
   // poset. By set_difference.

   fitness_as_genes fg = zero_fitness_as_genes();

   // fitness_as_genes fg;
   fg.flGenes = fe.fitnessLandscape.NumID;
   if(fg.flGenes.size()) {
     return fg;
   }

   fg.noInt = fe.genesNoInt.NumID;
   // //zz: debugging
   // for(auto const &o : fe.genesNoInt.NumID) {
   //   DP2(o);
   // }
   // for(auto const &o : fe.genesNoInt.names) {
   //   DP2(o);
   // }
   // // 

   std::multimap<int, int> MG;
   for( auto const &mt : fe.Gene_Module_tabl) {
     MG.insert({mt.ModuleNumID, mt.GeneNumID});
   }
   for (auto const &o : fe.allOrderG) {
     for(auto pos = MG.lower_bound(o); pos != MG.upper_bound(o); ++pos) 
       fg.orderG.push_back(pos->second);
   }
   sort(fg.orderG.begin(), fg.orderG.end());
 
   std::vector<int> tmpv = fg.orderG;
   tmpv.insert(tmpv.end(),fg.noInt.begin(), fg.noInt.end());
   sort(tmpv.begin(), tmpv.end()); // should not be needed
   
   set_difference(fe.allGenes.begin(), fe.allGenes.end(),
 		 tmpv.begin(), tmpv.end(),
 		 back_inserter(fg.posetEpistG));
   // fg.posetEpistG.sort(fg.posetEpistG.begin(),
   // 		      fg.posetEpistG.end());

 
   // // //zz: debugging
   // DP1("order");
   // for(auto const &o : fg.orderG) {
   //   DP2(o);
   // }
   // DP1("posetEpist");
   // for(auto const &o : fg.posetEpistG) {
   //   DP2(o);
   // }
   // DP1("noint");
   //   for(auto const &o : fg.noInt) {
   //   DP2(o);
   // }
   // DP1("fl");
   // for(auto const &o : fg.flGenes) {
   //   DP2(o);
   // }
 
   

   return fg;
 }
 
 
 std::map<int, std::string> mapGenesIntToNames(const fitnessEffectsAll& fe) {
   // This is a convenience, used in the creation of output.
   // Sure, we could do this when reading the data in.
   // The noInt in convertNoInts.
   std::map<int, std::string> gg;
 
   for(auto const &mt : fe.Gene_Module_tabl) {
     gg.insert({mt.GeneNumID, mt.GeneName});
   }
   // this is pedantic, as what is the size_type of NumID and of names? 
   // for(decltype(fe.genesNoInt.s.size()) i = 0;
   //     i != fe.genesNoInt.s.size(); ++i)
 
   for(size_t i = 0;
       i != fe.genesNoInt.NumID.size(); ++i){
     gg.insert({fe.genesNoInt.NumID[i], fe.genesNoInt.names[i]});
   }

 
   // zz
   for(size_t i = 0;
       i != fe.fitnessLandscape.NumID.size(); ++i){
     gg.insert({fe.fitnessLandscape.NumID[i], fe.fitnessLandscape.names[i]});
   }
 
   

   return gg;
 }
 
 // It is simple to write specialized functions for when
 // there are no restrictions or no order effects , etc.
 Genotype createNewGenotype(const Genotype& parent,
 			   const std::vector<int>& mutations,
 			   const fitnessEffectsAll& fe,

 			   std::mt19937& ran_gen,

 			   //randutils::mt19937_rng& ran_gen

 			   bool random = true) {
   // random: if multiple mutations, randomly shuffle the ordered ones?
   // This is the way to go if chromothripsis, but not if we give an
   // initial mutant
   

   Genotype newGenot = parent;
   std::vector<int> tempOrder; // holder for multiple muts if order.
   bool sort_rest = false;
   bool sort_epist = false;

   bool sort_flgenes = false;

 
   // FIXME: we need to create the mutations!
 
   // Order of ifs: I suspect order effects rare. No idea about
   // non-interaction genes, but if common the action is simple.

 
   // A gene that is involved both in order effects and epistasis, only
   // ends up in the orderEff container, for concision (even if a gene can
   // be involved in both orderEff and epistasis and rT). But in the
   // genotype evaluation, in evalGenotypeFitness, notice that we create
   // the vector of genes to be checked against epistais and order effects
   // using also those from orderEff:
   //   std::vector<int> mutG (ge.epistRtEff);
   //   mutG.insert( mutG.end(), ge.orderEff.begin(), ge.orderEff.end());

   for(auto const &g : mutations) {

     // If we are dealing with a fitness landscape, that is as far as we go here
     // at least for now. No other genes affect fitness.
     // But this can be easily fixed in the future; like this?
     // if(g <= (fe.fitnessLandscape.NumID.size() + 1)) {
     // and restructure the else logic for the noInt
     if(fe.fitnessLandscape.NumID.size()) {
       newGenot.flGenes.push_back(g);
       sort_flgenes = true;
     } else {
       if( (fe.genesNoInt.shift < 0) || (g < fe.genesNoInt.shift) ) { // Gene with int
 	// We can be dealing with modules
 	int m; 
 	if(fe.gMOneToOne) {
 	  m = g; 
 	} else {
 	  m = fe.Gene_Module_tabl[g].ModuleNumID;
 	}
 	if( !binary_search(fe.allOrderG.begin(), fe.allOrderG.end(), m) ) {
 	  newGenot.epistRtEff.push_back(g);
 	  sort_epist = true;
 	} else {
 	  tempOrder.push_back(g);
 	}

       } else {

 	// No interaction genes so no module stuff
 	newGenot.rest.push_back(g);
 	sort_rest = true;

       }

     }    
   }
   

   // If there is order but multiple simultaneous mutations
   // (chromothripsis), we randomly insert them

 
   // FIXME: initMutant cannot use this!! we give the order!!!
   if( (tempOrder.size() > 1) && random)

     shuffle(tempOrder.begin(), tempOrder.end(), ran_gen);

   // The new randutils engine:
   // if(tempOrder.size() > 1)
   //   ran_gen.shuffle(tempOrder.begin(), tempOrder.end());
 
   

   for(auto const &g : tempOrder)
     newGenot.orderEff.push_back(g);
 
   // Sorting done at end, in case multiple mutations
   if(sort_rest)
     sort(newGenot.rest.begin(), newGenot.rest.end());
   if(sort_epist)
     sort(newGenot.epistRtEff.begin(), newGenot.epistRtEff.end());

   if(sort_flgenes)
     sort(newGenot.flGenes.begin(), newGenot.flGenes.end());

   return newGenot;
 }
 
 // FIXME: Prepare specialized functions: 
 // Specialized functions:
 // Never interactions: push into rest and sort. Identify by shift == 1.
 // Never no interactions: remove the if. shift == -9.
 
 
 
 

 // A paranoid check in case R code has bugs.

 void breakingGeneDiff(const vector<int>& genotype,
 		      const vector<int>& fitness) {

   std::vector<int> diffg;

   set_difference(genotype.begin(), genotype.end(),
 		 fitness.begin(), fitness.end(),
 		 back_inserter(diffg));
   if(diffg.size()) {
     Rcpp::Rcout << "Offending genes :";
     for(auto const &gx : diffg) {
       Rcpp::Rcout << " " << gx;
     }

     Rcpp::Rcout << "\t Genotype: ";
     for(auto const &g1 : genotype) Rcpp::Rcout << " " << g1;
     Rcpp::Rcout << "\t Fitness: ";
     for(auto const &g1 : fitness) Rcpp::Rcout << " " << g1;
 

     Rcpp::Rcout << "\n ";

     throw std::logic_error("\n At least one gene in the genotype not in fitness effects."
 			   " Bug in R code.");

   }
 }
 
 void checkNoNegZeroGene(const vector<int>& ge) {
   if( ge[0] == 0 )

     throw std::logic_error("\n Genotype cannot contain 0. Bug in R code.");

   else if(ge[0] < 0)

     throw std::logic_error("\n Genotype cannot contain negative values. Bug in R code.");

 }
 
 void checkLegitGenotype(const Genotype& ge,
 			const fitnessEffectsAll& F) {

   if((ge.orderEff.size() + ge.epistRtEff.size() + ge.rest.size()) == 0) {
     // An empty genotype is always legitimate, even if silly
     return;
   }
   // vector<int> g0 = allGenesinFitness(F);
   vector<int> g0 = F.allGenes;

   vector<int> allgG = allGenesinGenotype(ge);
   checkNoNegZeroGene(allgG);
   breakingGeneDiff(allgG, g0);
 }
 
 void checkLegitGenotype(const vector<int>& ge,
 			const fitnessEffectsAll& F) {

   if (ge.size() == 0) {
     // An empty genotype is always legitimate, even if silly
     return;
   }
   // std::vector<int> g0 = allGenesinFitness(F);
   vector<int> g0 = F.allGenes;

   std::vector<int> allgG (ge);
   sort(allgG.begin(), allgG.end());
   checkNoNegZeroGene(allgG);
   breakingGeneDiff(allgG, g0);
 }
 
 
 

 Genotype convertGenotypeFromInts(const std::vector<int>& gg,
 				 const fitnessEffectsAll& fe) {

   // A genotype is of one kind or another depending on what genes are of
   // what type.
   Genotype newGenot;

   
   if(gg.size() != 0) { 
     // check_disable_later
     checkLegitGenotype(gg, fe);
 
     // Very similar to logic in createNewGenotype for placing each gene in
     // its correct place, which needs to look at module mapping.
     for(auto const &g : gg) {

       if(fe.fitnessLandscape.NumID.size()) {
 	newGenot.flGenes.push_back(g);
       } else {
 	if( (fe.genesNoInt.shift < 0) || (g < fe.genesNoInt.shift) ) { // Gene with int
 	  // We can be dealing with modules
 	  int m; 
 	  if(fe.gMOneToOne) {
 	    m = g; 
 	  } else {
 	    m = fe.Gene_Module_tabl[g].ModuleNumID;
 	  }
 	  if( !binary_search(fe.allOrderG.begin(), fe.allOrderG.end(), m) ) {
 	    newGenot.epistRtEff.push_back(g);
 	  } else {
 	    newGenot.orderEff.push_back(g);
 	  }

 	} else {

 	  // No interaction genes so no module stuff
 	  newGenot.rest.push_back(g);

 	}

       }

     }    

     sort(newGenot.flGenes.begin(), newGenot.flGenes.end());

     sort(newGenot.rest.begin(), newGenot.rest.end());
     sort(newGenot.epistRtEff.begin(), newGenot.epistRtEff.end());
   } else {
     newGenot = wtGenotype(); // be explicit!!
   }

   return newGenot;
 }
 
 

 Genotype convertGenotypeFromR(Rcpp::IntegerVector rG,
 			      const fitnessEffectsAll& fe) {
   std::vector<int> gg = Rcpp::as<std::vector<int> > (rG);
   return convertGenotypeFromInts(gg, fe);
 }
 
 
 // Previous version
 // Genotype convertGenotypeFromR(Rcpp::IntegerVector rG,
 // 			      const fitnessEffectsAll& fe) {
 //   // A genotype is of one kind or another depending on what genes are of
 //   // what type.
 
 //   std::vector<int> gg = Rcpp::as<std::vector<int> > (rG);
 //   Genotype newGenot;
 
 //   // check_disable_later
 //   checkLegitGenotype(gg, fe);
 
 
 //   // Very similar to logic in createNewGenotype for placing each gene in
 //   // its correct place, which needs to look at module mapping.
 //   for(auto const &g : gg) {
 //     if( (fe.genesNoInt.shift < 0) || (g < fe.genesNoInt.shift) ) { // Gene with int
 //       // We can be dealing with modules
 //       int m; 
 //       if(fe.gMOneToOne) {
 // 	m = g; 
 //       } else {
 // 	m = fe.Gene_Module_tabl[g].ModuleNumID;
 //       }
 //       if( !binary_search(fe.allOrderG.begin(), fe.allOrderG.end(), m) ) {
 // 	newGenot.epistRtEff.push_back(g);
 //       } else {
 // 	newGenot.orderEff.push_back(g);
 //       }
 //     } else {
 //       // No interaction genes so no module stuff
 //       newGenot.rest.push_back(g);
 //     }
 //   }    
 
 //   sort(newGenot.rest.begin(), newGenot.rest.end());
 //   sort(newGenot.epistRtEff.begin(), newGenot.epistRtEff.end());
   
 //   return newGenot;
 // }
 
 

 // Genotype convertGenotypeFromR(Rcpp::List rGE) {
 
 //   Genotype g;
 			
 //   Rcpp::IntegerVector oe = rGE["orderEffGenes"];
 //   Rcpp::IntegerVector ert = rGE["epistRTGenes"];
 //   Rcpp::IntegerVector rest = rGE["noInteractionGenes"];
 
 //   g.orderEff = Rcpp::as<std::vector<int> > (oe);
 //   g.epistRtEff = Rcpp::as<std::vector<int> > (ert);
 //   g.rest = Rcpp::as<std::vector<int> > (rest);
 //   sort(g.epistRtEff.begin(), g.epistRtEff.end());
 //   sort(g.rest.begin(), g.rest.end());
 
 //   return g;
 // }
 
 
 
 

 bool match_order_effects(const std::vector<int>& O,
 			 const std::vector<int>& G) {

   //As the name says: we check if the order effect is matched

   if(G.size() < O.size()) return false;
   

   std::vector<int>::const_iterator p;

   std::vector<size_t> vdist;
   
   auto itb = G.begin();
   
   for(auto const &o : O) {
     p = find(G.begin(), G.end(), o);
     if( p == G.end() ) {
       return false;
     } else {
       vdist.push_back(std::distance( itb, p ));
     }
   }
   // Rcpp::Rcout << "     ";
   // for(auto vv : vdist ) {
   //   Rcpp::Rcout << vv << " ";
   // }
   // Rcpp::Rcout << std:: endl;
   if( is_sorted(vdist.begin(), vdist.end()) ) {
     return true;
   } else {
     return false;
   }
 }
 
 std::vector<double> evalOrderEffects(const std::vector<int>& mutatedM,
 				     const std::vector<epistasis>& OE) {
   std::vector<double> s;
   for(auto const &o : OE) {
     if(match_order_effects(o.NumID, mutatedM))
       s.push_back(o.s);
   }
   return s;
 }
 
 

 bool match_negative_epist(const std::vector<int>& E,
 			  const std::vector<int>& G) {

   // When we have things like -1, 2 in epistasis. We need to check 2 is
   // present and 1 is not present. E is the vector of epistatic coeffs,
   // and G the sorted genotype.
 
   if(G.size() < 1) return false;
    
   for(auto const &e : E) {
     if(e < 0) {
       if(binary_search(G.begin(), G.end(), -e))
 	return false;
     } else {
       if(!binary_search(G.begin(), G.end(), e))
 	return false;
     }
   }
   return true;
 }
 
 
 std::vector<double> evalEpistasis(const std::vector<int>& mutatedModules,
 				  const std::vector<epistasis>& Epistasis) {
   std::vector<double> s;
 
   // check_disable_later
   if(! is_sorted(mutatedModules.begin(), mutatedModules.end()))

     throw std::logic_error("mutatedModules not sorted in evalEpistasis."
 			   " Bug in R code.");

   
   for(auto const &p : Epistasis ) {
     if(p.NumID[0] > 0 ) {
       if(includes(mutatedModules.begin(), mutatedModules.end(),
 		   p.NumID.begin(), p.NumID.end()))
 	s.push_back(p.s);
     } else {
       if(match_negative_epist(p.NumID, mutatedModules))
 	s.push_back(p.s);
     }
   }
   // An alternative, but more confusing way
   // for(auto const &p : Epistasis ) {
   //   if(p.NumID[0] > 0 ) {
   //     if(!includes(mutatedModules.begin(), mutatedModules.end(),
   // 		   p.NumID.begin(), p.NumID.end()))
   // 	continue;
   //   } else {
   //     if(!match_negative_epist(p.NumID, mutatedModules))
   // 	continue;
   //   }
   //   s.push_back(p.s);
   // }
   return s;
 }
 
 // FIXME: can we make it faster if we know each module a single gene?
 // FIXME: if genotype is always kept sorted, with drivers first, can it be
 // faster? As well, note that number of drivers is automatically known
 // from this table of constraints.
 
 // int getPosetByChild(const int child,
 // 		    const std::vector<Poset_struct>& Poset) {
   
 // }
 
 std::vector<double> evalPosetConstraints(const std::vector<int>& mutatedModules,
 					 const std::vector<Poset_struct>& Poset,
 					 const std::vector<int>& allPosetG) {
 
   // check_disable_later
   if(! is_sorted(mutatedModules.begin(), mutatedModules.end()))

     throw std::logic_error("mutatedModules not sorted in evalPosetConstraints."
 			   " Bug in R code.");

 
   size_t numDeps;
   size_t sumDepsMet = 0;
   Dependency deptype;
   std::vector<int> parent_matches;
   std::vector<double> s;
 
   //This works reverted w.r.t. to evalOrderEffects and evalEpistasis:
   //there, I examine if the effect is present in the genotype. Here, I
   //examine if the genotype satisfies the constraints.
 
 
   // Since the genotype can contain genes not in the poset, first find
   // those that are mutated in the genotype AND are in the poset. Then
   // check if the mutated have restrictions satisfied.
   
   std::vector<int> MPintersect;
 
   
   std::set_intersection(allPosetG.begin(), allPosetG.end(),
 			mutatedModules.begin(), mutatedModules.end(),
 			std::back_inserter(MPintersect));
     
   // We know MPintersect is sorted, so we can avoid an O(n*n) loop
   size_t i = 0;
   for(auto const &m : MPintersect) {
     while ( Poset[i].childNumID != m) ++i;
     // Not to catch the twisted case of an XOR with a 0 and something else
     // as parents
     if( (Poset[i].parentsNumID[0] == 0) &&
 	(Poset[i].parentsNumID.size() == 1)) {
       s.push_back(Poset[i].s);
     } else {
       parent_matches.clear();
       std::set_intersection(mutatedModules.begin(), mutatedModules.end(),
 			    Poset[i].parentsNumID.begin(),
 			    Poset[i].parentsNumID.end(),
 			    back_inserter(parent_matches));
       sumDepsMet = parent_matches.size();
       numDeps = Poset[i].parentsNumID.size();
       deptype = Poset[i].typeDep;
 
       if( ((deptype == Dependency::semimonotone) &&
 	   (sumDepsMet)) ||
 	  ((deptype == Dependency::monotone) &&
 	   (sumDepsMet == numDeps)) ||
 	  ((deptype == Dependency::xmpn) &&
 	   (sumDepsMet == 1)) ) {
 	s.push_back(Poset[i].s);
       } else {
 	s.push_back(Poset[i].sh);
       }
     }
   }
       // for(auto const &parent : Poset[i].parentsNumID ) {
       // 	parent_module_mutated = binary_search(mutatedModules.begin(), 
       // 					      mutatedModules.end(),
       // 					      parent);
       // 	if(parent_module_mutated) {
       // 	  ++sumDepsMet;
       // 	  if( deptype == Dependency::semimonotone ) {
       // 	    known = true;
       // 	    break;
       // 	  }
       // 	  if( (deptype == Dependency::xmpn && (sumDepsMet > 1))) {
       // 	    knwon = true;
       // 	    break;
       // 	  }
       // 	}
       // }
   return s;
 }
 
 
 
 
 std::vector<double> evalGenotypeFitness(const Genotype& ge,
 					const fitnessEffectsAll& F){

   // check_disable_later
   checkLegitGenotype(ge, F);

   std::vector<double> s;

   if( (ge.orderEff.size() + ge.epistRtEff.size() +
        ge.rest.size() + ge.flGenes.size()) == 0) {

     Rcpp::warning("WARNING: you have evaluated fitness of a genotype of length zero.");

     // s.push_back(1.0); //Eh??!! 1? or 0? FIXME It should be empty! and have prodFitness
     // deal with it.

     return s;
   }
 

   // If we are dealing with a fitness landscape, that is as far as we go here
   // at least for now. No other genes affect fitness.
   // But this can be easily fixed in the future; do not return
   // s below, but keep adding, maybe the noIntGenes.
   // Recall also  prodFitness uses, well, the prod of 1 + s
   // so we want an s s.t. 1 + s = birth rate passed, 
   // which is the value in the fitness landscape as interpreted now.
   // i.e., s = birth rate - 1;
   if(F.fitnessLandscape.NumID.size()) {
     std::string gs = concatIntsString(ge.flGenes);
     if(F.fitnessLandscape.flmap.find(gs) == F.fitnessLandscape.flmap.end()) {
       s.push_back(-1.0);
     } else {
       s.push_back(F.fitnessLandscape.flmap.at(gs) - 1);
     }
     return s;
   }
    
   

   // Genes without any restriction or epistasis are just genes. No modules.
   // So simple we do it here.
   if(F.genesNoInt.shift > 0) {
     int shift = F.genesNoInt.shift;
     for(auto const & r : ge.rest ) {
       s.push_back(F.genesNoInt.s[r - shift]);
     }
   }
 
   // For the rest, there might be modules. Three different effects on
   // fitness possible: as encoded in Poset, general epistasis, order effects.
   
   // Epistatis and poset are checked against all mutations. Create single
   // sorted vector with all mutations and map to modules, if needed. Then

   // eval. 

   // Why not use a modified genotypeSingleVector without the no ints? We
   // could, but not necessary. And you can place genes in any order you
   // want, since this is not for order restrictions. That goes below.

   // Why do I put the epist first? See previous answer.

   // Why do I sort if one to one? binary searches. Not done below for order.

   std::vector<int> mutG (ge.epistRtEff);

   // A gene can be involved in epistasis and order. This gene would only
   // be in the orderEff vector, as seen in "createNewGenotype" or
   // "convertGenotypeFromInts"

   mutG.insert( mutG.end(), ge.orderEff.begin(), ge.orderEff.end());
   std::vector<int> mutatedModules;
   if(F.gMOneToOne) {
     sort(mutG.begin(), mutG.end()); 
     mutatedModules = mutG;
   } else {
     mutatedModules = GeneToModule(mutG, F.Gene_Module_tabl, true, true);
   }
   std::vector<double> srt =
     evalPosetConstraints(mutatedModules, F.Poset, F.allPosetG);
   std::vector<double> se =
     evalEpistasis(mutatedModules, F.Epistasis);
 
   // For order effects we need a new vector of mutatedModules:
   if(F.gMOneToOne) {
     mutatedModules = ge.orderEff;
   } else {
     mutatedModules = GeneToModule(ge.orderEff, F.Gene_Module_tabl, false, true);
   }
   
   std::vector<double> so =
     evalOrderEffects(mutatedModules, F.orderE);
 
   // I keep s, srt, se, so separate for now for debugging.
   s.insert(s.end(), srt.begin(), srt.end());
   s.insert(s.end(), se.begin(), se.end());
   s.insert(s.end(), so.begin(), so.end());
   
   return s;
 }
 
 
 
 
 
 
 vector<int> getGenotypeDrivers(const Genotype& ge, const vector<int>& drv) {
   // Returns the actual mutated drivers in a genotype.
   // drv comes from R, and it is the vector with the
   // numbers of the genes, not modules.
   vector<int> presentDrv;
   vector<int> og = allGenesinGenotype(ge);
   set_intersection(og.begin(), og.end(),
 		   drv.begin(), drv.end(),
 		   back_inserter(presentDrv));
   return presentDrv;
 }
 
 

 double evalMutator(const Genotype& fullge,
 		  const std::vector<int>& full2mutator,
 		  const fitnessEffectsAll& muEF,
 		  bool verbose = false) {
   // In contrast to nr_fitness, that sets birth and death, this simply
   // returns the multiplication factor for the mutation rate. This is used
   // by mutationFromParent and mutationFromScratch
 
   // Remember that the fitnessEffectsAll struct for mutator does not use
   // the same mapping from gene names to gene numerical IDs as for
   // fitness. fitnessEffectsAll and its associated algorithms expects the
   // present genes to be indexed as successive integers. That is not
   // necessarily the case if only some of the genes are in the mutator
   // fitnessEffectsAll. So we need to remap the gene numerical IDs.  We
   // could try remapping by gene inside the struct, but painful and
   // error-prone. Much simpler at least for now to do:
 
   // full genotype -> vector of ints (preserving order)
   //     -> convert the ints to the ints for the genotype of mutator
   //     -> genotype in terms of mutator
 
   // the "genotype in terms of mutator" is never preserved. It is just a
   // transient mapping.
 
   // This will NOT work if we ever have order effects for mutator as we do
   // not record order for those that matter for mutator if they do not matter for
   // fitness.
 
   vector<int> g1 = genotypeSingleVector(fullge);
   vector<int> g2;
   int tmp;
   for (auto const & i : g1) {
     tmp = full2mutator[i - 1]; //gives a -9 if no correspondence
     if( tmp > 0 ) g2.push_back(tmp);
   }
   if(g2.size() == 0) {
     return 1.0;
   } else {
     Genotype newg = convertGenotypeFromInts(g2, muEF);
     vector<double> s = evalGenotypeFitness(newg, muEF);
     
     // just for checking
     if(verbose) {
       std::string sprod = "mutator product";
       Rcpp::Rcout << "\n Individual " << sprod << " terms are :";
       for(auto const &i : s) Rcpp::Rcout << " " << i;
       Rcpp::Rcout << std::endl;
     }
     return prodMuts(s);
   }
 }
 
 

 // [[Rcpp::export]]
 double evalRGenotype(Rcpp::IntegerVector rG, Rcpp::List rFE,

 		     bool verbose, bool prodNeg,
 		     Rcpp::CharacterVector calledBy_) {
   // Can evaluate both ONLY fitness or ONLY mutator. Not both at the same
   // time. Use evalRGenotypeAndMut for that.
   const std::string calledBy = Rcpp::as<std::string>(calledBy_);
   

   if(rG.size() == 0) {

     // Why don't we evaluate it?
     Rcpp::warning("WARNING: you have evaluated fitness/mutator status of a genotype of length zero.");

     return 1;
   }

   //const Rcpp::List rF(rFE);
   fitnessEffectsAll F = convertFitnessEffects(rFE);

   Genotype g = convertGenotypeFromR(rG, F);
   vector<double> s = evalGenotypeFitness(g, F);
   if(verbose) {

     std::string sprod;
     if(calledBy == "evalGenotype") {
       sprod = "s";
     } else { // if (calledBy == "evalGenotypeMut") {
       sprod = "mutator product";
     }
     Rcpp::Rcout << "\n Individual " << sprod << " terms are :";

     for(auto const &i : s) Rcpp::Rcout << " " << i;
     Rcpp::Rcout << std::endl;
   }

   if(calledBy == "evalGenotype") {
     if(!prodNeg)
       return prodFitness(s);
     else 
       return prodDeathFitness(s);
   } else { //if (calledBy == "evalGenotypeMut") {
     return prodMuts(s);
   }
 }
 
 
 // [[Rcpp::export]]
 Rcpp::NumericVector evalRGenotypeAndMut(Rcpp::IntegerVector rG,
 					Rcpp::List rFE,
 					Rcpp::List muEF,
 					Rcpp::IntegerVector full2mutator_,
 					bool verbose, bool prodNeg) {
   // Basically to test evalMutator. We repeat the conversion to genotype,
   // but that is unavoidable here.
 
 
   NumericVector out(2);
 
   // For fitness. Except for "evalGenotypeFromR", all is done as in the
   // rest of the internal code for evaluating a genotype.
   fitnessEffectsAll F = convertFitnessEffects(rFE);
   fitnessEffectsAll muef = convertFitnessEffects(muEF);
   Genotype g = convertGenotypeFromR(rG, F);
   vector<double> s = evalGenotypeFitness(g, F);

   if(!prodNeg)

     out[0] = prodFitness(s);

   else 

     out[0] = prodDeathFitness(s);
   if(verbose) {
     std::string sprod = "s";
     Rcpp::Rcout << "\n Individual " << sprod << " terms are :";
     for(auto const &i : s) Rcpp::Rcout << " " << i;
     Rcpp::Rcout << std::endl;
   }
   // out[0] = evalRGenotype(rG, rFE, verbose, prodNeg, "evalGenotype");
   // Genotype fullge = convertGenotypeFromR(rG, F);
   
   const std::vector<int> full2mutator = Rcpp::as<std::vector<int> >(full2mutator_);
   out[1] = evalMutator(g, full2mutator, muef, verbose);
   
   return out;

 }
 
 

 
 double mutationFromScratch(const std::vector<double>& mu,
 			   const spParamsP& spP,
 			   const Genotype& g,
 			   const fitnessEffectsAll& fe,
 			   const int mutationPropGrowth,
 			   const std::vector<int> full2mutator,
 			   const fitnessEffectsAll& muEF) {
   double mumult;
   if(full2mutator.size() > 0) { // so there are mutator effects
     mumult = evalMutator(g, full2mutator, muEF);
   } else mumult = 1.0;
 
   if(mu.size() == 1) {
     if(mutationPropGrowth)
       return(mumult * mu[0] * spP.numMutablePos * spP.birth);
     else
       return(mumult * mu[0] * spP.numMutablePos);
   } else {
     std::vector<int> sortedG = allGenesinGenotype(g);
     std::vector<int> nonmutated;
     set_difference(fe.allGenes.begin(), fe.allGenes.end(),
 		   sortedG.begin(), sortedG.end(),
 		   back_inserter(nonmutated));
     // std::vector<int> mutatedG = genotypeSingleVector(g);
     // Not worth it using an accumulator?
     // std::vector<int> gg = genotypeSingleVector(g);
     // accumulate(gg.begin(), gg.end(), 0.0,
     // 	       [](double x, int y) {return( x + mu[y - 1])});
     double mutrate = 0.0;
     for(auto const &nm : nonmutated) {
       mutrate += mu[nm - 1];
     }
     if(mutationPropGrowth)
       mutrate *= spP.birth;
     return(mumult * mutrate);
   }
 }
 
 

 std::vector < std::vector<int> > list_to_vector_of_int_vectors(Rcpp::List vlist) {
   // As it says. We check each vector is sorted!
   std::vector < std::vector<int> > vv(vlist.size());
   for(int i = 0; i != vlist.size(); ++i) {
     vv[i] = Rcpp::as<std::vector<int> >(vlist[i]);
     if( ! is_sorted(vv[i].begin(), vv[i].end()) )
       throw std::logic_error("Fixation genotypes not sorted. Bug in R code.");
   }
   return vv;
 }
 
 // // [[Rcpp::export]]
 // void wrap_list_to_vector_of_int_vectors(Rcpp::List vlist) {
 //   std::vector < std::vector<int> > vo(vlist.size());
 //   vo = list_to_vector_of_int_vectors(vlist);
 //   for(int ii = 0; ii != vo.size(); ++ii) {
 //     Rcpp::Rcout << "\n";
 //     Rcpp::Rcout << " list position " << ii + 1 << ": ";
 //     for(int jj = 0; jj != vo[ii].size(); ++jj ) {
 //       Rcpp::Rcout << vo[ii][jj] << " ";
 //     }
 //   }
 //   Rcpp::Rcout << "\n";
 // }
 
 
 

 // Wrong when/if there are mutator effects. For suppose there wre, and
 // they affected the parent, but no new mutator gene affects the child.
 // We will, however, multiply twice by the mutator effect.  Therefore, we
 // disable this for now.  We could fix this, checking if there are new
 // mutation effects, or mutliplying/subtracting only new, etc. But too
 // much of a mess.
 // double mutationFromParent(const std::vector<double>& mu,
 // 			  const spParamsP& newP,
 // 			  const spParamsP& parentP,
 // 			  const std::vector<int>& newMutations,
 // 			  // const std::vector<int>& nonmutated,
 // 			  const int mutationPropGrowth,
 // 			  const Genotype& fullge,
 // 			  const std::vector<int> full2mutator,
 // 			  const fitnessEffectsAll& muEF) {
 //   double mumult;
 //   if(full2mutator.size() > 0) { // so there are mutator effects
 //     mumult = evalMutator(fullge, full2mutator, muEF);
 //   } else mumult = 1.0;
   
 //   if(mu.size() == 1) {
 //     if(mutationPropGrowth)
 //       return(mumult * mu[0] * newP.numMutablePos * newP.birth);
 //     else
 //       return(mumult * mu[0] * newP.numMutablePos);
 //   } else {
 //     double mutrate = parentP.mutation;
 //     for(auto const mutated : newMutations) {
 //       mutrate -= mu[mutated - 1];
 //     }
 //     if(mutationPropGrowth)
 //       mutrate *= newP.birth;
 //     return(mumult * mutrate);
 //   }
 // }
 
 
   
 // About order of genes and their names, etc
 
 // We first read the R gene module. The $geneModule. The function is
 // R_GeneModuleToGeneModule
 
 // We also read the no interaction. They have their own number-name
 // correspondence, within the noInt genes part. See the struct
 // genesWithoutInt.  But that already comes order from R with numbers
 // starting after the last gene with interaction. See the R function
 // allFitnessEffects.