// 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 fitness landscape, a single s is passed, and that has a - 1 subtracted
// 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)));});
}
// // 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>());
}
// new cPDetect
double set_cPDetect(const double n2, const double p2,
const double PDBaseline) {
return ( -log(1.0 - p2) * (PDBaseline / (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)/PDBaseline ) ));
}
}
// Former cpDetect mechanism
// 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)));
// // }
// // }
// is this detected, 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 == "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_EFVMap(const std::map<std::string, double>& efv) {
Rcpp::Rcout << "\n Printing evalFVars_struct\n";
for(auto elem : efv) {
Rcpp::Rcout << elem.first << "\t " << elem.second << "\n";
}
}
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"]);
flS.flFDFmap.clear();//Set to 0 flFDFmap
flS.flfVarsmap.clear();//Set to 0 flfVarsmap
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;
}
// FIXME
// This is very limited: a symbol can only be a symbol if it is
// in the fitness landscape table. But we wanted to allow sparse fitness landscape tables
// where we just don't add genotypes with fitness 0.
// Or is it the logic that a fitness of 0 genotype will never appear in an equation?
// flg = rFE["fitnessLandscape_gene_id"];
// fl_df = rFE["fitnessLandscape_df"];
// fvars = rFE["fitnessLandscapeVariables"];
// FIXME: should use the full_FDF_spec dataframe in R. That shows the full mapping.
// And we could use the genotypes as strings? maybe not this
fitnessLandscape_struct convertFitnessLandscape_fdf(Rcpp::List flg,
Rcpp::List fl_df,
Rcpp::StringVector fvars) {
fitnessLandscape_struct flS;
flS.names = Rcpp::as<std::vector<std::string> >(flg["Gene"]);
flS.NumID = Rcpp::as<std::vector<int> >(flg["GeneNumID"]);
flS.flmap.clear();//Set to 0 flmap
std::vector<std::string> fvarsvect =
Rcpp::as<std::vector<std::string> > (fvars);
// From $fitnessLandscape_df
std::vector<std::string> genotNames =
Rcpp::as<std::vector<std::string> >(fl_df["Genotype"]);
// From $fitnessLandscape_df
std::vector<std::string> fitness =
Rcpp::as<std::vector<std::string> > (fl_df["Fitness"]);
if(fvarsvect.size() != genotNames.size() )
throw std::logic_error("fvarsvect (fitnessLandscapeVariables) and "
"genotNames (fitnessLandscape_df$Genotypes) "
"are of different lenght. Should have been caught in R");
// Fill up the map genotypes (as string) to fitness (as string)
// Length given by $fitnessLandscape_df
for(size_t i = 0; i != genotNames.size(); ++i) {
flS.flFDFmap.insert({genotNames[i], fitness[i]});
}
// Fill up the map genotypes (as string) 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.flfVarsmap.insert({genotNames[i], fvarsvect[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::StringVector fvars = rFE["fitnessLandscapeVariables"];
bool fdf = as<bool>(rFE["frequencyDependentFitness"]);
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()) {
//New block to deal with fdf stuff
if(fdf){
fe.fitnessLandscape = convertFitnessLandscape_fdf(flg, fl_df, fvars);
} else {
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
fe.fVars = as<std::vector<std::string> > (fvars); //new line to insert fVars
//fe.fVars = fvariables;
fe.frequencyDependentFitness = fdf;//new line to insert frequencyDependentFitness
//fe.frequencyType = as<std::string> (fType);
fe.frequencyType = fType;
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!!!
// 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;
}
// 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;
}
double frequency(const int& pos, const std::vector<spParamsP>& popParams){
double fqc;
double numerator = popParams[pos].popSize;
double denominator = 0.0;
for(size_t i = 0; i < popParams.size(); i++){
denominator += popParams[i].popSize;
}
fqc = numerator / denominator;
return fqc;
}
//This function returns 0 if genotype is not in Genotypes and position+1 otherwise
int findPositionInGenotypes(const std::vector<Genotype>& Genotypes,
const std::vector<int> genotype){
int index;
std::vector<std::vector<int> > flGenesInGenotypes;
for(size_t i = 0; i < Genotypes.size(); i++){
flGenesInGenotypes.push_back(Genotypes[i].flGenes);
}
int pos = std::find(flGenesInGenotypes.begin(),
flGenesInGenotypes.end(), genotype) - flGenesInGenotypes.begin();
int size = flGenesInGenotypes.size();
if(pos < size)
index = pos + 1;
else
index = 0;
return index;
}
//This function works oposite to concatIntsString #include <sstream>
std::vector<int> stringVectorToIntVector(const std::string str){
std::vector<int> vector;
std::stringstream ss(str);
int i;
while (ss >> i) {
vector.push_back(i);
if(ss.peek() == ','|| ss.peek() == ' '){
ss.ignore();
}
}
return vector;
}
//This function produce the map (structure) that links fVars (keys) to its frequencies (values)
evalFVars_struct evalFVars(const fitnessEffectsAll& F,
const std::vector<Genotype>& Genotypes,
const std::vector<spParamsP>& popParams){
evalFVars_struct efvs;
std::map<std::string, std::string> fvarsmap = F.fitnessLandscape.flfVarsmap;
std::string freqType = F.frequencyType;
for(const auto& iterator : fvarsmap) {
std::vector<int> genotype = stringVectorToIntVector(iterator.first);//genotype (as int vector)
std::string var = iterator.second;//variable associated to genotype
int position = findPositionInGenotypes(Genotypes, genotype);
if(position != 0){
int realPos = position - 1;
if(freqType == "abs"){
double freqAbs = popParams[realPos].popSize;
efvs.evalFVarsmap.insert({var, freqAbs});
} else {
double freqRel = frequency(realPos, popParams);
efvs.evalFVarsmap.insert({var, freqRel});
}
} else {
double freq = 0.0;
efvs.evalFVarsmap.insert({var, freq});
}
}
return efvs;
}
double totalPop(const std::vector<spParamsP>& popParams){
double sum = 0.0;
for(size_t i = 0; i < popParams.size(); i++){
sum += popParams[i].popSize;
}
return sum;
}
double evalGenotypeFDFitnessEcuation(const Genotype& ge,
const fitnessEffectsAll& F,
const std::vector<Genotype>& Genotypes,
const std::vector<spParamsP>& popParams,
const double& currentTime){
double f;
evalFVars_struct symbol_table_struct = evalFVars(F, Genotypes, popParams);
std::map<std::string, double> EFVMap = symbol_table_struct.evalFVarsmap;
// print_EFVMap(EFVMap); debugging
std::string gs = concatIntsString(ge.flGenes);
std::string expr_string = F.fitnessLandscape.flFDFmap.at(gs);
double N = totalPop(popParams);
double T = currentTime;
typedef exprtk::symbol_table<double> symbol_table_t;
typedef exprtk::expression<double> expression_t;
typedef exprtk::parser<double> parser_t;
symbol_table_t symbol_table;
for(auto& iterator : EFVMap) {
symbol_table.add_variable(iterator.first, iterator.second);
}
symbol_table.add_constant("N", N);//We reserve N to total population size
symbol_table.add_constant("T", T); //Pass current time to exprtk
symbol_table.add_constants();
expression_t expression;
expression.register_symbol_table(symbol_table);
parser_t parser;
if (!parser.compile(expr_string, expression)){
Rcpp::Rcout << "\nexprtk parser error: \n" << std::endl;
for (std::size_t i = 0; i < parser.error_count(); ++i){
typedef exprtk::parser_error::type error_t;
error_t error = parser.get_error(i);
// RDU: FIXME?
// Rcpp::Rcout <<
REprintf("Error[%02zu] Position: %02zu Type: [%14s] Msg: %s Expression: %s\n",
i,
error.token.position,
exprtk::parser_error::to_str(error.mode).c_str(),
error.diagnostic.c_str(),
expr_string.c_str());
// << std::endl;
}
std::string errorMessage1 = "Wrong evalGenotypeFDFitnessEcuation evaluation, ";
std::string errorMessage2 = "probably bad fitness columm especification.";
std::string errorMessage = errorMessage1 + errorMessage2;
throw std::invalid_argument(errorMessage);
}
f = expression.value();
return f;
}
std::vector<double> evalGenotypeFitness(const Genotype& ge,
const fitnessEffectsAll& F,
const std::vector<Genotype>& Genotypes,
const std::vector<spParamsP>& popParams,
const double& currentTime){
// check_disable_later
checkLegitGenotype(ge, F);
std::vector<double> s;
if( ((ge.orderEff.size() + ge.epistRtEff.size() +
ge.rest.size() + ge.flGenes.size() ) == 0) && !F.frequencyDependentFitness ) {
// Rcpp::Rcout << "NOTE: you have evaluated fitness of a genotype of length zero (WT?) when non-fdf-fitness. It is 1 by decree. \n";
// Empty and have prodFitness deal with it.
// Why not return a 0 (1 - 1)? FIXME: zz5
// s.push_back(0.0); return s;
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.frequencyDependentFitness){ //possible also with Genotype.size()==0 and popParams.size==0 ?
if(F.fitnessLandscape.flFDFmap.find(gs) == F.fitnessLandscape.flFDFmap.end()) {
s.push_back(-1.0);
} else {
s.push_back(evalGenotypeFDFitnessEcuation(ge, F, Genotypes, popParams, currentTime) - 1);
}
} else {
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);
}
}
}
// 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;
}
// FIXME: refactor?
// Why does this function take as arguments popParams and
// Genotypes? It shouldn't, since its behavior does not depend on those.
// It does because we call evalGenotypeFitness which calls in turn
// evalGenotypeFDFitnessEcuation, etc. But we probably want to do this
// differently. From creating a new function, using overloading, to do s =
// evalGenotypeFitness, that would be double newfun(A, B) =
// evalGenotypeFitness(A, B); and maybe having two different
// evalGenotypeFitness.
// This same problem happens with currentTime; it is used by evalGenotypeFDFitnessEcuation,
// etc
double evalMutator(const Genotype& fullge,
const std::vector<int>& full2mutator,
const fitnessEffectsAll& muEF,
const std::vector<Genotype>& Genotypes,
const std::vector<spParamsP>& popParams,
const double& currentTime, // because we use it in evalGenotypeFitness
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, Genotypes, popParams, currentTime);
// 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);
}
}
std::vector<Genotype> genotypesFromScratch(const std::vector<std::string>& genotNames){
std::vector<Genotype> genotypesVector;
for(size_t i = 0; i != genotNames.size(); ++i){
//Genotype genotype = wtGenotype();
Genotype genotype;
genotype.flGenes = stringVectorToIntVector(genotNames[i]);
genotypesVector.push_back(genotype);
}
return genotypesVector;
}
std::vector<spParamsP> popParamsFromScratch(const std::vector<int>& spPopSizes){
std::vector<spParamsP> popParamsVector;
for(size_t i = 0; i != spPopSizes.size(); ++i){
spParamsP spparams;
spparams.popSize = spPopSizes[i];
popParamsVector.push_back(spparams);
}
return popParamsVector;
}
// [[Rcpp::export]]
double evalRGenotype(Rcpp::IntegerVector rG,
Rcpp::List rFE,
Rcpp::IntegerVector spPop,
bool verbose,
bool prodNeg,
Rcpp::CharacterVector calledBy_,
double currentTime) {
// 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_);
const bool fdf = as<bool>(rFE["frequencyDependentFitness"]);
if(rG.size() == 0 && fdf == false) {
// Why don't we evaluate it?
Rcpp::Rcout << "NOTE: you have evaluated fitness/mutator status of a genotype of length zero (WT?) in non fdf fitness. It is 1 by decree. \n";
return 1;
}
std::vector<Genotype> Genotypes;
std::vector<spParamsP> popParams;
if(fdf){
//std::vector<int> spPopSizes;
//spPopSizes = as<std::vector<int> > (rFE["spPopSizes"]);
std::vector<int> spPopSizes = as<std::vector<int> > (spPop);
Rcpp::List fl_df = rFE["fitnessLandscape_df"];
std::vector<std::string> genotNames = Rcpp::as<std::vector<std::string> >(fl_df["Genotype"]);
Genotypes = genotypesFromScratch(genotNames);
popParams = popParamsFromScratch(spPopSizes);
}//else{
//const std::vector<Genotype> Genotypes(0);
//const std::vector<spParamsP> popParams(0);
//}
//const Rcpp::List rF(rFE);
fitnessEffectsAll F = convertFitnessEffects(rFE);
Genotype g = convertGenotypeFromR(rG, F);
vector<double> s = evalGenotypeFitness(g, F, Genotypes, popParams, currentTime);
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 spPop,
Rcpp::IntegerVector full2mutator_,
bool verbose,
bool prodNeg,
double currentTime) {
// Basically to test evalMutator. We repeat the conversion to genotype,
// but that is unavoidable here.
const bool fdf = as<bool>(rFE["frequencyDependentFitness"]);
/*
if(rG.size() == 0 && fdf == false) {
// Why don't we evaluate it?
Rcpp::warning("WARNING: you have evaluated fitness/mutator status of a genotype of length zero.");
return 1;
}
*/
NumericVector out(2);
std::vector<Genotype> Genotypes;
std::vector<spParamsP> popParams;
if(fdf){
//std::vector<int> spPopSizes;
//spPopSizes = as<std::vector<int> > (rFE["spPopSizes"]);
std::vector<int> spPopSizes = as<std::vector<int> > (spPop);
Rcpp::List fl_df = rFE["fitnessLandscape_df"];
std::vector<std::string> genotNames = Rcpp::as<std::vector<std::string> >(fl_df["Genotype"]);
Genotypes = genotypesFromScratch(genotNames);
popParams = popParamsFromScratch(spPopSizes);
}
// 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, Genotypes, popParams, currentTime);
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, Genotypes, popParams, currentTime, verbose);
return out;
}
// FIXME refactor
// as in evalMutator, that this takes Genotypes and popParams is arguably
// bad design and very confusing (those arguments have nothing to do with mutations)
// Same problem happens with currentTime. Yes, this sucks.
// return mutation rate for a new genotype; set to dummyMutationRate in specific cases
// do not use max w.r.t dummyMutationRate
// indiscriminately: we want to see
// dummyMutationRate being set
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,
const std::vector<Genotype>& Genotypes,
const std::vector<spParamsP>& popParams,
const double& currentTime,
const double& dummyMutationRate) {
double tmp;
double mumult;
if(full2mutator.size() > 0) { // so there are mutator effects
mumult = evalMutator(g, full2mutator, muEF, Genotypes, popParams, currentTime);
} else mumult = 1.0;
//FIXME: here the code for altering mutation rate
// with a procedure like ExprTk for fitness??
// Nope: alter directly spParams.mutation.??
// where the updateRatesFDF... are called.
// In BNB_nr.cpp, in nr_innerBNB function
// when we are sampling.
if(mu.size() == 1) {
if(spP.numMutablePos == 0) {
// Standard modus operandi. No need to flood the user with messages.
// Rcpp::Rcout << "mFS-message-1: "
// << "No mutable positions. Mutation set to dummyMutationRate "
// << dummyMutationRate << "\n";
return(dummyMutationRate);
}
if(mutationPropGrowth) {
tmp = mumult * mu[0] * spP.numMutablePos * spP.birth;
if(tmp <= dummyMutationRate) {
Rcpp::Rcout << "mFS-messagesMPL: Mutable positions left: ";
if(mumult == 1.0) {
// letters match codes for varmutrate
Rcpp::Rcout << "mFS-message-2-B: constant mut rate "
<< "no mutator and mutationPropGrowth "
<< ". birth rate = " << spP.birth << "\n";
} else {
Rcpp::Rcout << "mFS-message-2-C: constant mut rate "
<< " mutator and mutationPropGrowth "
<< ". birth rate = " << spP.birth
<< ". mumult = " << mumult << "\n";
}
Rcpp::Rcout << "\n mutation rate = " << tmp << " < dummyMutationRate "
<< dummyMutationRate
<< ". Expect numerical problems.\n";
// Should set mutation rate to dummyMutationRate? Remember genotypes of birth rate 0 never created
}
return(tmp);
} else {
tmp = mumult * mu[0] * spP.numMutablePos;
if(tmp <= dummyMutationRate) {
Rcpp::Rcout << "mFS-messagesMPL: Mutable positions left: ";
if(mumult == 1.0) {
// letters match codes for varmutrate
Rcpp::Rcout << "mFS-message-2-A: constant mut rate "
<< "no mutator and no mutationPropGrowth ";
} else {
Rcpp::Rcout << "mFS-message-2-D: constant mut rate"
<< " mutator and no mutationPropGrowth "
<< ". mumult = " << mumult << "\n";
}
Rcpp::Rcout << "\n mutation rate = " << tmp << " < dummyMutationRate "
<< dummyMutationRate
<< ". Expect numerical problems.\n";
}
return(tmp);
}
} 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(mutrate == 0) {
Rcpp::Rcout << "mFS-message-4 . No mutable positions? Mutation set to dummyMutationRate "
<< dummyMutationRate << "\n";
return(dummyMutationRate);
}
if(mutationPropGrowth)
mutrate *= spP.birth;
tmp = mumult * mutrate;
if(tmp <= dummyMutationRate) {
Rcpp::Rcout << "mFS-messagesMPL: Mutable positions left: ";
if( (mumult == 1.0) && (!mutationPropGrowth) ) {
Rcpp::Rcout << "mFS-message-5-A: variable mut rate "
<< "no mutator and no mutationPropGrowth\n ";
} else if ((mumult == 1.0) && (mutationPropGrowth) ) {
Rcpp::Rcout << "mFS-message-5-B: variable mut rate "
<< "no mutator and mutationPropGrowth "
<< ". birth rate = " << spP.birth << "\n";
} else if ( (mumult != 1.0) && (mutationPropGrowth) ) {
Rcpp::Rcout << "mFS-message-5-C: variable mut rate "
<< "mutator and mutationPropGrowth "
<< ". birth rate = " << spP.birth
<< ". mumult = " << mumult << "\n";
} else if ( (mumult != 1.0) && (!mutationPropGrowth) ) {
Rcpp::Rcout << "mFS-message-5-D: variable mut rate "
<< "mutator and no mutationPropGrowth "
<< ". mumult = " << mumult << "\n";
} else {
throw std::logic_error("\n This case should not exist: mFS-except\n");
}
Rcpp::Rcout << "\n mutation rate = " << tmp << " < dummyMutationRate "
<< dummyMutationRate << ". Expect numerical problems.\n";
}
return(tmp);
}
}
// list -> vector of int vectors.
// If input is of length zero, return object of length 0
// No checks are made of whether this makes sense (undefined behav
// if this is a list of floats, strings, etc)
// Check for order is necessary for fixation genotypes but not for init mutant.
std::vector < std::vector<int> > list_to_vector_of_int_vectors(Rcpp::List vlist,
bool check_ordered = true) {
// As it says. We check each vector is sorted!
std::vector < std::vector<int> > vv(vlist.size());
if (vlist.size() != 0) {
for(int i = 0; i != vlist.size(); ++i) {
vv[i] = Rcpp::as<std::vector<int> >(vlist[i]);
if( check_ordered && (! is_sorted(vv[i].begin(), vv[i].end())) )
throw std::logic_error("Fixation genotypes not sorted. Bug in R code.");
}
} else {
vv.resize(0);
}
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.
// FIXME in new code with time, why is there this line?
// double T = popParams[0].timeLastUpdate;
// if (T == std::numeric_limits<double>::infinity()
// or T == -std::numeric_limits<double>::infinity()) {
// T = 0;}
