deleted file mode 100644

@@ -1,1933 +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/>.

-

-

-// #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 "exprtk.h"

-#include <Rcpp.h>

-#include <iomanip>

-#include <algorithm>

-#include <random>

-#include <string>

-#include <sstream>

-#include <limits>

-#include <regex>

-

-

-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;

-}

-

-// Cumulative product of (1 + s).

-//      If using fitness landscape, a single s is passed, and that already

-//      has a "-1" subtracted. So when a single s is passed, we just

-//      return s + 1, which is the fitness (birth rate) as specified in

-//      the fitness landscape. This applies to FDF too.

-//      See function evalGenotypeFitness.

-//      Whenever we call prodFitness (e.g., nr_fitness or

-//      initMutantInitialization) it is called on the output of

-//      evalGenotypeFitness

-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)));});

-}

-

-

-// Plays same role as prodFitness, but for Bozic models

-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)));});

-}

-

-

-// To compute cumulative product of mutator effects

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

-  return accumulate(s.begin(), s.end(), 1.0,

-		    std::multiplies<double>());

-}

-

-// Compute cPDetect for Detection/stopping mechanism

-double set_cPDetect(const double n2, const double p2,

-		    const double PDBaseline) {

-  return ( -log(1.0 - p2) *  (PDBaseline / (n2 - PDBaseline)) );

-}

-

-// Probability of detection given size

-double probDetectSize(const double n, const double cPDetect,

-		      const double PDBaseline) {

-  if(n <= PDBaseline) {

-    return 0;

-  } else {

-    return (1 - exp( -cPDetect * ( (n - PDBaseline)/PDBaseline ) ));

-  }

-}

-

-

-// return true if this is detected (as given by the probability of detection as a function of size)

-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 == "mcfarlandlogd")

-    return TypeModel::mcfarlandlog_d;

-  else if (mod == "arbitrary")

-    return TypeModel::arbitrary;

-  else if (mod == "constant")

-    return TypeModel::constant;

-  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.

-}

-

-

-

-// genotype struct -> genotype as a vector of ints

-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(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 &b: F.genesNoInt.NumID) {

-    g0.push_back(b);

-  }

-  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.

-

-  // 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"]);

-

-    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])) );

-    }

-

-    // 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()); // 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.");

-    // remove the next two 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::List fl_df, Rcpp::List full_FDF_spec,

