`-wslr`/`-wslm` `.sitelh` violates its own `sum(exp(LnLW_k)) == exp(LnL)` invariant for codon/morphology models

[famulare]

Summary

The per-category site-log-likelihood output (-wslr, -wslm, -wslmr → .sitelh)
documents this invariant in its own file header:

#   LnLW_k: Logarithm of (category-k site likelihood times category-k weight)
#           Thus, sum of exp(LnLW_k) is equal to exp(LnL)

On numerically-extreme sites the invariant is violated whenever the likelihood kernel
runs with per-category scaling but the .sitelh writer takes its per-pattern
scaling branch. The written LnLW_k come out systematically too low — by exactly one
scaling step, LOG_SCALING_THRESHOLD ≈ 177.4457 — so logsumexp_k(LnLW_k) falls far
below LnL. The error is one-signed (LnL too high relative to logsumexp_k(LnLW_k)).

This affects alphabets that are not 4-state DNA or 20-state protein — codon,
morphology, and other custom-state models. DNA and protein output is correct.
.siteprob (-wspr) is not affected (its per-site posteriors are normalized, so the
constant per-site shift cancels), so on the affected sites .sitelh and .siteprob are
mutually inconsistent in log space.

Version

IQ-TREE 2.4.0 (current master, tip a00094e0, 2025-04-10). Reproduced with the
official 2.4.0 release binary.

Root cause

The likelihood kernel decides its scaling mode via PhyloTree::safe_numeric, set in
PhyloTree::setLikelihoodKernel() (tree/phylotreesse.cpp:95):

safe_numeric = (params && (params->lk_safe_scaling || leafNum >= params->numseq_safe_scaling)) ||
    (aln && aln->num_states != 4 && aln->num_states != 20);

When safe_numeric is true the kernel scales per category and lays the scale counts
out per category (scale_num[ptn*ncat_mix + c]).

PhyloTree::computePatternLikelihood() (tree/phylotree.cpp, ~L1542 and ~L1568), which
produces the per-category values written to .sitelh, chose its per-category-vs-per-pattern
branch with a different, narrower condition that omits the num_states clause:

if (params->lk_safe_scaling || leafNum >= params->numseq_safe_scaling) {
    // per-category scaling — reads scale_num[ptn*ncat + c]
} else {
    // per-pattern scaling  — reads scale_num[ptn], one scale for all categories
}

For e.g. a 61-state codon model this condition is false while safe_numeric is true, so
the writer takes the per-pattern branch and reads a per-category scale_num array as flat
per-pattern. It mis-indexes the scale counts and applies a single, wrong scale step
uniformly to every category. On sites deep enough that the categories underflow at
different scale counts, every LnLW_k ends up exactly LOG_SCALING_THRESHOLD too low and
the invariant breaks.

When it triggers (and when it does not)

The two conditions agree except on the num_states clause, so the bug requires the writer
condition to be false while safe_numeric is true:

• Triggers: non-4-state-DNA / non-20-state-protein alphabet (codon, morphology,
  custom), run without -safe, on a tree below numseq_safe_scaling (default
  2000) tips.
• Does not trigger: DNA/protein (the clause is inert); or any run where -safe is set
  or the tree has ≥ numseq_safe_scaling tips — those make the writer's condition true as
  well, so it matches the kernel. (I.e. -safe and large trees mask this bug rather than
  cause it.)

Fingerprint (how to identify affected output)

• Deviation LnL − logsumexp_k(LnLW_k) is quantized to integer multiples of
  LOG_SCALING_THRESHOLD ≈ 177.4457, always positive (never a small rounding value).
• Confined to numerically-extreme sites where per-category scale counts diverge (one
  dominant category, others underflowed) — deep / long-branch sites.

Minimal reproduction (self-contained, no external data)

Simulate a codon alignment with widely-separated FreeRate categories and long branches (so
some sites underflow the slow categories), then re-run with the true model and -wslr:

IQTREE=iqtree2   # 2.4.0

# 16-taxon codon alignment, deep branches, 5 well-separated rate categories
$IQTREE --alisim sim -t "RANDOM{yh/16}" \
  -m "MG{2.0}+F3X4+R5{0.4/0.001/0.3/0.05/0.2/1.0/0.08/8.0/0.02/60.0}" \
  --length 2001 --seqtype CODON --branch-scale 3.0 -seed 7

