[alpae]

estimate_pam (and with it estimate_many_pam) now uses a new C function EstimatePamEigen instead of EstimatePam:

• the aligned pairs of symbols are counted once, so scores and likelihoods are sums over the distinct pairs (at most 26²) instead of over all alignment positions;
• logPAM1 is the rate matrix of a reversible process, so it is decomposed once per MutipleAlEnv (symmetrized with its stationary distribution, then eigh). In the Newton iterations of the variance, exp(pam·logPAM1) and its first two derivatives are then evaluated from the eigenvalues, only for the pairs that occur, instead of with matrix exponentials and products.

Results

• Score and PAM distance are bit-identical to EstimatePam, and the variance agrees to ≤ 1e-9 relative (worst case at the PAM 1000 saturation). test_darwin passes unchanged.
• Speed: estimate_pam takes 11 µs instead of 378 µs per pair (300 aa, 30 % substitutions). estimate_many_pam on 100 sequences with 1 thread takes 0.09 s instead of 1.47 s.

Fallback

If a custom log_pam1 is not reversible (its symmetrized form isn't symmetric) or has no positive stationary distribution, MutipleAlEnv falls back to the original EstimatePam. It stays available as the private _estimate_pam_reference for tests.

Tests

test/test_estimate_pam.py compares the new code with the reference on 300 random alignments (10 to 1000 positions, all levels of divergence and gaps), on edge cases (identical sequences, saturated, gaps only, X, B/Z), and through estimate_many_pam. It also checks the fallback with a perturbed non-reversible logPAM1.