-                bool fdb, bool fdd){

-

-

-  // We assume this:

-  // A symbol can only be a symbol if it is

-  // in the fitness landscape table. A genotype not in the fitness landscape table

-  // always has fitness zero, and noting can depend on it.

-  // A genotype of 0 genotype will never appear in an equation.

-

-

-  fitnessLandscape_struct flS;

-

-  flS.names = Rcpp::as<std::vector<std::string> >(flg["Gene"]);

-  flS.NumID = Rcpp::as<std::vector<int> >(flg["GeneNumID"]);

-

-  if (fdb) {

-

-    flS.flbmap.clear(); //Set to 0 flmap

-    std::vector<std::string> genotNames =

-      Rcpp::as<std::vector<std::string> >(full_FDF_spec["Genotype_as_numbers"]);

-

-    std::vector<std::string> fvarsbvect =

-    Rcpp::as<std::vector<std::string> > (full_FDF_spec["Genotype_as_fvarsb"]);

-

-    std::vector<std::string> birth =

-      Rcpp::as<std::vector<std::string> > (full_FDF_spec["Birth_as_fvars"]);

-

-    if(fvarsbvect.size() != genotNames.size() )

-      throw std::logic_error("fvarsbvect (fitnessLandscapeVariables) and "

-          "genotNames (fitnessLandscape_df$Genotypes) "

-          "are of different lenght. Should have been caught in R");

-    

-    // Fill up the map genotypes (as string of ints) to fitness (as string:

-    // the expression that exprTk will evaluate).

-    // Length given by number of genotypes in fitness landscape

-    for(size_t i = 0; i != genotNames.size(); ++i) {

-      flS.flFDBmap.insert({genotNames[i], birth[i]});

-    }

-

-    // Fill up the map genotypes (as string of ints) to fVars (as string)

-    // If a mismatch in $fitnessLandscape_df and

-    // $fitnessLandscapeVariables, this does silly things.

-    for(size_t i = 0; i != genotNames.size(); ++i) {

-      flS.flfVarsBmap.insert({genotNames[i], fvarsbvect[i]});

-    }

-  } else {

-

-    flS.flFDBmap.clear();//Set to 0 flFDBmap

-    flS.flfVarsBmap.clear();//Set to 0 flfVarsBmap

-

-    std::vector<std::string> genotNames =

-      Rcpp::as<std::vector<std::string> >(fl_df["Genotype"]);

-    Rcpp::NumericVector birth = fl_df["Birth"];

-

-    for(size_t i = 0; i != genotNames.size(); ++i) {

-      flS.flbmap.insert({genotNames[i], birth[i]});

-    }

-  }

-

-  // Check if death is specified

-  if (fl_df.containsElementNamed("Death")) {

-

-    if(fdd) {

-

-      flS.fldmap.clear(); //Set to 0 fldmap

-

-      std::vector<std::string> genotNames =

-      Rcpp::as<std::vector<std::string> >(full_FDF_spec["Genotype_as_numbers"]);

-

-      std::vector<std::string> fvarsdvect =

-      Rcpp::as<std::vector<std::string> > (full_FDF_spec["Genotype_as_fvarsd"]);

-

-      std::vector<std::string> death =

-        Rcpp::as<std::vector<std::string> > (full_FDF_spec["Death_as_fvars"]);

-

-

-      if(fvarsdvect.size() != genotNames.size() )

-        throw std::logic_error("fvarsdvect (fitnessLandscapeVariables) and "

-            "genotNames (fitnessLandscape_df$Genotypes) "

-            "are of different lenght. Should have been caught in R");

-      

-      // Fill up the map genotypes (as string of ints) to fitness (as string:

-      // the expression that exprTk will evaluate).

-      // Length given by number of genotypes in fitness landscape

-      for(size_t i = 0; i != genotNames.size(); ++i) {

-        flS.flFDDmap.insert({genotNames[i], death[i]});

-      }

-

-      // Fill up the map genotypes (as string of ints) to fVars (as string)

-      // If a mismatch in $fitnessLandscape_df and

-      // $fitnessLandscapeVariables, this does silly things.

-      for(size_t i = 0; i != genotNames.size(); ++i) {

-        flS.flfVarsDmap.insert({genotNames[i], fvarsdvect[i]});

-      }

-    } else {

-

-      flS.flFDDmap.clear();//Set to 0 flFDDmap

-      flS.flfVarsDmap.clear();//Set to 0 flfVarsDmap

-

-      std::vector<std::string> genotNames =

-        Rcpp::as<std::vector<std::string> >(fl_df["Genotype"]);

-

-      Rcpp::NumericVector death = fl_df["Death"];

-

-      for(size_t i = 0; i != genotNames.size(); ++i) {

-        flS.fldmap.insert({genotNames[i], death[i]});

-      }

-    }

-

-  } else {

-    flS.flFDDmap.clear();//Set to 0 flFDDmap

-    flS.flfVarsDmap.clear();//Set to 0 flfVarsDmap

-    flS.fldmap.clear(); //Set to 0 fldmap

-  }

-

-  return flS;

-}

-

-genesWithoutInt convertNoInts(Rcpp::List nI) {

-  genesWithoutInt genesNoInt;

-  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]; // 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"];

-  bool fdb = as<bool>(rFE["frequencyDependentBirth"]);

-  bool fdd = as<bool>(rFE["frequencyDependentDeath"]);

-  std::string fType = as<std::string>(rFE["frequencyType"]);

-  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.

-  //<std::vector<std::string> > fvariables = as<std::vector<std::string> > (fvars);

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

-  if(fl_df.size()) {

-    

-    Rcpp::List full_FDF_spec;

-    if (fdb || fdd) {

-      // Should not be used in convertFitnessLandscape if not fdb and not fdd

-      full_FDF_spec = rFE["full_FDF_spec"];

-    }

-

-    fe.fitnessLandscape = convertFitnessLandscape(flg, fl_df, full_FDF_spec, fdb, fdd); 

-    if (fdb) {

-      fe.fVarsb = as<std::vector<std::string> > (full_FDF_spec["Genotype_as_fvarsb"]);

-      fe.frequencyType = fType;

-    }

-

-    if (fdd) {

-      fe.fVarsd = as<std::vector<std::string> > (full_FDF_spec["Genotype_as_fvarsd"]);

-      fe.frequencyType = fType;

-    }

-  }

-

-  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.size()) == 0) {

-      throw std::logic_error("\n Nothing inside this fitnessEffects; why are you here?"

-			     "  Bug in R code.");

-  }

-

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

-  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");

-  }

-

-  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

-

-  fe.frequencyDependentBirth = fdb; //new line to insert frequencyDependentBirth

-  fe.frequencyDependentDeath = fdd; //new line to insert frequencyDependentDeath

-

-  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

-    //chromothripsis 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.

-    // FIXMEmaybe: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;

-  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));

-  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]});

-  }

-

-  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;

-

-  // 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

-

-  // initMutant cannot use this: we give the order.

-  // That is why the call from initMutant uses random = false

-  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;

-

-  // 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;