$IQTREE -s sim.phy -te sim.treefile \
  -m "MG{2.0}+F3X4+R5{0.4/0.001/0.3/0.05/0.2/1.0/0.08/8.0/0.02/60.0}" \
  -st CODON -wslr -seed 7 -pre chk

Then check each site's logsumexp_k(LnLW_k) against LnL with this script:

#!/usr/bin/env python3
"""Check the -wslr .sitelh invariant: logsumexp_k(LnLW_k) == LnL per site."""
import sys, math

def logsumexp(xs):
    xs = [x for x in xs if x != float('-inf') and not math.isnan(x)]
    if not xs:
        return float('-inf')
    m = max(xs)
    return m + math.log(sum(math.exp(x - m) for x in xs))

def parse(path):
    rows = []
    with open(path) as f:
        for line in f:
            line = line.rstrip('\n')
            if not line or line.startswith('#'):
                continue
            parts = line.split('\t')
            if parts[0] in ('Site', 'Part'):
                header = parts
                continue
            rows.append(parts)
    return header, rows

def main(path):
    header, rows = parse(path)
    has_part = header[0] == 'Part'
    lnl_idx = header.index('LnL')
    cat_start = lnl_idx + 1
    maxdev = 0.0
    nviol = 0
    worst = []
    for parts in rows:
        try:
            lnl = float(parts[lnl_idx])
        except (ValueError, IndexError):
            continue
        cats = []
        for p in parts[cat_start:]:
            if p == '' or p is None:
                continue
            try:
                v = float(p)
            except ValueError:
                continue
            cats.append(v)
        lse = logsumexp(cats)
        dev = lnl - lse  # positive => LnL too high (invariant violated one-signed)
        if abs(dev) > abs(maxdev):
            maxdev = dev
        if abs(dev) > 1e-3:
            nviol += 1
            site = parts[1] if has_part else parts[0]
            part = parts[0] if has_part else '-'
            srt = sorted(cats, reverse=True)
            gap = (srt[0] - srt[1]) if len(srt) > 1 else float('nan')
            worst.append((abs(dev), part, site, lnl, lse, dev, gap, cats))
    worst.sort(reverse=True)
    print(f"file: {path}")
    print(f"rows: {len(rows)}  violations(|dev|>1e-3): {nviol}  max|dev|: {maxdev:.6f}")
    print(f"{'part':>4} {'site':>6} {'LnL':>12} {'logsumexp':>12} {'dev':>12} {'top-gap':>10}")
    for _, part, site, lnl, lse, dev, gap, cats in worst[:15]:
        print(f"{part:>4} {site:>6} {lnl:12.6f} {lse:12.6f} {dev:12.6f} {gap:10.3f}")
    return nviol, maxdev

if __name__ == '__main__':
    main(sys.argv[1])

On the buggy 2.4.0 build:

$ python3 check_sitelh_invariant.py chk.sitelh
rows: 667  violations(|dev|>1e-3): 18  max|dev|: 177.445710
part   site          LnL    logsumexp          dev    top-gap
   -    560   -56.767000  -234.212710   177.445710      9.318
   -    563   -74.797800  -252.243500   177.445700     32.959
   -     42   -11.026200  -188.471877   177.445677      1.973
   ...

Example violating row (sim.phy site 85), 2.4.0:

Site   LnL       LnLW_1     LnLW_2     LnLW_3     LnLW_4     LnLW_5
85    -3.8118   -181.9381  -182.2600  -183.3261  -188.9624  -215.2784

LnL = -3.8118, but logsumexp(LnLW_k) ≈ -181.26 — every LnLW_k is 177.4457 too low.
The same site's .siteprob is correct: 0.506 0.367 0.126 0.00045 1.7e-15 (sums to 1).

Fix

Select the writer's branch with the tree's own safe_numeric flag — the authoritative
record of the scaling layout the kernel actually used — instead of re-deriving a narrower
condition. Two-line change in PhyloTree::computePatternLikelihood() (both branches).
DNA/protein output is byte-identical; codon/morphology now satisfy the invariant to print
precision. PR to follow.

---

AI-assistance disclosure: this report was diagnosed and drafted with the help
of Claude Opus 4.8. All findings were reviewed and verified by Mike Famulare
before submission.