// Copyright 2013, 2014, 2015, 2016 Ramon Diaz-Uriarte
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#include <Rcpp.h>
// using namespace Rcpp;
inline int HammingDistance(const Rcpp::IntegerVector& x, const Rcpp::IntegerVector& y) {
Rcpp::NumericVector diff = Rcpp::abs( x - y );
return std::accumulate(diff.begin(), diff.end(), 0);
}
// eventually remove this. Left now for testing
// [[Rcpp::export]]
Rcpp::IntegerVector accessibleGenotypes_former(Rcpp::IntegerMatrix y,
Rcpp::NumericVector f,
Rcpp::IntegerVector numMut, //
double th) {
// Return just the indices. Could preserve fitness, but would need
// another matrix.
int ng = y.nrow(); //it counts the wt
Rcpp::IntegerMatrix adm(ng, ng);
int numMutdiff = 0;
// I would have thought this would be faster. It ain't.
// The last genotype never accesses anything.
// for(int i = 0; i < (ng - 1); ++i) {
// // Candidate genotypes to be accessed from i are always of larger
// // mutation by 1. And candidates can thus not have smaller index
// for(int j = (i + 1); j < ng; ++j) {
// if( (numMut(j) == (numMut(i) + 1)) &&
// ( (f(j) - f(i)) >= th) &&
// (HammingDistance(y(i, _), y(j, _)) == 1) ) {
// adm(i, j) = 1;
// } else if( (numMut(j) > (numMut(i) + 1)) ) {
// break;
// }
// }
// }
// The last genotype never accesses anything.
for(int i = 0; i < (ng - 1); ++i) {
// Candidate genotypes to be accessed from i are always of larger
// mutation by 1. And candidates can thus not have smaller index
for(int j = (i + 1); j < ng; ++j) {
numMutdiff = numMut(j) - numMut(i);
if( numMutdiff > 1) { // no more to search
break;
} else if(numMutdiff == 1) {
// f(j) - f(i) is faster than HammingDistance
// but might lead to more evals?
// or fewer, depending on landscape
if( ( (f(j) - f(i)) >= th) &&
(HammingDistance(y(i, Rcpp::_), y(j, Rcpp::_)) == 1)
) {
adm(i, j) = 1;
// Rcpp::Rcout << "i = " << i << " j = " << j << " adm " << adm(i,j) << "\n";
}
}
}
}
// Slightly different logic from R: Do not resize object; set the row to
// 0.
int colsum = 0;
// int indicator = 0; // indicator != 0 means we set one row to 0
// so we need to iterate at least once more.
// accessible is the genotype number, not the column! WT is 1,
// etc. This makes it easy to keep track of which are accessible.
Rcpp::IntegerVector accessible = Rcpp::seq_len(ng);
// This is doable in one pass
// while (true) {
// indicator = 0;
for(int k = 1; k < ng; ++k) {
if(accessible(k) > 0) {
colsum = std::accumulate(adm(Rcpp::_, k).begin(),
adm(Rcpp::_, k).end(), 0);
if(colsum == 0) { // This genotype ain't reachable
// Nothing can be reached from this genotype; fill with 0.
adm(k, Rcpp::_) = Rcpp::IntegerVector(ng);
accessible(k) = -9;
// indicator = 1;
}
}
}
// if(indicator == 0) break;
// }
return accessible;
}
// [[Rcpp::export]]
Rcpp::NumericMatrix genot2AdjMat(Rcpp::IntegerMatrix y,
Rcpp::NumericVector f,
Rcpp::IntegerVector numMut) {
// Return just the indices. Could preserve fitness, but would need
// another matrix.
int ng = y.nrow(); //it counts the wt
Rcpp::NumericMatrix adm(ng, ng);
// fill with NAs: https://stackoverflow.com/a/23753626
// Filling with NAs and in general having NAs might lead to performance
// penalties. But I use the NAs in a lot of the code for accessible
// genotypes, etc.
std::fill( adm.begin(), adm.end(), Rcpp::NumericVector::get_na() ) ;
int numMutdiff = 0;
// I would have thought this would be faster. It ain't.
// The last genotype never accesses anything.
// for(int i = 0; i < (ng - 1); ++i) {
// // Candidate genotypes to be accessed from i are always of larger
// // mutation by 1. And candidates can thus not have smaller index
// for(int j = (i + 1); j < ng; ++j) {
// if( (numMut(j) == (numMut(i) + 1)) &&
// ( (f(j) - f(i)) >= th) &&
// (HammingDistance(y(i, _), y(j, _)) == 1) ) {
// adm(i, j) = 1;
// } else if( (numMut(j) > (numMut(i) + 1)) ) {
// break;
// }
// }
// }
// The last genotype never accesses anything.
for(int i = 0; i < (ng - 1); ++i) {
// Candidate genotypes to be accessed from i are always of larger
// mutation by 1. And candidates can thus not have smaller index
for(int j = (i + 1); j < ng; ++j) {
numMutdiff = numMut(j) - numMut(i);
if( numMutdiff > 1) { // no more to search
break;
} else if(numMutdiff == 1) {
if( HammingDistance(y(i, Rcpp::_), y(j, Rcpp::_)) == 1) {
adm(i, j) = (f(j) - f(i));
}
}
}
}
return adm;
}
Rcpp::IntegerMatrix integerAdjMat(Rcpp::IntegerMatrix y,
Rcpp::NumericVector f,
Rcpp::IntegerVector numMut, //
double th) {
// Return a genotype adjacency matrix with a 1 if genotype j is
// accessible (fitness >, within th) from i.
int ng = y.nrow(); //it counts the wt
Rcpp::IntegerMatrix adm(ng, ng);
int numMutdiff = 0;
// I would have thought this would be faster. It ain't.
// The last genotype never accesses anything.
// for(int i = 0; i < (ng - 1); ++i) {
// // Candidate genotypes to be accessed from i are always of larger
// // mutation by 1. And candidates can thus not have smaller index
// for(int j = (i + 1); j < ng; ++j) {
// if( (numMut(j) == (numMut(i) + 1)) &&
// ( (f(j) - f(i)) >= th) &&
// (HammingDistance(y(i, _), y(j, _)) == 1) ) {
// adm(i, j) = 1;
// } else if( (numMut(j) > (numMut(i) + 1)) ) {
// break;
// }
// }
// }
// The last genotype never accesses anything.
for(int i = 0; i < (ng - 1); ++i) {
// Candidate genotypes to be accessed from i are always of larger
// mutation by 1. And candidates can thus not have smaller index
for(int j = (i + 1); j < ng; ++j) {
numMutdiff = numMut(j) - numMut(i);
if( numMutdiff > 1) { // no more to search
break;
} else if(numMutdiff == 1) {
// f(j) - f(i) is faster than HammingDistance
// but might lead to more evals?
// or fewer, depending on landscape
if( ( (f(j) - f(i)) >= th) &&
(HammingDistance(y(i, Rcpp::_), y(j, Rcpp::_)) == 1)
) {
adm(i, j) = 1;
// Rcpp::Rcout << "i = " << i << " j = " << j << " adm " << adm(i,j) << "\n";
}
}
}
}
return adm;
}
// used in both peaks and accessible genotypes
Rcpp::IntegerVector accessibleGenotypesPeaksLandscape(Rcpp::IntegerMatrix y,
Rcpp::NumericVector f,
Rcpp::IntegerVector numMut, //
double th,
bool returnpeaks) {
// Return the indices. This is like accessibleGenotypes, but we do an
// extra loop
int ng = y.nrow(); //it counts the wt
Rcpp::IntegerMatrix adm(ng, ng);
adm = integerAdjMat(y, f, numMut, th);
int numMutdiff = 0;
// Slightly different logic from R: Do not resize object; set the row to
// 0.
int colsum = 0;
// int indicator = 0; // indicator != 0 means we set one row to 0
// so we need to iterate at least once more.
// accessible is the genotype number, not the column! WT is 1,
// etc. This makes it easy to keep track of which are accessible.
Rcpp::IntegerVector accessible = Rcpp::seq_len(ng);
// This is doable in one pass
// while (true) {
// indicator = 0;
for(int k = 1; k < ng; ++k) {
if(accessible(k) > 0) {
colsum = std::accumulate(adm(Rcpp::_, k).begin(),
adm(Rcpp::_, k).end(), 0);
if(colsum == 0) { // This genotype ain't reachable
// Nothing can be reached from this genotype; fill with 0.
adm(k, Rcpp::_) = Rcpp::IntegerVector(ng);
accessible(k) = -9;
// indicator = 1;
}
}
}
// if(indicator == 0) break;
// }
if(!returnpeaks) {
return accessible;
} else {
// BEWARE: this will not work if several connected genotypes
// have the same fitness and are maxima
int rowsum = 0;
Rcpp::IntegerVector peaks;
for(int k = 0; k < ng; ++k) {
if(accessible(k) > 0) {
rowsum = std::accumulate(adm(k, Rcpp::_).begin(),
adm(k, Rcpp::_).end(), 0);
if(rowsum == 0) { // This genotype doesn't have children
peaks.push_back(k + 1); // k is index. But in R, WT is in pos 1
}
}
}
return peaks;
}
}
// [[Rcpp::export]]
Rcpp::IntegerVector accessibleGenotypes(Rcpp::IntegerMatrix y,
Rcpp::NumericVector f,
Rcpp::IntegerVector numMut, //
double th) {
return accessibleGenotypesPeaksLandscape(y, f, numMut, th, false);
}
// [[Rcpp::export]]
Rcpp::IntegerVector peaksLandscape(Rcpp::IntegerMatrix y,
Rcpp::NumericVector f,
Rcpp::IntegerVector numMut, //
double th) {
return accessibleGenotypesPeaksLandscape(y, f, numMut, th, true);
}
// // This would make it easier returning the actual accessible genotypes easily
// // preserving the fitness if needed
// // Not being used now
// // [[Rcpp::export]]
// IntegerVector acc_ge(Rcpp::IntegerMatrix y, Rcpp::NumericVector f,
// Rcpp::IntegerVector numMut,
// int ng, //it counts the wt
// double th) {
// IntegerMatrix adm(ng, ng);
// int numMutdiff = 0;
// for(int i = 0; i < (ng - 1); ++i) {
// // Candidates are always of larger mutation by 1
// for(int j = (i + 1); j < ng; ++j) {
// numMutdiff = numMut(j) - numMut(i);
// if(numMutdiff > 1) { // no more to search
// break;
// } else if(numMutdiff == 1) {
// if( ( (f(j) - f(i)) >= th) &&
// (HammingDistance(y(i, _), y(j, _)) == 1) ) {
// adm(i, j) = 1;
// }
// }
// }
// }
// // Keeps root in Rows
// IntegerMatrix admtmp = adm(Range(0, ng - 1), Range(1, ng - 1));
// // Slightly different logic from R: Do not resize object; set the row to
// // 0.
// int colsum = 0;
// int indicator = 0; // indicator != 0 means we set one row to 0
// // so we need to iterate at least once more.
// // accessible is the genotype number, not the column! WT is 1,
// // etc. This makes it easy to keep track of which are accessible.
// IntegerVector accessible = seq_len(ng - 1) + 1;
// while (true) {
// indicator = 0;
// for(int k = 0; k < (ng - 1); ++k) {
// if(accessible(k) > 0) {
// colsum = std::accumulate(admtmp(_, k).begin(),
// admtmp(_, k).end(), 0);
// if(colsum == 0) { // This genotype ain't reachable
// // Recall row keeps Root.
// // Nothing can be reached from this genotype; fill with 0.
// admtmp(k + 1, _) = IntegerVector(ng - 1);
// accessible(k) = -9;
// indicator = 1;
// }
// }
// }
// if(indicator == 0) break;
// }
// return accessible;
// }
// // [[Rcpp::export]]
// Rcpp::IntegerVector accessibleGenotypes(Rcpp::IntegerMatrix y,
// Rcpp::NumericVector f,
// Rcpp::IntegerVector numMut, //
// double th) {
// // Return just the indices. Could preserve fitness, but would need
// // another matrix.
// int ng = y.nrow(); //it counts the wt
// Rcpp::IntegerMatrix adm(ng, ng);
// adm = integerAdjMat(y, f, numMut, th);
// int numMutdiff = 0;
// // Slightly different logic from R: Do not resize object; set the row to
// // 0.
// int colsum = 0;
// // int indicator = 0; // indicator != 0 means we set one row to 0
// // so we need to iterate at least once more.
// // accessible is the genotype number, not the column! WT is 1,
// // etc. This makes it easy to keep track of which are accessible.
// Rcpp::IntegerVector accessible = Rcpp::seq_len(ng);
// // This is doable in one pass
// // while (true) {
// // indicator = 0;
// for(int k = 1; k < ng; ++k) {
// if(accessible(k) > 0) {
// colsum = std::accumulate(adm(Rcpp::_, k).begin(),
// adm(Rcpp::_, k).end(), 0);
// if(colsum == 0) { // This genotype ain't reachable
// // Nothing can be reached from this genotype; fill with 0.
// adm(k, Rcpp::_) = Rcpp::IntegerVector(ng);
// accessible(k) = -9;
// // indicator = 1;
// }
// }
// }
// // if(indicator == 0) break;
// // }
// return accessible;
// }
// // [[Rcpp::export]]
// Rcpp::IntegerVector peaksLandscape(Rcpp::IntegerMatrix y,
// Rcpp::NumericVector f,
// Rcpp::IntegerVector numMut, //
// double th) {
// // Return the indices. This is like accessibleGenotypes, but we do an
// // extra loop
// int ng = y.nrow(); //it counts the wt
// Rcpp::IntegerMatrix adm(ng, ng);
// adm = integerAdjMat(y, f, numMut, th);
// int numMutdiff = 0;
// // Slightly different logic from R: Do not resize object; set the row to
// // 0.
// int colsum = 0;
// // int indicator = 0; // indicator != 0 means we set one row to 0
// // so we need to iterate at least once more.
// // accessible is the genotype number, not the column! WT is 1,
// // etc. This makes it easy to keep track of which are accessible.
// Rcpp::IntegerVector accessible = Rcpp::seq_len(ng);
// // This is doable in one pass
// // while (true) {
// // indicator = 0;
// for(int k = 1; k < ng; ++k) {
// if(accessible(k) > 0) {
// colsum = std::accumulate(adm(Rcpp::_, k).begin(),
// adm(Rcpp::_, k).end(), 0);
// if(colsum == 0) { // This genotype ain't reachable
// // Nothing can be reached from this genotype; fill with 0.
// adm(k, Rcpp::_) = Rcpp::IntegerVector(ng);
// accessible(k) = -9;
// // indicator = 1;
// }
// }
// }
// // if(indicator == 0) break;
// // }
// int rowsum = 0;
// Rcpp::IntegerVector peaks;
// for(int k = 1; k < ng; ++k) {
// if(accessible(k) > 0) {
// rowsum = std::accumulate(adm(k, Rcpp::_).begin(),
// adm(k, Rcpp::_).end(), 0);
// if(rowsum == 0) { // This genotype doesn't have children
// peaks.push_back(k);
// }
// }
// }
// return peaks;
// }
