From dccf6c4a61b498d98551d40635d92e2fe4dac261 Mon Sep 17 00:00:00 2001 From: StefanFlaumberg Date: Wed, 9 Sep 2026 04:42:58 +0300 Subject: [PATCH] Use ptn_rate_mat similarly to ptn_state_freq, ensure they are copied by alignment-copying functions --- alignment/alignment.cpp | 109 ++++++++++++++++++++++++++++----------- alignment/alignment.h | 15 +++--- model/modelfactory.cpp | 21 +++----- utils/mutsel_wrapper.cpp | 58 ++++++++++++--------- 4 files changed, 129 insertions(+), 74 deletions(-) diff --git a/alignment/alignment.cpp b/alignment/alignment.cpp index 9a4a9526b..78bbcffdb 100644 --- a/alignment/alignment.cpp +++ b/alignment/alignment.cpp @@ -75,6 +75,7 @@ Alignment::Alignment() { non_stop_codon = nullptr; seq_type = SEQ_UNKNOWN; STATE_UNKNOWN = 126; + num_rates = 0; // pars_lower_bound = nullptr; // now a local variable in orderPatternByNumChars() } @@ -85,6 +86,11 @@ Alignment::~Alignment() { non_stop_codon = nullptr; // delete [] pars_lower_bound; // now a local variable in orderPatternByNumChars() // pars_lower_bound = nullptr; + for (vector::reverse_iterator rit = ptn_rate_mat.rbegin(); rit != ptn_rate_mat.rend(); ++rit) { + delete [] (*rit); + (*rit) = nullptr; + } + ptn_rate_mat.clear(); for (vector::reverse_iterator rit = ptn_state_freq.rbegin(); rit != ptn_state_freq.rend(); ++rit) { delete [] (*rit); (*rit) = nullptr; @@ -1589,6 +1595,13 @@ void Alignment::regroupSitePattern(const IntVector &site_group) { } } // refill the existing pattern-specific parameters + if (isSSM()) { + vector stored_ptn_rate_mat = ptn_rate_mat; + ptn_rate_mat.clear(); + for (size_t ptn = 0; ptn < getNPattern(); ++ptn) { + ptn_rate_mat.push_back(stored_ptn_rate_mat[new_to_old_pattern[ptn]]); + } + } if (isSSF()) { vector stored_ptn_state_freq = ptn_state_freq; ptn_state_freq.clear(); @@ -3623,6 +3636,7 @@ Alignment *Alignment::initAlignmentCopy() const { // Alignment members aln->seq_names = seq_names; aln->seq_type = seq_type; + aln->num_rates = num_rates; aln->num_states = num_states; aln->STATE_UNKNOWN = STATE_UNKNOWN; if (aln->seq_type == SEQ_CODON) { @@ -3665,14 +3679,23 @@ Alignment* Alignment::extractSubAlignment(const IntVector &seq_id, size_t true_chars = total_chars - gap_chars; if (true_chars >= min_true_chars) { bool added = aln->addPattern(pat); + if (isSSM() && added) { + // a new pattern is added, copy its rate matrix + double *rate_mat = nullptr; + if (ptn_rate_mat[ptn]) { + rate_mat = new double[num_rates]; + memcpy(rate_mat, ptn_rate_mat[ptn], num_rates*sizeof(double)); + } + aln->ptn_rate_mat.push_back(rate_mat); + } if (isSSF() && added) { // a new pattern is added, copy its state frequency vector - double *state_freqs = nullptr; + double *state_freq = nullptr; if (ptn_state_freq[ptn]) { - state_freqs = new double[num_states]; - memcpy(state_freqs, ptn_state_freq[ptn], num_states*sizeof(double)); + state_freq = new double[num_states]; + memcpy(state_freq, ptn_state_freq[ptn], num_states*sizeof(double)); } - aln->ptn_state_freq.push_back(state_freqs); + aln->ptn_state_freq.push_back(state_freq); } } // site is examined, add to progress @@ -3702,14 +3725,23 @@ Alignment *Alignment::extractPatterns(const IntVector &ptn_id) const { int ptn = *it; Pattern pat = at(ptn); bool added = aln->addPattern(pat); + if (isSSM() && added) { + // a new pattern is added, copy its rate matrix + double *rate_mat = nullptr; + if (ptn_rate_mat[ptn]) { + rate_mat = new double[num_rates]; + memcpy(rate_mat, ptn_rate_mat[ptn], num_rates*sizeof(double)); + } + aln->ptn_rate_mat.push_back(rate_mat); + } if (isSSF() && added) { // a new pattern is added, copy its state frequency vector - double *state_freqs = nullptr; + double *state_freq = nullptr; if (ptn_state_freq[ptn]) { - state_freqs = new double[num_states]; - memcpy(state_freqs, ptn_state_freq[ptn], num_states*sizeof(double)); + state_freq = new double[num_states]; + memcpy(state_freq, ptn_state_freq[ptn], num_states*sizeof(double)); } - aln->ptn_state_freq.push_back(state_freqs); + aln->ptn_state_freq.push_back(state_freq); } } aln->countConstSites(); @@ -3727,14 +3759,23 @@ Alignment *Alignment::extractPatternFreqs(const IntVector &ptn_freq) const { Pattern pat = at(ptn); pat.frequency = ptnf; bool added = aln->addPattern(pat); + if (isSSM() && added) { + // a new pattern is added, copy its rate matrix + double *rate_mat = nullptr; + if (ptn_rate_mat[ptn]) { + rate_mat = new double[num_rates]; + memcpy(rate_mat, ptn_rate_mat[ptn], num_rates*sizeof(double)); + } + aln->ptn_rate_mat.push_back(rate_mat); + } if (isSSF() && added) { // a new pattern is added, copy its state frequency vector - double *state_freqs = nullptr; + double *state_freq = nullptr; if (ptn_state_freq[ptn]) { - state_freqs = new double[num_states]; - memcpy(state_freqs, ptn_state_freq[ptn], num_states*sizeof(double)); + state_freq = new double[num_states]; + memcpy(state_freq, ptn_state_freq[ptn], num_states*sizeof(double)); } - aln->ptn_state_freq.push_back(state_freqs); + aln->ptn_state_freq.push_back(state_freq); } } } @@ -3752,14 +3793,23 @@ Alignment *Alignment::extractSites(const IntVector &site_id) const { Pattern pat = at(ptn); pat.frequency = 1; bool added = aln->addPattern(pat); + if (isSSM() && added) { + // a new pattern is added, copy its rate matrix + double *rate_mat = nullptr; + if (ptn_rate_mat[ptn]) { + rate_mat = new double[num_rates]; + memcpy(rate_mat, ptn_rate_mat[ptn], num_rates*sizeof(double)); + } + aln->ptn_rate_mat.push_back(rate_mat); + } if (isSSF() && added) { // a new pattern is added, copy its state frequency vector - double *state_freqs = nullptr; + double *state_freq = nullptr; if (ptn_state_freq[ptn]) { - state_freqs = new double[num_states]; - memcpy(state_freqs, ptn_state_freq[ptn], num_states*sizeof(double)); + state_freq = new double[num_states]; + memcpy(state_freq, ptn_state_freq[ptn], num_states*sizeof(double)); } - aln->ptn_state_freq.push_back(state_freqs); + aln->ptn_state_freq.push_back(state_freq); } } aln->countConstSites(); @@ -4051,6 +4101,7 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq position_spec = aln->position_spec; aln_file = aln->aln_file; seq_names.insert(seq_names.begin(), aln->seq_names.begin(), aln->seq_names.end()); + num_rates = aln->num_rates; num_states = aln->num_states; seq_type = aln->seq_type; genetic_code = aln->genetic_code; @@ -4075,7 +4126,7 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq pattern_freq->resize(0); pattern_freq->resize(aln->getNPattern(), 0); } - if (aln->isSSF() && spec) { + if ((aln->isSSM() || aln->isSSF()) && spec) { // resampling also the per-site state frequency vector outError("Unsupported bootstrap feature, pls contact the developers"); } @@ -4094,23 +4145,23 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq Pattern pat = aln->at(ptn); pat.frequency = 1; bool added = addPattern(pat); + if (aln->isSSM() && added) { + // a new pattern is added, copy its rate matrix + double *rate_mat = nullptr; + if (aln->ptn_rate_mat[ptn]) { + rate_mat = new double[num_rates]; + memcpy(rate_mat, aln->ptn_rate_mat[ptn], num_rates*sizeof(double)); + } + ptn_rate_mat.push_back(rate_mat); + } if (aln->isSSF() && added) { - // a new pattern is added, copy state frequency vector + // a new pattern is added, copy its state frequency vector double *state_freq = nullptr; if (aln->ptn_state_freq[ptn]) { state_freq = new double[num_states]; memcpy(state_freq, aln->ptn_state_freq[ptn], num_states*sizeof(double)); } ptn_state_freq.push_back(state_freq); - if (!aln->site_rate_matrices.empty()) { - /* Minh/Thomas: Better change 190 to num_states*(numstates-1)/2 so that if you want to - extend the model in the future, no change is needed here - Also: This is only for reversible models. For non-rev models - you need to store the full matrix, i.e., num_states*num_states entries - */ - const double *rate_matrix = aln->site_rate_matrices.data() + ptn * 190; - site_rate_matrices.insert(site_rate_matrices.end(), rate_matrix, rate_matrix + 190); - } } if (pattern_freq) { ((*pattern_freq)[ptn])++; @@ -4195,8 +4246,8 @@ void Alignment::createBootstrapAlignment(Alignment *aln, IntVector* pattern_freq out_site += site_vec[part+1]; } } - if (!aln->site_rate_matrices.empty()) { - ASSERT(aln->site_rate_matrices.size() == aln->getNPattern() * 190); + if (aln->isSSM()) { + ASSERT(ptn_rate_mat.size() == getNPattern()); } if (aln->isSSF()) { ASSERT(ptn_state_freq.size() == getNPattern()); diff --git a/alignment/alignment.h b/alignment/alignment.h index f42233a1c..0866ea39f 100644 --- a/alignment/alignment.h +++ b/alignment/alignment.h @@ -552,6 +552,7 @@ class Alignment : public vector, public CharSet, public StateSpace { */ bool isGapOnlySeq(int seq) const; + bool isSSM() const { return !ptn_rate_mat.empty(); } bool isSSF() const { return !ptn_state_freq.empty(); } virtual bool isSuperAlignment() const { return false; } @@ -944,15 +945,17 @@ class Alignment : public vector, public CharSet, public StateSpace { vector pomo_sampled_states; IntIntMap pomo_sampled_states_index; // indexing, to quickly find if a PoMo-2-state is already present - /* for site-specific state frequency model with Huaichun, Edward, Andrew */ + /* for site-specific models */ - /** pattern index to state frequency vector map */ + /** the size of a rate matrix in ptn_rate_mat */ + int num_rates; + + /** pattern ID to rate matrix map */ + vector ptn_rate_mat; + + /** pattern ID to state frequency vector map */ vector ptn_state_freq; - /** site to rate matrix. Stored in row-major order [num_sites, 190] used for MUTSEL */ - vector site_rate_matrices; - // Minh/Thomas: TODO rename this to ptn_rate_matrices to avoid confusions - /** * @return true if data type is SEQ_CODON and state is a stop codon */ diff --git a/model/modelfactory.cpp b/model/modelfactory.cpp index d677312e6..2555f208b 100644 --- a/model/modelfactory.cpp +++ b/model/modelfactory.cpp @@ -639,8 +639,7 @@ ModelFactory::ModelFactory(Params ¶ms, string &model_name, PhyloTree *tree, // fused_mix_rate &= model->isMixture() && site_rate->getNRate() > 1; } else { // site-specific model - - if (tree->aln->site_rate_matrices.empty()) { + if (!tree->aln->isSSM()) { // PMSF if (model_str == "JC" || model_str == "POISSON") outError("JC is not suitable for site-specific model"); @@ -660,29 +659,26 @@ ModelFactory::ModelFactory(Params ¶ms, string &model_name, PhyloTree *tree, modeli->setStateFrequency(state_freq); modeli->setRateMatrix(rates); } - if (tree->aln->ptn_state_freq[i]) - modeli->setStateFrequency (tree->aln->ptn_state_freq[i]); - + if (tree->aln->ptn_state_freq[i]) { + modeli->setStateFrequency(tree->aln->ptn_state_freq[i]); + } modeli->init(FREQ_USER_DEFINED); models->push_back(modeli); } delete [] rates; delete [] state_freq; - models->joinEigenMemory(); models->decomposeRateMatrix(); } else { // MUTSEL - ModelSet *models = new ModelSet(model_str.c_str(), tree); - + model = new ModelSet(model_str.c_str(), tree); + ModelSet *models = (ModelSet*)model; // assign pointer for convenience models->init((params.freq_type != FREQ_UNKNOWN) ? params.freq_type : FREQ_EMPIRICAL); models->fixParameters(true); - for (size_t i = 0; i < tree->aln->ptn_state_freq.size(); ++i) { ModelMarkov *modeli = new ModelMarkov(tree, true, true); // dummy model for getting num_states and num_rate_entries modeli->setStateFrequency(tree->aln->ptn_state_freq[i]); - modeli->setRateMatrix(tree->aln->site_rate_matrices.data() + i * 190); - + modeli->setRateMatrix(tree->aln->ptn_rate_mat[i]); // Minh/Thomas: Important to note: site Q matrices are not normalised // But make sure to normalise across all sites, i.e. // (1/nsites) * sum_i mu_i = 1.0, where mu_i = -sum_j pi^i_j * Q^i_{jj} @@ -691,11 +687,8 @@ ModelFactory::ModelFactory(Params ¶ms, string &model_name, PhyloTree *tree, modeli->fixParameters(true); models->push_back(modeli); } - models->joinEigenMemory(); models->decomposeRateMatrix(); - - model = models; } } diff --git a/utils/mutsel_wrapper.cpp b/utils/mutsel_wrapper.cpp index 29aa07bf9..05acc2a36 100644 --- a/utils/mutsel_wrapper.cpp +++ b/utils/mutsel_wrapper.cpp @@ -20,6 +20,7 @@ void rust_mutsel(int32_t *parents, std::cout << "Mutsel support not compiled in!" << std::endl; exit(1); } + void rust_set_rayon_threads(int32_t num_threads) { return; @@ -219,8 +220,13 @@ std::string read_binary_site_model_file_internal(std::string &filename, std::vec void write_site_models_to_alignment(Alignment &alignment, const double *site_freq, const double *rate_matrices, int len) { - alignment.ptn_state_freq.clear(); - alignment.site_rate_matrices.clear(); + ASSERT(alignment.ptn_rate_mat.empty() && + alignment.ptn_state_freq.empty()); + + // currently we only support 20 states for mutsel model, + // so this function should only be called for protein alignments + ASSERT(alignment.num_states == 20); + alignment.num_rates = 190; size_t nsite = alignment.getNSite(); if (len != static_cast(nsite)) @@ -239,16 +245,15 @@ void write_site_models_to_alignment(Alignment &alignment, const double *site_fre } bool aln_changed = false; - ASSERT(alignment.num_states == 20); // currently we only support 20 states for mutsel model, so this function should only be called for protein alignments vector models_freq; - vector models_rate; + vector models_rate; for (size_t site = 0; site < nsite; ++site) { site_model[site] = models_freq.size(); - const double *freq = site_freq + site * 20; - const double *rate_para = rate_matrices + site * 190; + const double *state_freq_ptr = site_freq + site * 20; + const double *rate_mat_ptr = rate_matrices + site * 190; bool add = true; int first_site = pattern_first_site[alignment.getPatternID(site)]; @@ -258,7 +263,7 @@ void write_site_models_to_alignment(Alignment &alignment, const double *site_fre bool matched_freq_and_rate = true; for (int i = 0; i < 20; ++i) { - if (freq[i] != models_freq[first_model][i]) + if (state_freq_ptr[i] != models_freq[first_model][i]) { matched_freq_and_rate = false; break; @@ -268,7 +273,7 @@ void write_site_models_to_alignment(Alignment &alignment, const double *site_fre { for (int i = 0; i < 190; ++i) { - if (rate_para[i] != models_rate[first_model][i]) + if (rate_mat_ptr[i] != models_rate[first_model][i]) { matched_freq_and_rate = false; break; @@ -290,9 +295,11 @@ void write_site_models_to_alignment(Alignment &alignment, const double *site_fre if (add) { double *site_freq_entry = new double[20]; - memcpy(site_freq_entry, freq, sizeof(double) * 20); + memcpy(site_freq_entry, state_freq_ptr, sizeof(double) * 20); models_freq.push_back(site_freq_entry); - models_rate.emplace_back(rate_para, rate_para + 190); + double *site_rate_entry = new double[190]; + memcpy(site_rate_entry, rate_mat_ptr, sizeof(double) * 190); + models_rate.push_back(site_rate_entry); } } @@ -310,28 +317,28 @@ void write_site_models_to_alignment(Alignment &alignment, const double *site_fre } } - vector used_model(models_freq.size(), false); + size_t used_models = 0; + vector model_used(models_freq.size(), false); for (size_t ptn = 0; ptn < alignment.getNPattern(); ++ptn) { int first_site = pattern_first_site[ptn]; int model_id = site_model[first_site]; - used_model[model_id] = true; + used_models++; + model_used[model_id] = true; + alignment.ptn_rate_mat.push_back(models_rate[model_id]); alignment.ptn_state_freq.push_back(models_freq[model_id]); - alignment.site_rate_matrices.insert( - alignment.site_rate_matrices.end(), - models_rate[model_id].begin(), - models_rate[model_id].end()); } for (size_t model_id = 0; model_id < models_freq.size(); ++model_id) { - if (!used_model[model_id]) + if (!model_used[model_id]) { delete[] models_freq[model_id]; + delete[] models_rate[model_id]; } } - cout << models_freq.size() << " distinct per-site state frequency vectors detected" << endl; + cout << used_models << " distinct per-site models detected" << endl; } void read_site_model_file(const std::string &filename, Alignment &alignment) @@ -348,8 +355,9 @@ void write_binary_site_model_file(const std::string &filename, Alignment &alignm { size_t nsites = alignment.getNSite(); size_t nstates = alignment.num_states; + size_t nrates = alignment.num_rates; ASSERT(nstates == 20); - + ASSERT(nrates == 190); try { ofstream out; @@ -368,13 +376,13 @@ void write_binary_site_model_file(const std::string &filename, Alignment &alignm for (size_t i = 0; i < nsites; ++i) { double *state_freq = alignment.ptn_state_freq[pattern_index[i]]; - out.write(reinterpret_cast(state_freq), 20 * sizeof(double)); + out.write(reinterpret_cast(state_freq), nstates * sizeof(double)); } for (size_t i = 0; i < nsites; ++i) { - double *rate_para_ptr = alignment.site_rate_matrices.data() + pattern_index[i] * 190; - out.write(reinterpret_cast(rate_para_ptr), 190 * sizeof(double)); + double *rate_mat = alignment.ptn_rate_mat[pattern_index[i]]; + out.write(reinterpret_cast(rate_mat), nrates * sizeof(double)); } cout << "Site mutsel model printed to " << filename << endl; @@ -388,15 +396,15 @@ void write_binary_site_model_file(const std::string &filename, Alignment &alignm DoubleVector computeMutselSiteRates(Alignment &alignment) { ASSERT(alignment.num_states == 20); + ASSERT(alignment.ptn_rate_mat.size() == alignment.getNPattern()); ASSERT(alignment.ptn_state_freq.size() == alignment.getNPattern()); - ASSERT(alignment.site_rate_matrices.size() == alignment.getNPattern() * 190); size_t npattern = alignment.getNPattern(); DoubleVector pattern_rates(npattern); for (size_t ptn = 0; ptn < npattern; ++ptn) { + double *R = alignment.ptn_rate_mat[ptn]; double *pi = alignment.ptn_state_freq[ptn]; - double *R = alignment.site_rate_matrices.data() + ptn * 190; double rate = 0.0; int idx = 0; for (int i = 0; i < 20; ++i) @@ -484,4 +492,4 @@ void computeMutselSiteFrequencyModel(Params ¶ms, Alignment *alignment) cout << endl << "===> CONTINUE ANALYSIS USING THE INFERRED MUTSEL MODEL" << endl; -} \ No newline at end of file +}