PMSR model implementation

[StefanFlaumberg]

This pull request adapts the iqtree/iqtree2#391 PR for IQ-Tree 3.

Commits of the original PR were directly transferred, with only minor changes introduced for the sake of compatibility and small bug fixing.

[bqminh]

Claude Code:

PR #40 — PMSR (Posterior Mean Site Rate) Model

Status: Draft | 2,386 additions / 1,475 deletions across 42 files

What It Does

Extends the existing PMSF (Posterior Mean Site Frequency) workflow to also support per-site rates. Three new modes:

1. -rs — load pre-computed per-site rates from a .siterate file
2. -rt <guide_tree> — infer per-site rates from a guide tree, then fix them for the main analysis
3. -frt — combine PMSF + PMSR simultaneously (per-site frequencies and per-site rates)

The core change is that Alignment gains a ptn_rate_scaler vector (one rate multiplier per pattern) alongside the existing ptn_state_freq, and ModelSet is extended to use both.

Key Changes

• alignment/alignment.h/.cpp: site_state_freq renamed to ptn_state_freq (pattern-indexed, not site-indexed — correct design), ptn_rate_scaler added, readSiteStateFreq() replaced by generic readSiteParamFile(file, "freq"|"rate")
• model/modelset.cpp: Checkpoint methods (setCheckpoint, saveCheckpoint, restoreCheckpoint) were missing — now added. This is a genuine bug fix that benefits all SSF users.
• model/rateheterogeneity.cpp: Constructor now requires PhyloTree* to detect SSR and set name = "+SSR"
• tree/phylotree.cpp: computePatternRate() rewritten to correctly handle +I offset and posterior mean/max modes
• tree/iqtree.cpp: BNNI bug fix — tree topology now established before model copy in bootstrap NNI loop
• main/phyloanalysis.cpp: computeSiteFrequencyModel() → computeSiteSpecificModel() extended for both freq and rate
• main/treetesting.cpp: printSiteRate() output now includes CatPP column (posterior probability of best category)

Issues That Should Block Merge

1. Memory aliasing / double-free in ptn_state_freq
ungroupSitePattern() and regroupSitePattern() copy the vector<double*> by value without deep-copying the pointed-to arrays. Multiple entries can point to the same buffer. When the old alignment is destroyed, those pointers are freed; the new alignment then double-frees them on destruction.

2. ptn_rate_scaler not checkpointed
ptn_rate_scaler lives on Alignment, which has no checkpoint serialisation for it. If a long run is interrupted and resumed from checkpoint, the site-specific rates are silently lost and must be re-read from the original file (which may not be available). Should be documented as a known limitation at minimum.

3. Tree mutation inside ModelSet constructor
normalizePtnRateScaler() is called inside the constructor and mutates the tree's branch lengths as a side effect. This is architecturally wrong — constructors should not mutate external objects. It also means restoreCheckpoint() does not redo the normalisation, leaving branch lengths in an inconsistent state after resume.

4. Partition model + SSR — no guard
--partition -rs is not explicitly blocked. The code reads site rates for the concatenated super-alignment, not per-partition, producing silently wrong results.

5. RateHeterogeneity(nullptr) null dereference
The new constructor evaluates phylo_tree->aln->isSSR() unconditionally. Any call site that passes nullptr will crash. A null check is needed.

Issues Worth Addressing

• Stringly-typed dispatch: readSiteParamFile(file, "freq"|"rate") should use an enum class SiteParamType instead of a string to catch typos at compile time.
• isSSF() / isSSR() lack const: prevents use on const Alignment& references.
• Non-deterministic tie-breaking in computePatternRate: random_double() < 0.5 for equal-likelihood categories produces non-reproducible output — worth documenting.
• Pairwise ML distances don't fully account for SSR (acknowledged limitation from iqtree2 PR #391 — should be documented).
• Comments in modelfactory.cpp still say "site-specific state frequency model" in places that now cover both SSF and SSR.

Positives

• The ModelSet checkpoint fix is a genuine bug fix benefiting all existing SSF users, independent of PMSR.
• Renaming from site-indexed to pattern-indexed storage is the correct design direction.
• The BNNI ordering bug fix in iqtree.cpp is correct.
• The CatPP column addition to .siterates output is a useful improvement.
• Scientific value is clear — the feature has been requested and positively received.

Verdict

Not ready to merge in current draft state. The core logic is correct, but the five structural issues (memory aliasing, missing checkpoint, constructor side-effects, partition incompatibility, null dereference) need to be resolved first. The 42-file scope and architectural changes to Alignment, ModelSet, and RateHeterogeneity warrant careful review of memory management and checkpoint correctness paths.

[StefanFlaumberg]

Dear @bqminh,

Thank you for the neuro-review) I'll try to use this tech myself.

Apart from the issues mentioned here, there've been some new changes I'd like to include, and the commits evidently fall behind the current IQ-Tree version quite significantly, which would require major merging or rebasing.
I'm going to tackle these problems in two steps:
First, as soon as I find time, I will include new changes and resolve the issues from the review. I'll write a short summary of what will have been made.
Second, upon receiving initial approval, I will readily make the update to comply with the master version, so that merging the PR becomes possible. Since it's quite hard to constantly follow all the changes happening to the master branch, I would like to make the update only once, right before everyone is ready to merge the PR.

Hope we will finish this PR soon!

[bqminh]

Thank you. And much sorry for not processing your pull request earlier, so that you now have to merge master branch into it. It's not easy to review a big PR, but AI now helps me to review it. When merging master, if there are conflicts that you can't resolve, please let me know. In that case we can have a coding session to resolve it together. Or I can also try to resolve it, if I have permission to write to your branch.

The IQ-TREE team is now extensively using Claude Code to speed up development significantly (say months can go down to a week or two). I highly recommend it.

[StefanFlaumberg]

@bqminh,
Alright!

A slight digression on code formatting:
It is evident the current code has a totally free and piecewise formatting style: countless extra or missing indentations, some files using tabs, some using spaces, yet most using both of them, lines as long as one could wish, and various other bad practices. On the one hand, it is always a pain to look at and asks for fixing, but, on the other hand, mixing format and actual code fixing in a single commit significantly decreases readability of the proposed changes, which encourages to leave the formatting alone, as bad as it is.
That said, have you ever considered using code auto-formatting tools, such as clang-format? It would take only creating a simple format-definition file, a single commit to fix all the existing code, and, well, reformatting of the changes before publishing every new commit. In turn, it would give a clear code that is easier to read and develop.

[bqminh]

@trongnhanuit @thomaskf Could this indentation be caused by the softwipe refactoring and different IDE settings?