Skip to content

Commit

Permalink
Made xi equivalent to other rates, set docs to use markdown
Browse files Browse the repository at this point in the history
  • Loading branch information
lucasnell committed Jun 11, 2018
1 parent 96f3ac1 commit b37e018
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 12 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -16,3 +16,4 @@ man
gemino_*.tgz
*.pdf
methods
_crash_test.R
1 change: 1 addition & 0 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@ Imports:
LinkingTo: Rcpp, RcppArmadillo, RcppProgress
SystemRequirements: C++11, GNU make
RoxygenNote: 6.0.1
Roxygen: list(markdown = TRUE)
Suggests:
knitr,
vcfR,
Expand Down
33 changes: 21 additions & 12 deletions src/mevo_rate_matrices.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ arma::mat TN93_rate_matrix(const std::vector<double>& pi_tcag,
// Filling in diagonals
Q.diag().fill(0.0); // reset to zero so summing by row works
arma::vec rowsums = arma::sum(Q, 1);
rowsums += xi;
rowsums += xi * 0.25;
rowsums *= -1;
Q.diag() = rowsums;

Expand Down Expand Up @@ -157,7 +157,7 @@ arma::mat GTR_rate_matrix(const std::vector<double>& pi_tcag,

// Filling in diagonals
arma::vec rowsums = arma::sum(Q, 1);
rowsums += xi;
rowsums += xi * 0.25;
rowsums *= -1;
Q.diag() = rowsums;

Expand All @@ -170,13 +170,18 @@ arma::mat GTR_rate_matrix(const std::vector<double>& pi_tcag,



/*
Estimates equilibrium nucleotide frequencies from an input rate matrix.
It does this by solving for πQ = 0 by finding the left eigenvector of Q that corresponds
to the eigenvalue closest to zero.
This is only needed for the UNREST model.
*/

//' Estimates equilibrium nucleotide frequencies from an input rate matrix.
//'
//' It does this by solving for πQ = 0 by finding the left eigenvector of Q that
//' corresponds to the eigenvalue closest to zero.
//' This is only needed for the UNREST model.
//'
//' @inheritParams Q UNREST_rate_matrix
//' @inheritParams pi_tcag UNREST_rate_matrix
//'
//' @noRd
//'
inline void est_pi_tcag(const arma::mat& Q, std::vector<double>& pi_tcag) {

arma::cx_vec eigvals;
Expand Down Expand Up @@ -210,8 +215,8 @@ inline void est_pi_tcag(const arma::mat& Q, std::vector<double>& pi_tcag) {
//' creating the matrix.
//'
//'
//' @param Q Matrix of rates for "T", "C", "A", and "G", respectively.
//' Diagonal values are ignored.
//' @param Q Matrix of substitution rates for "T", "C", "A", and "G", respectively.
//' Do not include indel rates here! Diagonal values are ignored.
//' @param pi_tcag Empty vector of equilibrium frequencies for for "T", "C", "A", and "G",
//' respectively. This vector will be filled in by this function.
//' @param xi Overall rate of indels.
Expand All @@ -232,8 +237,12 @@ void UNREST_rate_matrix(arma::mat& Q, std::vector<double>& pi_tcag, const double
// Estimate pi_tcag before incorporating indels:
est_pi_tcag(Q, pi_tcag);

// Now include indel rates
Q.diag() -= xi;
/*
Now include indel rates
(I'm subtracting bc diagonals are negative)
(Also, I'm doing this here so it doesn't affect pi_tcag estimates)
*/
Q.diag() -= xi * 0.25;

return;
}
Expand Down

0 comments on commit b37e018

Please sign in to comment.