Skip to content

Fix MutSel potential bugs related to alignment copying #211

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
bqminh merged 1 commit into from
Sep 9, 2026
Merged

Fix MutSel potential bugs related to alignment copying #211

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
109 changes: 80 additions & 29 deletions alignment/alignment.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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()
}

Expand All @@ -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<double *>::reverse_iterator rit = ptn_rate_mat.rbegin(); rit != ptn_rate_mat.rend(); ++rit) {
delete [] (*rit);
(*rit) = nullptr;
}
ptn_rate_mat.clear();
for (vector<double *>::reverse_iterator rit = ptn_state_freq.rbegin(); rit != ptn_state_freq.rend(); ++rit) {
delete [] (*rit);
(*rit) = nullptr;
Expand Down Expand Up @@ -1589,6 +1595,13 @@ void Alignment::regroupSitePattern(const IntVector &site_group) {
}
}
// refill the existing pattern-specific parameters
if (isSSM()) {
vector<double*> 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<double*> stored_ptn_state_freq = ptn_state_freq;
ptn_state_freq.clear();
Expand Down Expand Up @@ -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) {
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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();
Expand All @@ -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);
}
}
}
Expand All @@ -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();
Expand Down Expand Up @@ -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;
Expand All @@ -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");
}
Expand All @@ -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])++;
Expand Down Expand Up @@ -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());
Expand Down
15 changes: 9 additions & 6 deletions alignment/alignment.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -552,6 +552,7 @@ class Alignment : public vector<Pattern>, 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; }
Expand Down Expand Up @@ -944,15 +945,17 @@ class Alignment : public vector<Pattern>, public CharSet, public StateSpace {
vector<uint32_t> 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<double*> ptn_rate_mat;

/** pattern ID to state frequency vector map */
vector<double*> ptn_state_freq;

/** site to rate matrix. Stored in row-major order [num_sites, 190] used for MUTSEL */
vector<double> 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
*/
Expand Down
21 changes: 7 additions & 14 deletions model/modelfactory.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -639,8 +639,7 @@ ModelFactory::ModelFactory(Params &params, 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");
Expand All @@ -660,29 +659,26 @@ ModelFactory::ModelFactory(Params &params, 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}
Expand All @@ -691,11 +687,8 @@ ModelFactory::ModelFactory(Params &params, string &model_name, PhyloTree *tree,
modeli->fixParameters(true);
models->push_back(modeli);
}

models->joinEigenMemory();
models->decomposeRateMatrix();

model = models;
}
}

Expand Down
Loading
Loading