diff --git a/alignment/alignmentpairwise.cpp b/alignment/alignmentpairwise.cpp index 8b39bfeac..0fb8b8196 100644 --- a/alignment/alignmentpairwise.cpp +++ b/alignment/alignmentpairwise.cpp @@ -18,7 +18,9 @@ * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * ***************************************************************************/ #include "alignmentpairwise.h" -#include "tree/phylosupertree.h" + +#include "tree/phylotree.h" +#include AlignmentPairwise::AlignmentPairwise() : Alignment(), Optimization() @@ -31,12 +33,11 @@ AlignmentPairwise::AlignmentPairwise() STATE_UNKNOWN = 0; trans_size = 0; trans_mat = nullptr; - sum_trans_mat = nullptr; trans_derv1 = nullptr; trans_derv2 = nullptr; + sum_trans = nullptr; sum_derv1 = nullptr; sum_derv2 = nullptr; - sum_trans = nullptr; pairCount = 0; derivativeCalculationCount = 0; costCalculationCount = 0; @@ -58,22 +59,21 @@ void AlignmentPairwise::setTree(PhyloTree* atree) { auto model = tree->getModel(); bool isModelSiteSpecific = (model==nullptr) ? false: model->isSiteSpecificModel(); if (model!=nullptr) { - trans_size = model->getTransMatrixSize(); + trans_size = num_states_squared; } + total_size = num_states_squared; if (!isModelSiteSpecific && !isRateSiteSpecific && rate!=nullptr && rate->getPtnCat(0) >= 0) { total_size *= rate->getNDiscreteRate(); } + pair_freq = new double[total_size]; trans_mat = new double[trans_size]; - sum_trans_mat = new double[trans_size]; + trans_derv1 = new double[trans_size]; + trans_derv2 = new double[trans_size]; sum_trans = new double[trans_size]; sum_derv1 = new double[trans_size]; sum_derv2 = new double[trans_size]; - trans_derv1 = new double[trans_size]; - trans_derv2 = new double[trans_size]; - total_size = num_states_squared; - pair_freq = new double[total_size]; - + pairCount = 0; derivativeCalculationCount = 0; costCalculationCount = 0; @@ -87,58 +87,47 @@ void AlignmentPairwise::setSequenceNumbers(int seq1, int seq2) { ++pairCount; seq_id1 = seq1; seq_id2 = seq2; - auto rate = tree->getRate(); - bool isRateSiteSpecific = (rate==nullptr) ? false : rate->isSiteSpecificRate(); - auto model = tree->getModel(); - bool isModelSiteSpecific = (model==nullptr) ? false: model->isSiteSpecificModel(); - if (isRateSiteSpecific || isModelSiteSpecific) { + ModelSubst *model = tree->getModel(); + RateHeterogeneity *site_rate = tree->getRate(); + size_t nptn = tree->aln->getNPattern(); + bool isModelSiteSpecific = (model) ? model->isSiteSpecificModel() : false; + bool isRateSiteSpecific = (site_rate) ? site_rate->isSiteSpecificRate() : false; + bool isRateCategorized = (site_rate) ? (site_rate->getPtnCat(0) >= 0) : false; + if (isModelSiteSpecific || isRateSiteSpecific) { return; } - memset(pair_freq, 0, sizeof(double)*total_size); - if (tree->hasMatrixOfConvertedSequences() - && rate->getPtnCat(0) < 0 ) { - auto sequence1 = tree->getConvertedSequenceByNumber(seq1); - auto sequence2 = tree->getConvertedSequenceByNumber(seq2); - auto frequencies = tree->getConvertedSequenceFrequencies(); + std::fill_n(pair_freq, total_size, 0.0); + if (tree->hasMatrixOfConvertedSequences() && !isRateCategorized) { + const char *sequence1 = tree->getConvertedSequenceByNumber(seq1); + const char *sequence2 = tree->getConvertedSequenceByNumber(seq2); + const int *frequencies = tree->getConvertedSequenceFrequencies(); size_t sequenceLength = tree->getConvertedSequenceLength(); - for (size_t i=0; i= num_states || state2 >= num_states) { continue; } - if ( state1 != STATE_UNKNOWN && state2 != STATE_UNKNOWN ) { + double *pairRow = pair_freq + state1*num_states; + if (state1 != STATE_UNKNOWN && state2 != STATE_UNKNOWN) { pairRow[state2] += frequencies[i]; } } - //Add back the cumulative frequencies for any sites - //that have the same state in every sequence. - for (int state=0; stategetSumOfFrequenciesForSitesWithConstantState(state); } - //Todo: Handle the multiple category case here - return; - } else if (tree->getRate()->getPtnCat(0) >= 0) { - int i = 0; - for (auto it = tree->aln->begin(); it != tree->aln->end(); it++, i++) { - int state1 = tree->aln->convertPomoState((*it)[seq_id1]); - int state2 = tree->aln->convertPomoState((*it)[seq_id2]); - addPattern(state1, state2, it->frequency, rate->getPtnCat(i)); - } - return; - } else { - for (auto it = tree->aln->begin(); it != tree->aln->end(); it++) { - int state1 = tree->aln->convertPomoState((*it)[seq_id1]); - int state2 = tree->aln->convertPomoState((*it)[seq_id2]); - addPattern(state1, state2, it->frequency); - } return; } + for (size_t ptn = 0; ptn < nptn; ++ptn) { + const Pattern &pat = tree->aln->at(ptn); + int state1 = tree->aln->convertPomoState(pat[seq_id1]); + int state2 = tree->aln->convertPomoState(pat[seq_id2]); + double freq = double(pat.frequency); + addPattern(state1, state2, freq, site_rate->getPtnCat(ptn)); + } } AlignmentPairwise::AlignmentPairwise(PhyloTree *atree, int seq1, int seq2) @@ -191,263 +180,198 @@ bool AlignmentPairwise::addPattern(int state1, int state2, int freq, int cat) { double AlignmentPairwise::computeFunction(double value) { ++costCalculationCount; - RateHeterogeneity *site_rate = tree->getRate(); - int ncat = site_rate->getNDiscreteRate(); - ModelSubst *model = tree->getModel(); - int nptn = tree->aln->getNPattern(); - double lh = 0.0; - - if (tree->hasMatrixOfConvertedSequences()) { - auto sequence1 = tree->getConvertedSequenceByNumber(seq_id1); - auto sequence2 = tree->getConvertedSequenceByNumber(seq_id2); - auto frequencies = tree->getConvertedSequenceFrequencies(); - size_t sequenceLength = tree->getConvertedSequenceLength(); - - if (site_rate->isSiteSpecificRate()) { - for (int i = 0; i < sequenceLength; i++) { - int state1 = sequence1[i]; - int state2 = sequence2[i]; - if (state1 >= num_states || state2 >= num_states) { - continue; - } - double trans = tree->getModelFactory()->computeTrans(value * site_rate->getPtnRate(i), state1, state2); - lh -= log(trans) * frequencies[i]; - } - return lh; - } else if (tree->getModel()->isSiteSpecificModel()) { - for (int i = 0; i < nptn; i++) { - int state1 = sequence1[i]; - int state2 = sequence2[i]; - if (state1 >= num_states || state2 >= num_states) { - continue; - } - double trans = tree->getModelFactory()->computeTrans(value * site_rate->getPtnRate(i), state1, state2); - lh -= log(trans) * frequencies[i]; - } - return lh; - } - } - // site-specific rates - if (site_rate->isSiteSpecificRate()) { - for (int i = 0; i < nptn; i++) { - int state1 = tree->aln->at(i)[seq_id1]; - int state2 = tree->aln->at(i)[seq_id2]; - if (state1 >= num_states || state2 >= num_states) continue; - double trans = tree->getModelFactory()->computeTrans(value * site_rate->getPtnRate(i), state1, state2); - lh -= log(trans) * tree->aln->at(i).frequency; - } - return lh; - } - if (tree->getModel()->isSiteSpecificModel()) { - for (int i = 0; i < nptn; i++) { - int state1 = tree->aln->at(i)[seq_id1]; - int state2 = tree->aln->at(i)[seq_id2]; - if (state1 >= num_states || state2 >= num_states) continue; - double trans = tree->getModel()->computeTrans(value, model->getPtnModelID(i), state1, state2); - lh -= log(trans) * tree->aln->at(i).frequency; - } - return lh; - } - - // categorized rates - if (site_rate->getPtnCat(0) >= 0) { - for (int cat = 0; cat < ncat; cat++) { - tree->getModelFactory()->computeTransMatrix(value*site_rate->getRate(cat), trans_mat); - double *pair_pos = pair_freq + cat*trans_size; - for (int i = 0; i < trans_size; i++) - if (pair_pos[i] > Params::getInstance().min_branch_length) { - if (trans_mat[i] <= 0) { - throw "Negative transition probability"; - } - lh -= pair_pos[i] * log(trans_mat[i]); - } - } - return lh; - } - - if (tree->getModelFactory()->site_rate->getGammaShape() == 0.0) - tree->getModelFactory()->computeTransMatrix(value, sum_trans_mat); - else { - tree->getModelFactory()->computeTransMatrix(value * site_rate->getRate(0), sum_trans_mat); - for (int cat = 1; cat < ncat; cat++) { - tree->getModelFactory()->computeTransMatrix(value * site_rate->getRate(cat), trans_mat); - for (int i = 0; i < trans_size; i++) - sum_trans_mat[i] += trans_mat[i]; - } - } - for (int i = 0; i < trans_size; i++) { - lh -= pair_freq[i] * log(sum_trans_mat[i]); - } - // negative log-likelihood (for minimization) - return lh; + double lh = 0.0, df = 0.0, ddf = 0.0; + likelihoodKernelFunction(value, lh, df, ddf); + return -lh; } void AlignmentPairwise::computeFuncDerv(double value, double &df, double &ddf) { ++derivativeCalculationCount; - RateHeterogeneity *site_rate = tree->getRate(); - int ncat = site_rate->getNDiscreteRate(); + double lh = 0.0; + likelihoodKernelFunction(value, lh, df, ddf); + df = -df; + ddf = -ddf; +} + +template +void AlignmentPairwise::likelihoodKernelFunction(double value, double &lh, double &df, double &ddf) { ModelSubst *model = tree->getModel(); - int trans_size = tree->getModel()->getTransMatrixSize(); - int nptn = tree->aln->getNPattern(); - df = 0.0; - ddf = 0.0; - - auto sequence1 = tree->getConvertedSequenceByNumber(seq_id1); - auto sequence2 = tree->getConvertedSequenceByNumber(seq_id2); - auto frequencies = tree->getConvertedSequenceFrequencies(); + RateHeterogeneity *site_rate = tree->getRate(); + ModelFactory *model_factory = tree->getModelFactory(); + size_t nptn = tree->aln->getNPattern(); + size_t ncat = site_rate->getNRate(); // # rate categories + size_t mcat = (model_factory->fused_mix_rate) ? 1 : ncat; // # rate categories per mixture class + size_t nmix = model->getNMixtures(); // # mixture classes + size_t ncat_mix = mcat * nmix; // # rate-mixture categories + lh = df = ddf = 0.0; + const double MIN_FREQ = Params::getInstance().min_branch_length; + const char *sequence1 = tree->getConvertedSequenceByNumber(seq_id1); + const char *sequence2 = tree->getConvertedSequenceByNumber(seq_id2); + const int *frequencies = tree->getConvertedSequenceFrequencies(); size_t sequenceLength = tree->getConvertedSequenceLength(); - if (sequenceLength!=nptn) { - sequence1 = sequence2 = nullptr; - frequencies = nullptr; - } - - if (site_rate->isSiteSpecificRate()) { - if (sequence1!=nullptr && sequence2!=nullptr && frequencies!=nullptr) { - #pragma omp parallel for reduction(-:df,ddf) schedule(dynamic,100) - for (int i = 0; i < nptn; ++i) { - int state1 = sequence1[i]; - if (num_states<=state1) { - continue; - } - int state2 = sequence2[i]; - if (num_states<=state2) { - continue; - } - double freq = frequencies[i]; - double rate_val = site_rate->getPtnRate(i); - double rate_sqr = rate_val * rate_val; - double derv1, derv2; - double trans = tree->getModelFactory()->computeTrans(value * rate_val, state1, state2, derv1, derv2); - double d1 = derv1 / trans; - df -= rate_val * d1 * freq; - ddf -= rate_sqr * (derv2/trans - d1*d1) * freq; - } + bool use_converted = tree->hasMatrixOfConvertedSequences() && (sequenceLength == nptn); + auto getPtnStatesAndFreq = + [use_converted, sequence1, sequence2, frequencies, this](size_t ptn, int &state1, int &state2, double &freq) { + if (use_converted) { + state1 = sequence1[ptn]; + state2 = sequence2[ptn]; + freq = double(frequencies[ptn]); } else { - for (int i = 0; i < nptn; i++) { - int state1 = tree->aln->at(i)[seq_id1]; - if (num_states<=state1) { - continue; - } - int state2 = tree->aln->at(i)[seq_id2]; - if (num_states<=state2) { - continue; - } - double rate_val = site_rate->getPtnRate(i); - double rate_sqr = rate_val * rate_val; - double derv1, derv2; - double trans = tree->getModelFactory()->computeTrans(value * rate_val, state1, state2, derv1, derv2); - double d1 = derv1 / trans; - double freq = tree->aln->at(i).frequency; - df -= rate_val * d1 * freq; - ddf -= rate_sqr * (derv2/trans - d1*d1) * freq; - } + const Pattern &pat = tree->aln->at(ptn); + state1 = pat[seq_id1]; + state2 = pat[seq_id2]; + freq = double(pat.frequency); } - return; - } - - if (tree->getModel()->isSiteSpecificModel()) { - if (sequence1!=nullptr && sequence2!=nullptr && frequencies!=nullptr) { - #pragma omp parallel for reduction(-:df,ddf) schedule(dynamic,100) - for (int i = 0; i < nptn; i++) { - int state1 = sequence1[i]; - if (num_states<=state1) { - continue; - } - int state2 = sequence2[i]; - if (num_states<=state2) { - continue; - } - double freq = frequencies[i]; - double rate_val = site_rate->getPtnRate(i); - double rate_sqr = rate_val * rate_val; - double derv1, derv2; - double trans = tree->getModel()->computeTrans(value * rate_val,model->getPtnModelID(i), state1, state2, derv1, derv2); - double d1 = derv1 / trans; - df -= rate_val * d1 * freq; - ddf -= rate_sqr * (derv2/trans - d1*d1) * freq; + }; + // site-specific model or rates + // Covers all relevant combinations: + // - site-specific model + site-specific/categorized/usual rates + // - usual model + site-specific rates + if (model->isSiteSpecificModel() || site_rate->isSiteSpecificRate()) { +#ifdef _OPENMP +#pragma omp parallel for reduction(+:lh,df,ddf) schedule(dynamic,100) +#endif + for (size_t ptn = 0; ptn < nptn; ++ptn) { + double freq; + int state1, state2; + getPtnStatesAndFreq(ptn, state1, state2, freq); + if (state1 >= num_states || state2 >= num_states) { + continue; } - } else { - for (int i = 0; i < nptn; i++) { - int state1 = tree->aln->at(i)[seq_id1]; - if (num_states<=state1) { - continue; + double lh_ptn = 0.0, df_ptn = 0.0, ddf_ptn = 0.0; + int model_id = model->getPtnModelID(ptn); + double rate = site_rate->getPtnRate(ptn); + for (size_t cm = 0; cm < ncat_mix; ++cm) { + size_t m = cm/mcat; + size_t c = cm%ncat; + if (nmix > 1) { + model_id = m; } - int state2 = tree->aln->at(i)[seq_id2]; - if (num_states<=state2) { - continue; + if (ncat > 1) { + rate = site_rate->getRate(c); } - double rate_val = site_rate->getPtnRate(i); - double rate_sqr = rate_val * rate_val; - double derv1, derv2; - double trans = tree->getModel()->computeTrans(value * rate_val,model->getPtnModelID(i), state1, state2, derv1, derv2); - double d1 = derv1 / trans; - double freq = tree->aln->at(i).frequency; - df -= rate_val * d1 * freq; - ddf -= rate_sqr * (derv2/trans - d1*d1) * freq; + double prop = site_rate->getProp(c) * model->getMixtureWeight(m); + if (!COMPUTE_DERV) { + double trans = model_factory->computeTrans(value * rate, state1, state2, model_id); + lh_ptn += trans * prop; + } else { + double prop_rate = prop * rate; + double prop_rate2 = prop_rate * rate; + double derv1, derv2; + double trans = model_factory->computeTrans(value * rate, state1, state2, derv1, derv2, model_id); + lh_ptn += trans * prop; + df_ptn += derv1 * prop_rate; + ddf_ptn += derv2 * prop_rate2; + } + } + if (state1 == state2) { + lh_ptn += site_rate->getPInvar(); + } + if (!COMPUTE_DERV) { + lh += log(lh_ptn) * freq; + } else { + // df = log(lh)' = lh'/lh + df_ptn /= lh_ptn; + df += df_ptn * freq; + // ddf = log(lh)'' = (lh'/lh)' = lh''/lh - (lh'/lh)^2 + ddf_ptn /= lh_ptn; + ddf_ptn -= df_ptn * df_ptn; + ddf += ddf_ptn * freq; } } return; } - - // categorized rates + // usual model and categorized rates if (site_rate->getPtnCat(0) >= 0) { - for (int cat = 0; cat < ncat; cat++) { - double rate_val = site_rate->getRate(cat); - double derv1 = 0.0, derv2 = 0.0; - tree->getModelFactory()->computeTransDerv(value*rate_val, trans_mat, trans_derv1, trans_derv2); + ASSERT(site_rate->getPInvar() == 0.0); + for (size_t cat = 0; cat < site_rate->getNDiscreteRate(); ++cat) { + std::fill_n(sum_trans, trans_size, 0.0); + std::fill_n(sum_derv1, trans_size, 0.0); + std::fill_n(sum_derv2, trans_size, 0.0); + double rate = site_rate->getRate(cat); + for (size_t m = 0; m < nmix; ++m) { + double prop = model->getMixtureWeight(m); + if (!COMPUTE_DERV) { + model_factory->computeTransMatrix(value * rate, trans_mat, m); + for (int i = 0; i < trans_size; ++i) { + sum_trans[i] += trans_mat[i] * prop; + } + } else { + double prop_rate = prop * rate; + double prop_rate2 = prop_rate * rate; + model_factory->computeTransDerv(value * rate, trans_mat, trans_derv1, trans_derv2, m); + for (int i = 0; i < trans_size; ++i) { + sum_trans[i] += trans_mat[i] * prop; + sum_derv1[i] += trans_derv1[i] * prop_rate; + sum_derv2[i] += trans_derv2[i] * prop_rate2; + } + } + } double *pair_pos = pair_freq + cat*trans_size; - for (int i = 0; i < trans_size; i++) if (pair_pos[i] > 0) { - if (trans_mat[i] <= 0) { - throw "Negative transition probability"; + for (int i = 0; i < trans_size; ++i) { + if (pair_pos[i] > MIN_FREQ) { + ASSERT(sum_trans[i] > 0.0); + if (!COMPUTE_DERV) { + lh += log(sum_trans[i]) * pair_pos[i]; + } else { + // df = log(lh)' = lh'/lh + double df_pair = sum_derv1[i] / sum_trans[i]; + df += df_pair * pair_pos[i]; + // ddf = log(lh)'' = (lh'/lh)' = lh''/lh - (lh'/lh)^2 + double ddf_pair = sum_derv2[i] / sum_trans[i] - df_pair * df_pair; + ddf += ddf_pair * pair_pos[i]; + } } - double d1 = trans_derv1[i] / trans_mat[i]; - derv1 += pair_pos[i] * d1; - derv2 += pair_pos[i] * (trans_derv2[i]/trans_mat[i] - d1 * d1); } - df -= derv1 * rate_val; - ddf -= derv2 * rate_val * rate_val; } return; } - - memset(sum_trans, 0, sizeof(double) * trans_size); - memset(sum_derv1, 0, sizeof(double) * trans_size); - memset(sum_derv2, 0, sizeof(double) * trans_size); - - for (int cat = 0; cat < ncat; cat++) { - double rate_val = site_rate->getRate(cat); - double prop_val = site_rate->getProp(cat); - if (tree->getModelFactory()->site_rate->getGammaShape() == 0.0) - { - rate_val = 1.0; - } - double coeff1 = rate_val * prop_val; - double coeff2 = rate_val * coeff1; - //cout << "cat " << cat << "," << (intptr_t)trans_mat << ", " << (intptr_t)trans_derv1 << ", " << (intptr_t)trans_derv2 << endl; - tree->getModelFactory()->computeTransDerv(value * rate_val, trans_mat, trans_derv1, trans_derv2); - for (int i = 0; i < trans_size; i++) { - sum_trans[i] += trans_mat[i] * prop_val; - sum_derv1[i] += trans_derv1[i] * coeff1; - sum_derv2[i] += trans_derv2[i] * coeff2; + // usual model and rates + std::fill_n(sum_trans, trans_size, 0.0); + std::fill_n(sum_derv1, trans_size, 0.0); + std::fill_n(sum_derv2, trans_size, 0.0); + for (size_t cm = 0; cm < ncat_mix; ++cm) { + size_t m = cm/mcat; + size_t c = cm%ncat; + double rate = site_rate->getRate(c); + double prop = site_rate->getProp(c) * model->getMixtureWeight(m); + if (!COMPUTE_DERV) { + model_factory->computeTransMatrix(value * rate, trans_mat, m); + for (int i = 0; i < trans_size; ++i) { + sum_trans[i] += trans_mat[i] * prop; + } + } else { + double prop_rate = prop * rate; + double prop_rate2 = prop_rate * rate; + model_factory->computeTransDerv(value * rate, trans_mat, trans_derv1, trans_derv2, m); + for (int i = 0; i < trans_size; ++i) { + sum_trans[i] += trans_mat[i] * prop; + sum_derv1[i] += trans_derv1[i] * prop_rate; + sum_derv2[i] += trans_derv2[i] * prop_rate2; + } } } - - // 2019-07-03: incorporate p_invar double p_invar = site_rate->getPInvar(); if (p_invar > 0.0) { - for (int i = 0; i < num_states; i++) { - sum_trans[i*num_states+i] += p_invar; + for (int x = 0; x < num_states; ++x) { + sum_trans[x*num_states+x] += p_invar; } } - - for (int i = 0; i < trans_size; i++) { - if (pair_freq[i] > Params::getInstance().min_branch_length && sum_trans[i] > 0.0) { - double d1 = sum_derv1[i] / sum_trans[i]; - df -= pair_freq[i] * d1; - ddf -= pair_freq[i] * (sum_derv2[i]/sum_trans[i] - d1 * d1); + for (int i = 0; i < trans_size; ++i) { + if (pair_freq[i] > MIN_FREQ) { + ASSERT(sum_trans[i] > 0.0); + if (!COMPUTE_DERV) { + lh += log(sum_trans[i]) * pair_freq[i]; + } else { + // df = log(lh)' = lh'/lh + double df_pair = sum_derv1[i] / sum_trans[i]; + df += df_pair * pair_freq[i]; + // ddf = log(lh)'' = (lh'/lh)' = lh''/lh - (lh'/lh)^2 + double ddf_pair = sum_derv2[i] / sum_trans[i] - df_pair * df_pair; + ddf += ddf_pair * pair_freq[i]; + } } } - return; } double AlignmentPairwise::optimizeDist(double initial_dist, double &d2l) { @@ -529,7 +453,6 @@ AlignmentPairwise::~AlignmentPairwise() delete [] sum_trans; delete [] trans_derv2; delete [] trans_derv1; - delete [] sum_trans_mat; delete [] trans_mat; delete [] pair_freq; } diff --git a/alignment/alignmentpairwise.h b/alignment/alignmentpairwise.h index cae1c75b0..d922fc095 100644 --- a/alignment/alignmentpairwise.h +++ b/alignment/alignmentpairwise.h @@ -112,10 +112,22 @@ class AlignmentPairwise : public Alignment, public Optimization */ virtual ~AlignmentPairwise(); +private: + /** + * Implementation of both computeFunction() and computeFuncDerv() + * @param value A distance between the two sequences + * @param[out] lh The log-likelihood at this distance + * @param[out] df The log-likelihood first derivative at this distance + * @param[out] ddf The log-likelihood second derivative at this distance + */ + template + void likelihoodKernelFunction(double value, double &lh, double &df, double &ddf); + +public: size_t pairCount; size_t derivativeCalculationCount; size_t costCalculationCount; - + protected: PhyloTree* tree; //multi-species alignment tree from which sequences //to be aligned are to be drawn @@ -125,13 +137,12 @@ class AlignmentPairwise : public Alignment, public Optimization //size is num_states_squared times 1 (or by the number //of categories). int trans_size; //number of elements (rows x columns) in transition matrices - double* trans_mat; //used in computeFunction(), - double* sum_trans_mat; //used in computeFunction() + double* trans_mat; //used in computeFunction() and computeFuncDerv() double* trans_derv1; //used in computeFuncDerv() double* trans_derv2; //used in computeFuncDerv() + double* sum_trans; //used in computeFunction() and computeFuncDerv() double* sum_derv1; //used in computeFuncDerv() double* sum_derv2; //used in computeFuncDerv() - double* sum_trans; //used in computeFuncDerv() int seq_id1; int seq_id2; diff --git a/model/modelcodon.cpp b/model/modelcodon.cpp index 8c7df792e..59036671b 100644 --- a/model/modelcodon.cpp +++ b/model/modelcodon.cpp @@ -832,7 +832,7 @@ void ModelCodon::decomposeRateMatrix() { ModelMarkov::decomposeRateMatrix(); } -void ModelCodon::getQMatrix(double *q_mat, int mixture) { +void ModelCodon::getQMatrix(double *q_mat, int) { double **rate_matrix = (double**) new double[num_states]; int i, j; diff --git a/model/modelcodon.h b/model/modelcodon.h index 921dcbce7..67b318c6e 100644 --- a/model/modelcodon.h +++ b/model/modelcodon.h @@ -100,7 +100,7 @@ class ModelCodon: public ModelMarkov { * compute Q matrix * @param q_mat (OUT) Q matrix, assuming of size num_states * num_states */ - virtual void getQMatrix(double *q_mat, int mixture = 0); + virtual void getQMatrix(double *q_mat, int model_id = -1); /** * read codon model from a stream, modying rates and state_freq accordingly diff --git a/model/modelfactory.cpp b/model/modelfactory.cpp index a45f1a842..64836c7e5 100644 --- a/model/modelfactory.cpp +++ b/model/modelfactory.cpp @@ -1771,17 +1771,18 @@ void ModelFactory::stopStoringTransMatrix() { } -double ModelFactory::computeTrans(double time, int state1, int state2) { - return model->computeTrans(time, state1, state2); +double ModelFactory::computeTrans(double time, int state1, int state2, int model_id) { + return model->computeTrans(time, state1, state2, model_id); } -double ModelFactory::computeTrans(double time, int state1, int state2, double &derv1, double &derv2) { - return model->computeTrans(time, state1, state2, derv1, derv2); +double ModelFactory::computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id) { + return model->computeTrans(time, state1, state2, derv1, derv2, model_id); } -void ModelFactory::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { +void ModelFactory::computeTransMatrix(double time, double *trans_matrix, int model_id, int selected_row) { if (!store_trans_matrix || !is_storing || model->isSiteSpecificModel()) { - model->computeTransMatrix(time, trans_matrix, mixture, selected_row); + model->computeTransMatrix(time, trans_matrix, model_id, selected_row); return; } int mat_size = model->num_states * model->num_states; @@ -1790,20 +1791,19 @@ void ModelFactory::computeTransMatrix(double time, double *trans_matrix, int mix // allocate memory for 3 matricies double *trans_entry = new double[mat_size * 3]; trans_entry[mat_size] = trans_entry[mat_size+1] = 0.0; - model->computeTransMatrix(time, trans_entry, mixture, selected_row); + model->computeTransMatrix(time, trans_entry, model_id, selected_row); ass_it = insert(value_type(round(time * 1e6), trans_entry)).first; } else { //if (verbose_mode >= VB_MAX) //cout << "ModelFactory bingo" << endl; } - memcpy(trans_matrix, ass_it->second, mat_size * sizeof(double)); } void ModelFactory::computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture) { + double *trans_derv1, double *trans_derv2, int model_id) { if (!store_trans_matrix || !is_storing || model->isSiteSpecificModel()) { - model->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2, mixture); + model->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2, model_id); return; } int mat_size = model->num_states * model->num_states; @@ -1812,11 +1812,11 @@ void ModelFactory::computeTransDerv(double time, double *trans_matrix, // allocate memory for 3 matricies double *trans_entry = new double[mat_size * 3]; trans_entry[mat_size] = trans_entry[mat_size+1] = 0.0; - model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), mixture); + model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), model_id); ass_it = insert(value_type(round(time * 1e6), trans_entry)).first; } else if (ass_it->second[mat_size] == 0.0 && ass_it->second[mat_size+1] == 0.0) { double *trans_entry = ass_it->second; - model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), mixture); + model->computeTransDerv(time, trans_entry, trans_entry+mat_size, trans_entry+(mat_size*2), model_id); } memcpy(trans_matrix, ass_it->second, mat_size * sizeof(double)); memcpy(trans_derv1, ass_it->second + mat_size, mat_size * sizeof(double)); diff --git a/model/modelfactory.h b/model/modelfactory.h index 8b29cb062..022eda3fd 100644 --- a/model/modelfactory.h +++ b/model/modelfactory.h @@ -139,47 +139,29 @@ class ModelFactory : public unordered_map, public Optimization, pu */ void stopStoringTransMatrix(); - /** - Wrapper for computing the transition probability matrix from the model. It use ModelFactory - that stores matrix computed before for effiency purpose. - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); + /** + * Wrapper for the computeTransMatrix() function of the model. + * ModelFactory stores the matrix computed before for efficiency purposes + */ + void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); - /** - Wrapper for computing the transition probability between two states. - @param time time between two events - @param state1 first state - @param state2 second state - */ - double computeTrans(double time, int state1, int state2); + /** + * Wrapper for the computeTransDerv() function of the model. + * ModelFactory stores the matrix computed before for efficiency purposes + */ + void computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id = -1); - /** - Wrapper for computing the transition probability between two states - @param time time between two events - @param state1 first state - @param state2 second state - @param derv1 (OUT) 1st derivative - @param derv2 (OUT) 2nd derivative - */ - virtual double computeTrans(double time, int state1, int state2, double &derv1, double &derv2); + /** + * Wrapper for the computeTrans() function of the model + */ + double computeTrans(double time, int state1, int state2, int model_id = -1); - /** - Wrapper for computing the transition probability matrix and the derivative 1 and 2 from the model. - It use ModelFactory that stores matrix computed before for effiency purpose. - @param time time between two events - @param mixture (optional) class for mixture model - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - @param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. - @param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. - */ - void computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture = 0); + /** + * Wrapper for the computeTrans() function of the model + */ + double computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id = -1); /** destructor diff --git a/model/modelliemarkov.cpp b/model/modelliemarkov.cpp index 501cb0976..78d24c036 100644 --- a/model/modelliemarkov.cpp +++ b/model/modelliemarkov.cpp @@ -2269,9 +2269,11 @@ void ModelLieMarkov::decomposeRateMatrixClosedForm() { } } +// OBSOLETE +/* void ModelLieMarkov::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { return ModelMarkov::computeTransMatrix(time, trans_matrix, mixture, selected_row); - /* + MatrixExpTechnique technique = phylo_tree->params->matrix_exp_technique; if (technique == MET_SCALING_SQUARING || nondiagonalizable ) { Matrix4d A = Map(rate_matrix); @@ -2331,6 +2333,5 @@ void ModelLieMarkov::computeTransMatrix(double time, double *trans_matrix, int m } else ModelMarkov::computeTransMatrix(time, trans_matrix); - */ } - +*/ diff --git a/model/modelliemarkov.h b/model/modelliemarkov.h index 31c5b3083..d3c6cac75 100644 --- a/model/modelliemarkov.h +++ b/model/modelliemarkov.h @@ -95,15 +95,6 @@ class ModelLieMarkov: public ModelMarkov { /** decompose rate matrix using Eigen library */ virtual void decomposeRateMatrixEigen3lib(); - /** - compute the transition probability matrix. - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); // overrides Optimization::restartParameters bool restartParameters(double guess[], int ndim, double lower[], double upper[], bool bound_check[], int iteration); diff --git a/model/modelmarkov.cpp b/model/modelmarkov.cpp index c1cb4df0e..d5132e5b3 100644 --- a/model/modelmarkov.cpp +++ b/model/modelmarkov.cpp @@ -322,14 +322,14 @@ void ModelMarkov::init_state_freq(StateFreqType type) { if (phylo_tree->aln->seq_type == SEQ_CODON) { double ntfreq[12]; phylo_tree->aln->computeCodonFreq(freq_type, state_freq, ntfreq); -// phylo_tree->aln->computeCodonFreq(state_freq); } else if (phylo_tree->aln->seq_type != SEQ_POMO) { - double emp_state_freq[num_states]; - phylo_tree->aln->computeStateFreq(emp_state_freq); - setStateFrequency(emp_state_freq); - } for (i = 0; i < num_states; i++) - if (state_freq[i] > state_freq[highest_freq_state]) - highest_freq_state = i; + phylo_tree->aln->computeStateFreq(state_freq); + } + for (int x = 0; x < num_states; ++x) { + if (state_freq[x] > state_freq[highest_freq_state]) { + highest_freq_state = x; + } + } break; case FREQ_USER_DEFINED: { @@ -460,7 +460,7 @@ void ModelMarkov::report_state_freqs(ostream& out, double *custom_state_freq) { } } -void ModelMarkov::computeTransMatrixNonrev(double time, double *trans_matrix, int mixture) { +void ModelMarkov::computeTransMatrixNonrev(double time, double *trans_matrix) { auto technique = phylo_tree->params->matrix_exp_technique; if (technique == MET_SCALING_SQUARING || nondiagonalizable) { // scaling and squaring technique @@ -495,7 +495,7 @@ void ModelMarkov::computeTransMatrixNonrev(double time, double *trans_matrix, in cout << "INFO: Switch to scaling-squaring due to unstable eigen-decomposition rowsum: " << mincoeff << " to " << maxcoeff << endl; nondiagonalizable = true; - computeTransMatrixNonrev(time, trans_matrix, mixture); + computeTransMatrixNonrev(time, trans_matrix); nondiagonalizable = false; } } else { @@ -504,10 +504,10 @@ void ModelMarkov::computeTransMatrixNonrev(double time, double *trans_matrix, in } -void ModelMarkov::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { +void ModelMarkov::computeTransMatrix(double time, double *trans_matrix, int, int selected_row) { if (!is_reversible) { - computeTransMatrixNonrev(time, trans_matrix, mixture); + computeTransMatrixNonrev(time, trans_matrix); return; } @@ -567,40 +567,61 @@ void ModelMarkov::computeTransMatrix(double time, double *trans_matrix, int mixt // delete [] exptime; } -double ModelMarkov::computeTrans(double time, int state1, int state2) { - +double ModelMarkov::computeTrans(double time, int state1, int state2, int) { + double trans = 0.0; if (is_reversible) { + // reversible double evol_time = time / total_num_subst; - int i; - double trans_prob = 0.0; - for (i = 0; i < num_states; i++) { - trans_prob += eigenvectors[state1*num_states+i] * inv_eigenvectors[i*num_states+state2] * exp(evol_time * eigenvalues[i]); + for (int i = 0; i < num_states; ++i) { + double cof = eigenvalues[i]; + double lhval = eigenvectors[state1*num_states+i] + * inv_eigenvectors[i*num_states+state2] + * exp(evol_time * cof); + trans += lhval; } - return trans_prob; } else { // non-reversible + int addr = state1*num_states+state2; double *trans_matrix = new double[num_states*num_states]; computeTransMatrix(time, trans_matrix); - double trans = trans_matrix[state1*num_states+state2]; + trans = trans_matrix[addr]; delete [] trans_matrix; - return trans; } + return trans; } -double ModelMarkov::computeTrans(double time, int state1, int state2, double &derv1, double &derv2) { - double evol_time = time / total_num_subst; - double trans_prob = 0.0; +double ModelMarkov::computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int) { + double trans = 0.0; derv1 = derv2 = 0.0; - for (int i = 0; i < num_states; i++) { - double trans = eigenvectors[state1*num_states+i] - * inv_eigenvectors[i*num_states+state2] - * exp(evol_time * eigenvalues[i]); - double trans2 = trans * eigenvalues[i]; - trans_prob += trans; - derv1 += trans2; - derv2 += trans2 * eigenvalues[i]; + if (is_reversible) { + // reversible + double evol_time = time / total_num_subst; + for (int i = 0; i < num_states; ++i) { + double cof = eigenvalues[i]; + double lhval = eigenvectors[state1*num_states+i] + * inv_eigenvectors[i*num_states+state2] + * exp(evol_time * cof); + double dfval = cof*lhval; + trans += lhval; + derv1 += dfval; + derv2 += cof*dfval; + } + } else { + // non-reversible + int addr = state1*num_states+state2; + double *trans_matrix = new double[num_states*num_states]; + double *trans_derv1 = new double[num_states*num_states]; + double *trans_derv2 = new double[num_states*num_states]; + computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2); + trans = trans_matrix[addr]; + derv1 = trans_derv1[addr]; + derv2 = trans_derv2[addr]; + delete [] trans_matrix; + delete [] trans_derv1; + delete [] trans_derv2; } - return trans_prob; + return trans; } void ModelMarkov::calculateExponentOfScalarMultiply ( const double* source, int size @@ -729,9 +750,8 @@ void ModelMarkov::aTimesDiagonalBTimesTransposeOfC(const double* matrixA, const } } -void ModelMarkov::computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture) -{ +void ModelMarkov::computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int) { if (!is_reversible) { computeTransMatrix(time, trans_matrix); // First derivative = Q * e^(Qt) @@ -844,92 +864,80 @@ void ModelMarkov::computeTransDerv(double time, double *trans_matrix, // delete [] exptime; } -void ModelMarkov::getRateMatrix(double *rate_mat) { - int nrate = getNumRateEntries(); - memcpy(rate_mat, rates, nrate * sizeof(double)); +void ModelMarkov::getRateMatrix(double *rate_mat, int) { + ASSERT(rates); + int nrate = getNumRateEntries(); + std::copy_n(rates, nrate, rate_mat); } -void ModelMarkov::setRateMatrix(double* rate_mat) -{ - int nrate = getNumRateEntries(); - memcpy(rates, rate_mat, nrate * sizeof(double)); +void ModelMarkov::setRateMatrix(double *rate_mat) { + ASSERT(rates); + int nrate = getNumRateEntries(); + std::copy_n(rate_mat, nrate, rates); } -void ModelMarkov::setFullRateMatrix(double* rate_mat, double *freq) -{ - int i, j, k; +void ModelMarkov::setQMatrix(double *q_mat, double *freq_vec) { + // set rates if (isReversible()) { - for (i = 0, k = 0; i < num_states; i++) - for (j = i+1; j < num_states; j++) - rates[k++] = rate_mat[i*num_states+j] / freq[j]; - memcpy(state_freq, freq, sizeof(double)*num_states); + // reversible + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = i+1; j < num_states; ++j) { + rates[k++] = q_mat[i*num_states+j] / freq_vec[j]; + } + } } else { // non-reversible - for (i = 0, k = 0; i < num_states; i++) - for (j = 0; j < num_states; j++) - if (i != j) - rates[k++] = rate_mat[i*num_states+j]; + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = 0; j < num_states; ++j) { + if (i != j) { + rates[k++] = q_mat[i*num_states+j]; + } + } + } } + // set state_freq + setStateFrequency(freq_vec); } -void ModelMarkov::getStateFrequency(double *freq, int mixture) { - ASSERT(state_freq); - ASSERT(freq_type != FREQ_UNKNOWN); - memcpy(freq, state_freq, sizeof(double) * num_states); - // // DEBUG. - // cout << setprecision(8); - // cout << "State frequency reported by ModelMarkov: "; - // for (int i = 0; i < num_states; i++) { - // cout << state_freq[i] << " "; - // } - // cout << endl; +void ModelMarkov::getStateFrequency(double *freq_vec, int) { + ASSERT(state_freq); + ASSERT(freq_type != FREQ_UNKNOWN); + std::copy_n(state_freq, num_states, freq_vec); // 2015-09-07: relax the sum of state_freq to be 1, this will be done at the end of optimization double sum = 0.0; - int i; - for (i = 0; i < num_states; i++) sum += freq[i]; - sum = 1.0/sum; - for (i = 0; i < num_states; i++) freq[i] *= sum; + for (int x = 0; x < num_states; ++x) { + sum += freq_vec[x]; + } + for (int x = 0; x < num_states; ++x) { + freq_vec[x] /= sum; + } } -void ModelMarkov::setStateFrequency(double* freq) -{ - ASSERT(state_freq); - /* - if (!isReversible()) { - // integrate out state_freq from rate_matrix - int i, j, k = 0; - for (i = 0, k = 0; i < num_states; i++) - for (j = 0; j < num_states; j++) - if (i != j) { - rates[k] = (rates[k])*freq[j]; - if (state_freq[j] != 0.0) - rates[k] /= state_freq[j]; - k++; - } - } - */ - ModelSubst::setStateFrequency(freq); +void ModelMarkov::setStateFrequency(double *freq_vec) { + ASSERT(state_freq); + std::copy_n(freq_vec, num_states, state_freq); } -void ModelMarkov::adaptStateFrequency(double* freq) -{ +void ModelMarkov::adaptStateFrequency(double *freq_vec) { ASSERT(state_freq); if (!isReversible()) { - // integrate out state_freq from rate_matrix - int i, j, k = 0; - for (i = 0, k = 0; i < num_states; i++) - for (j = 0; j < num_states; j++) + // substitute state_freq with freq_vec in the Q matrix (rates) + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = 0; j < num_states; ++j) { if (i != j) { - rates[k] = (rates[k])*freq[j]; - if (state_freq[j] > ZERO_FREQ) + rates[k] *= freq_vec[j]; + if (state_freq[j] > ZERO_FREQ) { rates[k] /= state_freq[j]; + } k++; } + } + } } - ModelSubst::setStateFrequency(freq); + setStateFrequency(freq_vec); } -void ModelMarkov::getQMatrix(double *q_mat, int mixture) { +void ModelMarkov::getQMatrix(double *q_mat, int) { if (!is_reversible) { // non-reversible model @@ -1904,6 +1912,24 @@ void ModelMarkov::freeMem() internalFreeMem(); } +void ModelMarkov::multiplyWithInvEigenvector(double *state_lh) { + int nmix = getNMixtures(); + int nstates = get_safe_upper_limit(num_states); + double saved_state_lh[num_states]; + memcpy(saved_state_lh, state_lh, sizeof(double)*num_states); + memset(state_lh, 0, sizeof(double)*num_states*nmix); + const double *inv_eigenvectors = getInverseEigenvectors(); + for (int m = 0; m < nmix; ++m) { + double *this_state_lh = &state_lh[m*num_states]; + const double *inv_evec = &inv_eigenvectors[m * nstates * num_states]; + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = 0; j < num_states; ++j) { + this_state_lh[i] += inv_evec[k++] * saved_state_lh[j]; + } + } + } +} + double *ModelMarkov::getEigenvalues() const { return eigenvalues; diff --git a/model/modelmarkov.h b/model/modelmarkov.h index 0c93b28fa..2d28cdcba 100644 --- a/model/modelmarkov.h +++ b/model/modelmarkov.h @@ -178,86 +178,113 @@ class ModelMarkov : public ModelSubst, public EigenDecomposition */ void readParametersString(string &model_str, bool adapt_tree = true); - /** - compute the transition probability matrix. - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); + /** + * Compute the transition probability matrix on a branch. + * The default is the JC model, valid for all kinds of data + * @param time The branch length + * @param model_id ID of a submodel (for mixture and site-specific models) + * @param selected_row Only compute the entries for the selected row, + * the default is to compute entries for all rows + * @param[out] trans_matrix The transition matrix between all pairs of states, + * assumed to have the size of num_states*num_states + */ + virtual void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); + + /** Helper function of computeTransMatrix() for non-reversible models */ + void computeTransMatrixNonrev(double time, double *trans_matrix); /** - compute the transition probability matrix for non-reversible model - @param time time between two events - @param mixture (optional) class for mixture model - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. + * The same as computeTransMatrix() above, but also computes the + * 1st and 2nd derivative matrices with respect to the branch length + * @param[out] trans_matrix The transition matrix between all pairs of states, + * assumed to have the size of num_states*num_states + * @param[out] trans_derv1 The 1st derivative matrix between all pairs of states + * @param[out] trans_derv2 The 2nd derivative matrix between all pairs of states */ - virtual void computeTransMatrixNonrev(double time, double *trans_matrix, int mixture = 0); + virtual void computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id = -1); - /** - compute the transition probability between two states - @param time time between two events - @param state1 first state - @param state2 second state - */ - virtual double computeTrans(double time, int state1, int state2); + /** + * Compute the transition probability between the two states on a branch. + * The default is the JC model, valid for all kinds of data + * @param time The branch length between the two states + * @param model_id ID of a submodel (for mixture and site-specific models) + * @param state1 The start state + * @param state2 The end state + * @param[out] derv1 The 1st derivative + * @param[out] derv2 The 2nd derivative + * @return The transition probability + */ + virtual double computeTrans(double time, int state1, int state2, int model_id = -1); - /** - compute the transition probability between two states - @param time time between two events - @param state1 first state - @param state2 second state - @param derv1 (OUT) 1st derivative - @param derv2 (OUT) 2nd derivative - */ - virtual double computeTrans(double time, int state1, int state2, double &derv1, double &derv2); + /** + * The same as computeTrans() above, but also computes the + * 1st and 2nd derivatives with respect to the branch length + * @param[out] derv1 The 1st derivative + * @param[out] derv2 The 2nd derivative + * @return The transition probability + */ + virtual double computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id = -1); - /** - Get the rate matrix. - @param rate_mat (OUT) upper-triagle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 - */ - virtual void getRateMatrix(double *rate_mat); + /** + * Get the rate parameters, such as a,b,c,d,e,f for a DNA model. + * Get the above-diagonal entries of the rate matrix, assuming that + * the last element is 1. + * The default is equal rates of 1 (JC Model), valid for all kinds of data + * @param[out] rate_mat An upper-triangle rate matrix, assumed to have the + * size of num_states*(num_states-1)/2 + * @param model_id ID of a submodel (for mixture and site-specific models) + */ + virtual void getRateMatrix(double *rate_mat, int model_id = -1); - /** - Set the rate matrix. - @param rate_mat upper-triagle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 - */ - virtual void setRateMatrix(double *rate_mat); + /** + * Set the rate parameters + * @param rate_mat An upper-triangle rate matrix, assumed to have the + * size of num_states*(num_states-1)/2 + */ + virtual void setRateMatrix(double *rate_mat); /** - Set the full rate matrix of size num_states*num_states - @param rate_mat full rate matrix - @param freq state frequency + * Get the instantaneous rate matrix Q. + * The default is derived from equal rates and equal state frequencies + * @param[out] q_mat A full matrix: qij >= 0, qii = -sum_j qij (j != i), + * assumed to have the size of num_states*num_states + * @param model_id ID of a submodel (for mixture and site-specific models) */ - virtual void setFullRateMatrix(double *rate_mat, double *freq); + virtual void getQMatrix(double *q_mat, int model_id = -1); - /** - compute the state frequency vector - @param mixture (optional) class for mixture model - @param state_freq (OUT) state frequency vector. Assume state_freq has size of num_states - */ - virtual void getStateFrequency(double *state_freq, int mixture = 0); + /** + * Set the instantaneous rate matrix Q. + * @param q_mat A full matrix: qij >= 0, qii = -sum_j qij (j != i), + * assumed to have the size of num_states*num_states + * @param freq_vec Respective state frequency vector, assumed to have the + * size of num_states + */ + virtual void setQMatrix(double *q_mat, double *freq_vec); - /** - set the state frequency vector - @param state_freq (IN) state frequency vector. Assume state_freq has size of num_states - */ - virtual void setStateFrequency(double *state_freq); + /** + * Get the state frequency vector. + * The default is equal state frequencies, valid for all kinds of data + * @param[out] freq_vec A state frequency vector, assumed to have the + * size of num_states + * @param model_id ID of a submodel (for mixture and site-specific models) + */ + virtual void getStateFrequency(double *freq_vec, int model_id = -1); /** - set the state frequency vector - @param state_freq (IN) state frequency vector. Assume state_freq has size of num_states + * Set the state frequency vector + * @param freq_vec A state frequency vector, assumed to have the + * size of num_states */ - virtual void adaptStateFrequency(double *state_freq); + virtual void setStateFrequency(double *freq_vec); - /** - * compute Q matrix - * @param q_mat (OUT) Q matrix, assuming of size num_states * num_states - */ - virtual void getQMatrix(double *q_mat, int mixture = 0); + /** + * Set the state frequency vector and adjust the Q matrix accordingly + * @param freq_vec A state frequency vector, assumed to have the + * size of num_states + */ + virtual void adaptStateFrequency(double *freq_vec); /** rescale the state frequencies @@ -271,19 +298,6 @@ class ModelMarkov : public ModelSubst, public EigenDecomposition */ virtual StateFreqType getFreqType() { return freq_type; } - - /** - compute the transition probability matrix.and the derivative 1 and 2 - @param time time between two events - @param mixture (optional) class for mixture model - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - @param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. - @param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. - */ - virtual void computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture = 0); - /** @return the number of dimensions */ @@ -352,22 +366,21 @@ class ModelMarkov : public ModelSubst, public EigenDecomposition */ virtual void decomposeRateMatrix(); -// double *getEigenCoeff() const; - - virtual double *getEigenvalues() const; + /** + * For reversible models, multiply the partial likelihood vector with + * the matrix of inverse eigenvectors for the fast pruning algorithm + * @param[in/out] state_lh The partial likelihood vector + */ + virtual void multiplyWithInvEigenvector(double *state_lh); - virtual double *getEigenvectors() const; - virtual double *getInverseEigenvectors() const; + virtual double *getEigenvalues() const; + virtual double *getEigenvectors() const; + virtual double *getInverseEigenvectors() const; virtual double *getInverseEigenvectorsTransposed() const; -// void setEigenCoeff(double *eigenCoeff); - - void setEigenvalues(double *eigenvalues); - - void setEigenvectors(double *eigenvectors); - - void setInverseEigenvectors(double *inv_eigenvectors); - + void setEigenvalues(double *eigenvalues); + void setEigenvectors(double *eigenvectors); + void setInverseEigenvectors(double *inv_eigenvectors); void setInverseEigenvectorsTransposed(double *inv_eigenvectors); static void calculateExponentOfScalarMultiply(const double* source, int size diff --git a/model/modelmixture.cpp b/model/modelmixture.cpp index dab09e7d5..ea54af6dd 100644 --- a/model/modelmixture.cpp +++ b/model/modelmixture.cpp @@ -3803,35 +3803,6 @@ void ModelMixture::restoreCheckpoint() { phylo_tree->clearAllPartialLH(); } -void ModelMixture::getStateFrequency(double *state_freq, int mixture) { - ASSERT(mixture < getNMixtures()); - if (mixture >= 0) { - at(mixture)->getStateFrequency(state_freq); - return; - } - // special case: return weighted sum of state_freq across classes - // for mixture model, take the weighted sum of frequency vectors - double state_freq_class[num_states]; - int mix = getNMixtures(); - memset(state_freq, 0, sizeof(double)*num_states); - bool fused = isFused(); - for (int i = 0; i < mix; i++) { - at(i)->getStateFrequency(state_freq_class); - double weight = getMixtureWeight(i); - // fused model, take the weight from site_rate - if (fused) - weight = phylo_tree->getRate()->getProp(i) / (1.0 - phylo_tree->getRate()->getPInvar()); - for (int j = 0; j < num_states; j++) - state_freq[j] += weight*state_freq_class[j]; - } - // // DEBUG. - // cout << "Weighted state frequency of mixture component zero: "; - // for (int i = 0; i < num_states; i++) { - // cout << state_freq[i] << " "; - // } - // cout << endl; -} - // estimate the initial frequency vectors // method 1: given a set of classes in the mixture model, randomly assign each alignment position to one of the classes. // Then the nucleotide frequency array of each class is initialized according to the nucleotide frequencies among the positions assigned to the class. @@ -3947,29 +3918,78 @@ void ModelMixture::estimateInitFreq2() { } } -void ModelMixture::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { - ASSERT(mixture < getNMixtures()); - at(mixture)->computeTransMatrix(time, trans_matrix, 0, selected_row); +void ModelMixture::computeTransMatrix(double time, double *trans_matrix, int model_id, int selected_row) { + ASSERT(model_id > -1); // no default + at(model_id)->computeTransMatrix(time, trans_matrix, -1, selected_row); } -void ModelMixture::getQMatrix(double *q_mat, int mixture) -{ - ASSERT(mixture < getNMixtures()); - at(mixture)->getQMatrix(q_mat); +void ModelMixture::computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id) { + ASSERT(model_id > -1); // no default + at(model_id)->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2); } -void ModelMixture::computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture) { - ASSERT(mixture < getNMixtures()); - at(mixture)->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2); +double ModelMixture::computeTrans(double time, int state1, int state2, int model_id) { + ASSERT(model_id > -1); // no default + return at(model_id)->computeTrans(time, state1, state2); } -void ModelMixture::adaptStateFrequency(double* freq) -{ +double ModelMixture::computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id) { + ASSERT(model_id > -1); // no default + return at(model_id)->computeTrans(time, state1, state2, derv1, derv2); +} + +void ModelMixture::getRateMatrix(double *rate_mat, int model_id) { + ASSERT(model_id > -1); // no default + at(model_id)->getRateMatrix(rate_mat); +} + +void ModelMixture::setRateMatrix(double *rate_mat) { + ASSERT(false); // call for each submodel separately! +} + +void ModelMixture::getQMatrix(double *q_mat, int model_id) { + ASSERT(model_id > -1); // no default + at(model_id)->getQMatrix(q_mat); +} + +void ModelMixture::setQMatrix(double *q_mat, double *freq_vec) { + ASSERT(false); // call for each submodel separately! +} + +void ModelMixture::getStateFrequency(double *freq_vec, int model_id) { + ASSERT(model_id >= -1); + if (model_id >= 0) { + at(model_id)->getStateFrequency(freq_vec); + return; + } + // default: return the weighted sum of state_freq across classes + bool fused = isFused(); + double state_freq_class[num_states]; + std::fill_n(freq_vec, num_states, 0.0); + for (size_t m = 0; m < getNMixtures(); ++m) { + at(m)->getStateFrequency(state_freq_class); + double weight = getMixtureWeight(m); + if (fused) { + // fused model, take the weight from site_rate + weight = phylo_tree->getRate()->getProp(m) / (1.0 - phylo_tree->getRate()->getPInvar()); + } + for (int x = 0; x < num_states; ++x) { + freq_vec[x] += weight * state_freq_class[x]; + } + } +} + +void ModelMixture::setStateFrequency(double *freq_vec) { + ASSERT(false); // call for each submodel separately! +} + +void ModelMixture::adaptStateFrequency(double *freq_vec) { ASSERT(state_freq); for (iterator it = begin(); it != end(); it++) { if ((*it)->freq_type == FREQ_ESTIMATE || (*it)->freq_type == FREQ_EMPIRICAL) - (*it)->adaptStateFrequency(freq); + (*it)->adaptStateFrequency(freq_vec); } } diff --git a/model/modelmixture.h b/model/modelmixture.h index acbf81612..bc784c993 100644 --- a/model/modelmixture.h +++ b/model/modelmixture.h @@ -134,14 +134,6 @@ class ModelMixture: virtual public ModelMarkov, public vector { */ virtual void setMixtureClass(int cat, ModelSubst* m) { at(cat) = (ModelMarkov*)m; } - /** - compute the state frequency vector - @param mixture (optional) class for mixture model. - -1 to get weighted sum of class state frequency vector - @param state_freq (OUT) state frequency vector. Assume state_freq has size of num_states - */ - virtual void getStateFrequency(double *state_freq, int mixture = 0); - // estimate the initial frequence vector for the class // method 1: given a set of classes in the mixture model, randomly assign each alignment position to one of the classes. @@ -152,38 +144,29 @@ class ModelMixture: virtual public ModelMarkov, public vector { // The nucleotide frequency array of i-th class is initialized according to the nucleotide frequencies in the i-th partition void estimateInitFreq2(); - /** - compute the transition probability matrix. One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); - - /** - Get the rate matrix Q. One should override this function when defining new model. - The default is equal rate of 1 (JC Model), valid for all kind of data. - @param rate_mat (OUT) upper-triagle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 - */ - virtual void getQMatrix(double *q_mat, int mixture = 0); + virtual void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); + virtual void computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id = -1); - /** - compute the transition probability matrix.and the derivative 1 and 2 - @param time time between two events - @param mixture (optional) class for mixture model - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - @param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. - @param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. - */ - virtual void computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture = 0); + virtual double computeTrans(double time, int state1, int state2, int model_id = -1); + + virtual double computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id = -1); + + virtual void getRateMatrix(double *rate_mat, int model_id = -1); + + virtual void setRateMatrix(double *rate_mat); + + virtual void getQMatrix(double *q_mat, int model_id = -1); + + virtual void setQMatrix(double *q_mat, double *freq_vec); + + virtual void getStateFrequency(double *freq_vec, int model_id = -1); + + virtual void setStateFrequency(double *freq_vec); - virtual void adaptStateFrequency(double* freq); + virtual void adaptStateFrequency(double *freq_vec); /** @return the number of dimensions diff --git a/model/modelpomo.cpp b/model/modelpomo.cpp index c79a8d8e4..e7fd1360a 100644 --- a/model/modelpomo.cpp +++ b/model/modelpomo.cpp @@ -1090,9 +1090,7 @@ void ModelPoMo::set_heterozygosity_boundaries() { // generic like this one. However, I thought my problems may be related to the // implementation in ModelMarkov::computeTransMatrix(). Now I am not so sure // anymore. - -// TODO DS: The parameter mixture is unused at the moment. -void ModelPoMo::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { +void ModelPoMo::computeTransMatrix(double time, double *trans_matrix, int, int selected_row) { MatrixExpTechnique technique = phylo_tree->params->matrix_exp_technique; if (technique == MET_SCALING_SQUARING || !is_reversible) { // Do not change the object rate_matrix, but only trans_matrix. @@ -1142,7 +1140,7 @@ void ModelPoMo::computeTransMatrix(double time, double *trans_matrix, int mixtur } } - else ModelMarkov::computeTransMatrix(time, trans_matrix, 0, selected_row); + else ModelMarkov::computeTransMatrix(time, trans_matrix, -1, selected_row); } void ModelPoMo::computeTipLikelihood(PML::StateType state, double *lh) { diff --git a/model/modelpomo.h b/model/modelpomo.h index b26ef6e2d..8c170fcef 100644 --- a/model/modelpomo.h +++ b/model/modelpomo.h @@ -280,15 +280,7 @@ class ModelPoMo : virtual public ModelMarkov // empirical value. void set_heterozygosity_boundaries(); - /** - compute the transition probability matrix. - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); + virtual void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); /** * Set the scale factor of the mutation rates to NEW_SCALE. diff --git a/model/modelpomomixture.cpp b/model/modelpomomixture.cpp index 190fe9cd0..2f5cec251 100644 --- a/model/modelpomomixture.cpp +++ b/model/modelpomomixture.cpp @@ -265,7 +265,6 @@ bool ModelPoMoMixture::isUnstableParameters() { // version did not work for non-reversible substitution models. However, this // led to a clash because then computeTransMatrix is defined in both, // ModelMixture and ModelPoMo and inheritance is flawed. -void ModelPoMoMixture::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { - ASSERT(mixture < getNMixtures()); - at(mixture)->computeTransMatrix(time, trans_matrix, 0, selected_row); +void ModelPoMoMixture::computeTransMatrix(double time, double *trans_matrix, int model_id, int selected_row) { + ModelMixture::computeTransMatrix(time, trans_matrix, model_id, selected_row); } diff --git a/model/modelpomomixture.h b/model/modelpomomixture.h index c45fc6c90..1ebb341ab 100644 --- a/model/modelpomomixture.h +++ b/model/modelpomomixture.h @@ -156,16 +156,7 @@ class ModelPoMoMixture : public ModelPoMo, public ModelMixture { virtual void update_eigen_pointers(double *eval, double *evec , double *inv_evec, double *inv_evec_transposed); - - /** - compute the transition probability matrix. - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); + virtual void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); protected: diff --git a/model/modelset.cpp b/model/modelset.cpp index 038d1b08d..395435fbe 100644 --- a/model/modelset.cpp +++ b/model/modelset.cpp @@ -26,36 +26,28 @@ ModelSet::ModelSet(const char *model_name, PhyloTree *tree) : ModelMarkov(tree) full_name += "+site-specific state-frequency model (unpublished)"; } -void ModelSet::computeTransMatrix(double time, double* trans_matrix, int mixture, int selected_row) -{ - // TODO not working with vectorization - ASSERT(0); - for (iterator it = begin(); it != end(); it++) { - (*it)->computeTransMatrix(time, trans_matrix, mixture, selected_row); - trans_matrix += (num_states * num_states); - } +int ModelSet::getPtnModelID(size_t ptn) const { + ASSERT(ptn < size()); + return ptn; } -void ModelSet::computeTransDerv(double time, double* trans_matrix, double* trans_derv1, double* trans_derv2, int mixture) -{ +void ModelSet::computeTransMatrix(double time, double *trans_matrix, int model_id, int selected_row) { // TODO not working with vectorization ASSERT(0); - for (iterator it = begin(); it != end(); it++) { - (*it)->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2, mixture); - trans_matrix += (num_states * num_states); - trans_derv1 += (num_states * num_states); - trans_derv2 += (num_states * num_states); - } + ASSERT(model_id > -1); // no default + at(model_id)->computeTransMatrix(time, trans_matrix, -1, selected_row); } -int ModelSet::getPtnModelID(int ptn) -{ - ASSERT(ptn >= 0 && ptn < size()); - return ptn; +void ModelSet::computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id) { + // TODO not working with vectorization + ASSERT(0); + ASSERT(model_id > -1); // no default + at(model_id)->computeTransDerv(time, trans_matrix, trans_derv1, trans_derv2); } - -double ModelSet::computeTrans(double time, int model_id, int state1, int state2) { +double ModelSet::computeTrans(double time, int state1, int state2, int model_id) { + ASSERT(model_id > -1); // no default if (phylo_tree->vector_size == 1) { return at(model_id)->computeTrans(time, state1, state2); } @@ -77,7 +69,9 @@ double ModelSet::computeTrans(double time, int model_id, int state1, int state2) return trans_prob; } -double ModelSet::computeTrans(double time, int model_id, int state1, int state2, double &derv1, double &derv2) { +double ModelSet::computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id) { + ASSERT(model_id > -1); // no default if (phylo_tree->vector_size == 1) { return at(model_id)->computeTrans(time, state1, state2, derv1, derv2); } @@ -103,6 +97,42 @@ double ModelSet::computeTrans(double time, int model_id, int state1, int state2, return trans_prob; } +void ModelSet::getRateMatrix(double *rate_mat, int model_id) { + ASSERT(model_id > -1); // no default + at(model_id)->getRateMatrix(rate_mat); +} + +void ModelSet::setRateMatrix(double *rate_mat) { + ASSERT(false); // call for each submodel separately! +} + +void ModelSet::getQMatrix(double *q_mat, int model_id) { + ASSERT(model_id > -1); // no default + at(model_id)->getQMatrix(q_mat); +} + +void ModelSet::setQMatrix(double *q_mat, double *freq_vec) { + ASSERT(false); // call for each submodel separately! +} + +void ModelSet::getStateFrequency(double *freq_vec, int model_id) { + ASSERT(model_id >= -1); + if (model_id >= 0) { + at(model_id)->getStateFrequency(freq_vec); + return; + } + // default: return the +F freqs across all patterns + ModelMarkov::getStateFrequency(freq_vec); +} + +void ModelSet::setStateFrequency(double *freq_vec) { + ASSERT(false); // call for each submodel separately! +} + +void ModelSet::adaptStateFrequency(double *freq_vec) { + ASSERT(false); // call for each submodel separately! +} + int ModelSet::getNDim() { ASSERT(size()); @@ -125,16 +155,6 @@ void ModelSet::writeInfo(ostream& out) } } -void ModelSet::getStateFrequency(double *state_freq, int mixture) { - ASSERT(mixture >= -1); - if (mixture >= 0) { - at(mixture)->getStateFrequency(state_freq); - return; - } - // default: return the +F freqs across all patterns - ModelMarkov::getStateFrequency(state_freq); -} - void ModelSet::decomposeRateMatrix() { if (empty()) { diff --git a/model/modelset.h b/model/modelset.h index 413d11cf8..29ef8bfd2 100644 --- a/model/modelset.h +++ b/model/modelset.h @@ -31,93 +31,35 @@ class ModelSet : public ModelMarkov, public vector public: ModelSet(const char *model_name, PhyloTree *tree); - /** - * @return TRUE if this is a site-specific model, FALSE otherwise - */ - virtual bool isSiteSpecificModel() { return true; } - /** - * get the size of transition matrix, default is num_states*num_states. - * can be changed for e.g. site-specific model - */ - virtual int getTransMatrixSize() { return num_states * num_states * size(); } + virtual bool isSiteSpecificModel() { return true; } + virtual int getPtnModelID(size_t ptn) const; - /** - compute the transition probability matrix. - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); + virtual void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); - - /** - compute the transition probability matrix.and the derivative 1 and 2 - @param time time between two events - @param mixture (optional) class for mixture model - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - @param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. - @param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. - */ - virtual void computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture = 0); + virtual void computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id = -1); - /** - To AVOID 'hides overloaded virtual functions - compute the transition probability between two states - @param time time between two events - @param state1 first state - @param state2 second state - */ - virtual double computeTrans(double time, int state1, int state2) { return 0; } + virtual double computeTrans(double time, int state1, int state2, int model_id = -1); - /** - To AVOID 'hides overloaded virtual functions - compute the transition probability between two states - @param time time between two events - @param state1 first state - @param state2 second state - @param derv1 (OUT) 1st derivative - @param derv2 (OUT) 2nd derivative - */ - virtual double computeTrans(double time, int state1, int state2, double &derv1, double &derv2) { return 0; } + virtual double computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id = -1); + virtual void getRateMatrix(double *rate_mat, int model_id = -1); + virtual void setRateMatrix(double *rate_mat); - /** - compute the transition probability between two states at a specific site - One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param model_id model ID - @param state1 first state - @param state2 second state - */ - virtual double computeTrans(double time, int model_id, int state1, int state2); + virtual void getQMatrix(double *q_mat, int model_id = -1); - /** - compute the transition probability and its 1st and 2nd derivatives between two states at a specific site - One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param model_id model ID - @param state1 first state - @param state2 second state - @param derv1 (OUT) 1st derivative - @param derv2 (OUT) 2nd derivative - */ - virtual double computeTrans(double time, int model_id, int state1, int state2, double &derv1, double &derv2); + virtual void setQMatrix(double *q_mat, double *freq_vec); + + virtual void getStateFrequency(double *freq_vec, int model_id = -1); + + virtual void setStateFrequency(double *freq_vec); + + virtual void adaptStateFrequency(double *freq_vec); - /** - * @return pattern ID to model ID map, useful for e.g., partition model - * @param ptn pattern ID of the alignment - */ - virtual int getPtnModelID(int ptn); - /** return the number of dimensions */ @@ -130,13 +72,6 @@ class ModelSet : public ModelMarkov, public vector */ virtual void writeInfo(ostream &out); - /** - * Compute the state frequency vector - * @param mixture Submodel id of the state frequency vector, 0 for default frequencies - * @param[out] state_freq State frequency vector, must have the size of num_states - */ - virtual void getStateFrequency(double *state_freq, int mixture = -1); - /** decompose the rate matrix into eigenvalues and eigenvectors */ diff --git a/model/modelsubst.cpp b/model/modelsubst.cpp index ee209e011..3d9154ad2 100644 --- a/model/modelsubst.cpp +++ b/model/modelsubst.cpp @@ -1,7 +1,7 @@ // // C++ Implementation: substmodel // -// Description: +// Description: // // // Author: BUI Quang Minh, Steffen Klaere, Arndt von Haeseler , (C) 2008 @@ -10,19 +10,23 @@ // // #include "modelsubst.h" + #include "utils/tools.h" +#include -ModelSubst::ModelSubst(int nstates) : Optimization(), CheckpointFactory() -{ - num_states = nstates; - name = "JC"; - full_name = "JC (Juke and Cantor, 1969)"; - state_freq = new double[num_states]; - for (int i = 0; i < num_states; i++) - state_freq[i] = 1.0 / num_states; - freq_type = FREQ_EQUAL; +ModelSubst::ModelSubst(int nstates) { + num_states = nstates; fixed_parameters = false; -// linked_model = nullptr; + name = "JC"; + full_name = "JC (Jukes and Cantor, 1969)"; + freq_type = FREQ_EQUAL; + state_freq = new double[num_states]; + std::fill_n(state_freq, num_states, 1.0 / nstates); +} + +ModelSubst::~ModelSubst() { + delete [] state_freq; + state_freq = nullptr; } void ModelSubst::startCheckpoint() { @@ -31,11 +35,6 @@ void ModelSubst::startCheckpoint() { void ModelSubst::saveCheckpoint() { startCheckpoint(); -// CKP_SAVE(num_states); -// CKP_SAVE(name); -// CKP_SAVE(full_name); -// CKP_SAVE(freq_type); - // if (freq_type == FREQ_ESTIMATE && !fixed_parameters) // output the frequencies in any circumstances CKP_ARRAY_SAVE(num_states, state_freq); endCheckpoint(); @@ -45,189 +44,102 @@ void ModelSubst::saveCheckpoint() { void ModelSubst::restoreCheckpoint() { CheckpointFactory::restoreCheckpoint(); startCheckpoint(); -// CKP_RESTORE(num_states); -// CKP_RESTORE(name); -// CKP_RESTORE(full_name); -// int freq_type = this->freq_type; -// CKP_RESTORE(freq_type); -// this->freq_type = (StateFreqType)freq_type; - if (freq_type == FREQ_ESTIMATE && !fixed_parameters) + if (freq_type == FREQ_ESTIMATE && !fixed_parameters) { CKP_ARRAY_RESTORE(num_states, state_freq); + } endCheckpoint(); - decomposeRateMatrix(); } -// here the simplest Juke-Cantor model is implemented, valid for all kind of data (DNA, AA,...) -void ModelSubst::computeTransMatrix(double time, double *trans_matrix, int mixture, int selected_row) { - double non_diagonal = (1.0 - exp(-time*num_states/(num_states - 1))) / num_states; - double diagonal = 1.0 - non_diagonal * (num_states - 1); - int nstates_sqr = num_states * num_states; - - for (int i = 0; i < nstates_sqr; i++) - if (i % (num_states+1) == 0) - trans_matrix[i] = diagonal; - else - trans_matrix[i] = non_diagonal; -} - - -double ModelSubst::computeTrans(double time, int state1, int state2) { - double expt = exp(-time * num_states / (num_states-1)); - if (state1 != state2) { - return (1.0 - expt) / num_states; - } - return (1.0 + (num_states-1)*expt) / num_states; +// The following functions directly implement the simplest Jukes-Cantor model, +// which is valid for all kinds of data (DNA, AA, MORPH, etc) -/* double non_diagonal = (1.0 - exp(-time*num_states/(num_states - 1))) / num_states; - if (state1 != state2) - return non_diagonal; - return 1.0 - non_diagonal * (num_states - 1);*/ -} - -double ModelSubst::computeTrans(double time, int model_id, int state1, int state2) { - return computeTrans(time, state1, state2); -} - -double ModelSubst::computeTrans(double time, int state1, int state2, double &derv1, double &derv2) { - double coef = -double(num_states) / (num_states-1); - double expt = exp(time * coef); - if (state1 != state2) { - derv1 = expt / (num_states-1); - derv2 = derv1 * coef; - return (1.0 - expt) / num_states; - } - - derv1 = -expt; - derv2 = derv1 * coef; - return (1.0 + (num_states-1)*expt) / num_states; -} - -double ModelSubst::computeTrans(double time, int model_id, int state1, int state2, double &derv1, double &derv2) { - return computeTrans(time, state1, state2, derv1, derv2); +void ModelSubst::computeTransMatrix(double time, double *trans_matrix, int, int) { + double coef = -double(num_states) / (num_states-1); + double expt = exp(time * coef); + double lh_non_diag = (1.0 - expt) / num_states; + double lh_diag = 1.0 - (lh_non_diag * (num_states-1)); + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = 0; j < num_states; ++j, ++k) { + trans_matrix[k] = (i == j) ? lh_diag : lh_non_diag; + } + } } -void ModelSubst::getRateMatrix(double *rate_mat) { - int nrate = getNumRateEntries(); - for (int i = 0; i < nrate; i++) - rate_mat[i] = 1.0; +void ModelSubst::computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int) { + double coef = -double(num_states) / (num_states-1); + double expt = exp(time * coef); + double lh_non_diag = (1.0 - expt) / num_states; + double lh_diag = 1.0 - (lh_non_diag * (num_states-1)); + double derv1_non_diag = expt / (num_states-1); + double derv1_diag = -expt; + double derv2_non_diag = derv1_non_diag * coef; + double derv2_diag = derv1_diag * coef; + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = 0; j < num_states; ++j, ++k) { + trans_matrix[k] = (i == j) ? lh_diag : lh_non_diag; + trans_derv1[k] = (i == j) ? derv1_diag : derv1_non_diag; + trans_derv2[k] = (i == j) ? derv2_diag : derv2_non_diag; + } + } } -void ModelSubst::getQMatrix(double *q_mat, int mixture) { - int i, j, k; - for (i = 0, k = 0; i < num_states; i++) - for (j = 0; j < num_states; j++, k++) - if (i == j) q_mat[k] = -1.0; else q_mat[k] = 1.0/3; +double ModelSubst::computeTrans(double time, int state1, int state2, int) { + double coef = -double(num_states) / (num_states-1); + double expt = exp(time * coef); + if (state1 != state2) { + return (1.0 - expt) / num_states; + } + return (1.0 + (expt * (num_states-1))) / num_states; } -void ModelSubst::getStateFrequency(double *state_freq, int mixture) { - double freq = 1.0 / num_states; - for (int i = 0; i < num_states; i++) - state_freq[i] = freq; +double ModelSubst::computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int) { + double coef = -double(num_states) / (num_states-1); + double expt = exp(time * coef); + if (state1 != state2) { + derv1 = expt / (num_states-1); + derv2 = derv1 * coef; + return (1.0 - expt) / num_states; + } + derv1 = -expt; + derv2 = derv1 * coef; + return (1.0 + (expt * (num_states-1))) / num_states; } -void ModelSubst::setStateFrequency(double *state_freq) { - memcpy(this->state_freq, state_freq, sizeof(double)*num_states); +void ModelSubst::getRateMatrix(double *rate_mat, int) { + int nrate = getNumRateEntries(); + std::fill_n(rate_mat, nrate, 1.0); } -void ModelSubst::computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture) -{ - double expf = exp(-time*num_states/(num_states - 1)); - double non_diag = (1.0 - expf) / num_states; - double diag = 1.0 - non_diag * (num_states - 1); - double derv1_non_diag = expf / (num_states-1); - double derv1_diag = -expf; - double derv2_non_diag = -derv1_non_diag*num_states/(num_states-1); - double derv2_diag = -derv1_diag*num_states/(num_states-1); - - int nstates_sqr = num_states * num_states; - int i; - for (i = 0; i < nstates_sqr; i++) - if (i % (num_states+1) == 0) { - trans_matrix[i] = diag; - trans_derv1[i] = derv1_diag; - trans_derv2[i] = derv2_diag; - } else { - trans_matrix[i] = non_diag; - trans_derv1[i] = derv1_non_diag; - trans_derv2[i] = derv2_non_diag; - } - - // DEBUG - /*int j; - if (verbose_mode == VB_DEBUG) { - cout.precision(4); - cout << "time = " << time << endl; - for (i = 0; i < num_states; i++, cout << endl) { - for (j = 0; j < num_states; j++) { - cout.width(8); - cout << right << trans_matrix[i*num_states+j] << " "; - } - cout << "| "; - for (j = 0; j < num_states; j++) { - cout << right << trans_derv1[i*num_states+j] << " "; - cout.width(8); - } - cout << "| "; - for (j = 0; j < num_states; j++) { - cout.width(8); - cout << right << trans_derv2[i*num_states+j] << " "; - } - } - cout.precision(10); - }*/ - +void ModelSubst::getQMatrix(double *q_mat, int) { + double q_non_diag = 1.0 / (num_states-1); + for (int i = 0, k = 0; i < num_states; ++i) { + for (int j = 0; j < num_states; ++j, ++k) { + q_mat[k] = (i == j) ? -1.0 : q_non_diag; + } + } } -void ModelSubst::multiplyWithInvEigenvector(double *state_lk) { - int mnstates = get_safe_upper_limit(num_states); - int nmixtures = getNMixtures(); - double *inv_eigenvectors = getInverseEigenvectors(); - double saved_state_lk[num_states]; - memcpy(saved_state_lk, state_lk, sizeof(double)*num_states); - memset(state_lk, 0, sizeof(double)*num_states*nmixtures); - for (int m = 0; m < nmixtures; m++) { - double *inv_evec = &inv_eigenvectors[m * mnstates * num_states]; - double *this_state_lk = &state_lk[m*num_states]; - for (int i = 0; i < num_states; i++) - for (int j = 0; j < num_states; j++, inv_evec++) - this_state_lk[i] += (*inv_evec) * saved_state_lk[j]; - } +void ModelSubst::getStateFrequency(double *freq_vec, int) { + std::fill_n(freq_vec, num_states, 1.0 / num_states); } void ModelSubst::computeTipLikelihood(PML::StateType state, double *state_lk) { if (state < num_states) { // single state - memset(state_lk, 0, num_states*sizeof(double)); + std::fill_n(state_lk, num_states, 0.0); state_lk[state] = 1.0; } else { // unknown state - for (int i = 0; i < num_states; i++) - state_lk[i] = 1.0; + std::fill_n(state_lk, num_states, 1.0); } } -double *ModelSubst::newTransMatrix() { - return new double[num_states * num_states]; -} - -void ModelSubst::printMrBayesModelText(ofstream& out, string partition, string charset) { +void ModelSubst::printMrBayesModelText(ofstream &out, string partition, string charset) { out << "using MrBayes model GTR+G+I]" << endl; out << " [Model not supported by MrBayes, defaulting to GTR+G+I (DNA)]" << endl; outWarning("MrBayes output is not supported by model " + name + ", defaulting to GTR+G+I (DNA)!"); - out << " lset applyto=(" << partition << ") nucmodel=4by4 nst=" << 6 << " rates=" << "invgamma" << ";" << endl; } - -ModelSubst::~ModelSubst() -{ - // mem space pointing to target model and thus avoid double free here -// if (linked_model && linked_model != this) -// return; - - if (state_freq) delete [] state_freq; -} - - - diff --git a/model/modelsubst.h b/model/modelsubst.h index 51a761f96..02a7bb884 100644 --- a/model/modelsubst.h +++ b/model/modelsubst.h @@ -178,106 +178,94 @@ class ModelSubst: public Optimization, public CheckpointFactory return 0; } - /** - * get the size of transition matrix, default is num_states*num_states. - * can be changed for e.g. site-specific model - */ - virtual int getTransMatrixSize() { return num_states * num_states; } - - /** - compute the transition probability matrix. One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param mixture (optional) class for mixture model - @param selected_row (optional) only compute the entries of one selected row. By default, compute all rows - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - */ - virtual void computeTransMatrix(double time, double *trans_matrix, int mixture = 0, int selected_row = -1); - - /** - compute the transition probability between two states. - One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param state1 first state - @param state2 second state - */ - virtual double computeTrans(double time, int state1, int state2); + /** + * Map alignment patterns to submodels (for site-specific models) + * @param ptn ID of an alignment pattern + * @return ID of the corresponding model + */ + virtual int getPtnModelID(size_t ptn) const { return -1; } - /** - compute the transition probability between two states at a specific model ID, useful for partition model - One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param model_id model ID - @param state1 first state - @param state2 second state - */ - virtual double computeTrans(double time, int model_id, int state1, int state2); + /** + * Compute the transition probability matrix on a branch. + * The default is the JC model, valid for all kinds of data + * @param time The branch length + * @param model_id ID of a submodel (for mixture and site-specific models) + * @param selected_row Only compute the entries for the selected row, + * the default is to compute entries for all rows + * @param[out] trans_matrix The transition matrix between all pairs of states, + * assumed to have the size of num_states*num_states + */ + virtual void computeTransMatrix(double time, double *trans_matrix, int model_id = -1, int selected_row = -1); - /** - compute the transition probability and its 1st and 2nd derivatives between two states. - One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param state1 first state - @param state2 second state - @param derv1 (OUT) 1st derivative - @param derv2 (OUT) 2nd derivative - */ - virtual double computeTrans(double time, int state1, int state2, double &derv1, double &derv2); + /** + * The same as computeTransMatrix() above, but also computes the + * 1st and 2nd derivative matrices with respect to the branch length + * @param[out] trans_matrix The transition matrix between all pairs of states, + * assumed to have the size of num_states*num_states + * @param[out] trans_derv1 The 1st derivative matrix between all pairs of states + * @param[out] trans_derv2 The 2nd derivative matrix between all pairs of states + */ + virtual void computeTransDerv(double time, double *trans_matrix, + double *trans_derv1, double *trans_derv2, int model_id = -1); - /** - compute the transition probability and its 1st and 2nd derivatives between two states at a specific model ID - One should override this function when defining new model. - The default is the Juke-Cantor model, valid for all kind of data (DNA, AA, Codon, etc) - @param time time between two events - @param model_id model ID - @param state1 first state - @param state2 second state - @param derv1 (OUT) 1st derivative - @param derv2 (OUT) 2nd derivative - */ - virtual double computeTrans(double time, int model_id, int state1, int state2, double &derv1, double &derv2); + /** + * Compute the transition probability between the two states on a branch. + * The default is the JC model, valid for all kinds of data + * @param time The branch length between the two states + * @param model_id ID of a submodel (for mixture and site-specific models) + * @param state1 The start state + * @param state2 The end state + * @param[out] derv1 The 1st derivative + * @param[out] derv2 The 2nd derivative + * @return The transition probability + */ + virtual double computeTrans(double time, int state1, int state2, int model_id = -1); - /** - * @return pattern ID to model ID map, useful for e.g., partition model - * @param ptn pattern ID of the alignment - */ - virtual int getPtnModelID(int ptn) { return 0; } - + /** + * The same as computeTrans() above, but also computes the + * 1st and 2nd derivatives with respect to the branch length + * @param[out] derv1 The 1st derivative + * @param[out] derv2 The 2nd derivative + * @return The transition probability + */ + virtual double computeTrans(double time, int state1, int state2, + double &derv1, double &derv2, int model_id = -1); - /** - * Get the rate parameters like a,b,c,d,e,f for DNA model!!! - Get the above-diagonal entries of the rate matrix, assuming that the last element is 1. - ONE SHOULD OVERRIDE THIS FUNCTION WHEN DEFINING NEW MODEL!!! - The default is equal rate of 1 (JC Model), valid for all kind of data. - @param rate_mat (OUT) upper-triangle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 - */ - - virtual void getRateMatrix(double *rate_mat); + /** + * Get the rate parameters, such as a,b,c,d,e,f for a DNA model. + * Get the above-diagonal entries of the rate matrix, assuming that + * the last element is 1. + * The default is equal rates of 1 (JC Model), valid for all kinds of data + * @param[out] rate_mat An upper-triangle rate matrix, assumed to have the + * size of num_states*(num_states-1)/2 + * @param model_id ID of a submodel (for mixture and site-specific models) + */ + virtual void getRateMatrix(double *rate_mat, int model_id = -1); - /** - Get the rate matrix Q. One should override this function when defining new model. - The default is equal rate of 1 (JC Model), valid for all kind of data. - @param rate_mat (OUT) upper-triagle rate matrix. Assume rate_mat has size of num_states*(num_states-1)/2 - */ - virtual void getQMatrix(double *q_mat, int mixture = 0); + /** + * Get the instantaneous rate matrix Q. + * The default is derived from equal rates and equal state frequencies + * @param[out] q_mat A full matrix: qij >= 0, qii = -sum_j qij (j != i), + * assumed to have the size of num_states*num_states + * @param model_id ID of a submodel (for mixture and site-specific models) + */ + virtual void getQMatrix(double *q_mat, int model_id = -1); - /** - compute the state frequency vector. One should override this function when defining new model. - The default is equal state sequency, valid for all kind of data. - @param mixture (optional) class for mixture model - @param[out] state_freq state frequency vector. Assume state_freq has size of num_states - */ - virtual void getStateFrequency(double *state_freq, int mixture = 0); + /** + * Get the state frequency vector. + * The default is equal state frequencies, valid for all kinds of data + * @param[out] freq_vec A state frequency vector, assumed to have the + * size of num_states + * @param model_id ID of a submodel (for mixture and site-specific models) + */ + virtual void getStateFrequency(double *freq_vec, int model_id = -1); /** - set the state frequency vector. - @param state_freq state frequency vector. Assume state_freq has size of num_states + * Set the state frequency vector + * @param freq_vec A state frequency vector, assumed to have the + * size of num_states */ - virtual void setStateFrequency(double *state_freq); + virtual void setStateFrequency(double *freq_vec) {} /** get frequency type @@ -290,38 +278,19 @@ class ModelSubst: public Optimization, public CheckpointFactory @param tree the associated tree */ virtual void setTree(PhyloTree *tree) {} - - /** for reversible models, multiply likelihood with inverse eigenvectors for fast pruning algorithm - @param[in/out] state_lk state likelihood multiplied with inverse eigenvectors + + /** + * For reversible models, multiply the partial likelihood vector with + * the matrix of inverse eigenvectors for the fast pruning algorithm + * @param[in/out] state_lh The partial likelihood vector */ - void multiplyWithInvEigenvector(double *state_lk); + virtual void multiplyWithInvEigenvector(double *state_lh) {} /** compute the tip likelihood vector of a state for Felsenstein's pruning algorithm @param state character state @param[out] state_lk state likehood vector of size num_states */ virtual void computeTipLikelihood(PML::StateType state, double *state_lk); - - /** - allocate memory for a transition matrix. One should override this function when defining new model - such as Gamma model. The default is to allocate a double vector of size num_states * num_states. This - is equivalent to the memory needed by a square matrix. - @return the pointer to the newly allocated transition matrix - */ - virtual double *newTransMatrix(); - - - /** - compute the transition probability matrix.and the derivative 1 and 2 - @param time time between two events - @param mixture (optional) class for mixture model - @param trans_matrix (OUT) the transition matrix between all pairs of states. - Assume trans_matrix has size of num_states * num_states. - @param trans_derv1 (OUT) the 1st derivative matrix between all pairs of states. - @param trans_derv2 (OUT) the 2nd derivative matrix between all pairs of states. - */ - virtual void computeTransDerv(double time, double *trans_matrix, - double *trans_derv1, double *trans_derv2, int mixture = 0); /** decompose the rate matrix into eigenvalues and eigenvectors diff --git a/model/modelunrest.cpp b/model/modelunrest.cpp index 066e69f7a..5de8ce126 100644 --- a/model/modelunrest.cpp +++ b/model/modelunrest.cpp @@ -110,7 +110,7 @@ void ModelUnrest::setRates() { return; } -void ModelUnrest::setStateFrequency(double* freq) { +void ModelUnrest::setStateFrequency(double *freq_vec) { // DOES NOTHING } diff --git a/model/modelunrest.h b/model/modelunrest.h index ebe8cfcdb..56965af4c 100644 --- a/model/modelunrest.h +++ b/model/modelunrest.h @@ -30,7 +30,7 @@ class ModelUnrest: public ModelMarkov { set the state frequency vector. @param state_freq state frequency vector. Assume state_freq has size of num_states */ - virtual void setStateFrequency(double *state_freq); + virtual void setStateFrequency(double *freq_vec); /** start structure for checkpointing diff --git a/model/partitionmodel.cpp b/model/partitionmodel.cpp index 6d18757de..185e3a532 100644 --- a/model/partitionmodel.cpp +++ b/model/partitionmodel.cpp @@ -107,7 +107,7 @@ PartitionModel::PartitionModel(Params ¶ms, PhyloSuperTree *tree, ModelsBlock cout << "rowsum: " << Q.rowwise().sum() << endl; Map >(pair_freq, nstates, nstates) = Q; */ - ((ModelMarkov*)linked_models.begin()->second)->setFullRateMatrix(pair_freq, state_freq); + ((ModelMarkov*)linked_models.begin()->second)->setQMatrix(pair_freq, state_freq); ((ModelMarkov*)linked_models.begin()->second)->decomposeRateMatrix(); delete [] state_freq; delete [] pair_freq;