From a64577e0d77664c97483aec4db9fc110ef6a831e Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Thu, 29 Oct 2026 21:39:39 +0100 Subject: [PATCH 1/2] Estimate PAM distances with an eigendecomposition of logPAM1 EstimatePam spends most of its time computing exp(pam * logPAM1) and its derivatives with matrix exponentials and products in every Newton iteration of the variance estimation, and summing scores over all alignment positions. EstimatePamEigen counts the distinct aligned pairs of symbols once and evaluates scores and likelihoods over these pairs. logPAM1 is the rate matrix of a reversible process, so it is decomposed once per MutipleAlEnv (symmetrized with its stationary distribution, then eigh), and the needed entries of exp(pam * logPAM1) and its derivatives are computed from the eigenvalues only for the pairs that occur. Scores and PAM distances are identical, variances agree to about 1e-9 relative. estimate_pam is about 30x faster per pair. If a custom logPAM1 is not reversible, estimate_pam and estimate_many_pam fall back to the original EstimatePam, which stays available as _estimate_pam_reference for tests. --- MANIFEST.in | 4 +- c_source/EstimatePamEigen.c | 266 ++++++++++++++++++++++++++++++++++++ c_source/EstimatePamEigen.h | 26 ++++ pyopa/backend/cpyopa.pxd | 14 ++ pyopa/backend/pyopa.pyx | 113 +++++++++++++-- setup.py | 2 +- test/test_estimate_pam.py | 91 ++++++++++++ 7 files changed, 498 insertions(+), 18 deletions(-) create mode 100644 c_source/EstimatePamEigen.c create mode 100644 c_source/EstimatePamEigen.h create mode 100644 test/test_estimate_pam.py diff --git a/MANIFEST.in b/MANIFEST.in index bf15659..7cb787f 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1,9 +1,9 @@ include setup.py requirements.txt pyproject.toml README.rst CHANGELOG.md include c_source/Python_extension.c c_source/DynProgr_sse_short.c c_source/DynProgr_sse_byte.c -include c_source/DynProgr_scalar.c c_source/EstimatePam.c c_source/Page_size.c c_source/DynProgr_sse_double.c +include c_source/DynProgr_scalar.c c_source/EstimatePam.c c_source/EstimatePamEigen.c c_source/Page_size.c c_source/DynProgr_sse_double.c include c_source/Python_extension.h c_source/DynProgr_sse_short.h c_source/DynProgr_sse_byte.h include c_source/DynProgr_scalar.h c_source/matrix.h c_source/debug.h c_source/swps3.h c_source/EstimatePam.h -include c_source/Page_size.h c_source/DynProgr_sse_double.h +include c_source/Page_size.h c_source/DynProgr_sse_double.h c_source/EstimatePamEigen.h graft c_source/simde include pyopa/backend/cpyopa.pxd pyopa/backend/pyopa.pyx include pyopa/matrices/* diff --git a/c_source/EstimatePamEigen.c b/c_source/EstimatePamEigen.c new file mode 100644 index 0000000..5362dfa --- /dev/null +++ b/c_source/EstimatePamEigen.c @@ -0,0 +1,266 @@ +#include +#include "EstimatePamEigen.h" + +/* + # EstimatePam (see EstimatePam.c) with two changes that make it faster: + # - the aligned pairs of symbols are counted, so scores and likelihoods + # are sums over the (at most MATRIX_DIM^2) distinct pairs instead of + # over all positions of the alignment + # - exp(p*logPAM1) and its derivatives are evaluated from the + # eigendecomposition of logPAM1, only for the pairs that occur, + # instead of computing matrix exponentials and products in every + # Newton iteration + */ + +#define NPAIRS (MATRIX_DIM * MATRIX_DIM) + +typedef struct PairCounts { + int delF, /* number of indels */ + delI, /* number of amino acids deleted - delF */ + nmatch, /* number of aligned pairs of symbols */ + npairs, /* number of distinct pairs */ + idx[NPAIRS]; /* distinct pairs, as index into a Dayhoff matrix */ + double cnt[NPAIRS]; /* number of occurrences of each pair */ +} PairCounts; + +static void CountPairs(const char *o1, const char *o2, int len, PairCounts *pc) { + int counts[NPAIRS] = { 0 }; + int i, ind; + + pc->delF = pc->delI = pc->nmatch = pc->npairs = 0; + for (ind = i = 0; i < len; i++) { + if (o1[i] == '_') { + if (o2[i] == '_') + continue; /* do not count del==del */ + if (ind == 1) + pc->delI++; + else { + pc->delF++; + ind = 1; + } + } else if (o2[i] == '_') { + if (ind == 2) + pc->delI++; + else { + pc->delF++; + ind = 2; + } + } else { + if (o1[i] != 'X' && o2[i] != 'X') { + counts[o1[i] * MATRIX_DIM + o2[i]]++; + pc->nmatch++; + } + ind = 0; + } + } + for (i = 0; i < NPAIRS; i++) + if (counts[i] > 0) { + pc->idx[pc->npairs] = i; + pc->cnt[pc->npairs++] = counts[i]; + } +} + +static double ComputeScorePairs(int idms, const PairCounts *pc, const DayMatrix *DMS) { + int k; + const double *Simi = DMS[idms].Simi; + double cj = pc->delF * DMS[idms].FixedDel + pc->delI * DMS[idms].IncDel; + + for (k = 0; k < pc->npairs; k++) + cj += pc->cnt[k] * Simi[pc->idx[k]]; + return cj; +} + +/* Multi-point multiple maximum searching a la Brent, as FindMaxPamR */ +static int FindMaxPamPairs(int x[], double f[], int n, double *res, + const PairCounts *pc, const DayMatrix *DMS) { + int j, j1, j2, x2[3 * MAXPOINTS + 1]; + double cj, c1, c2, f2[3 * MAXPOINTS + 1]; + + restart: + /* end condition */ + for (j = 1; j < n && x[j] - x[j - 1] <= 1; j++) + ; + if (j >= n) { + for (j1 = j = 0; j < n; j++) + if (f[j] > f[j1]) + j1 = j; + *res = f[j1]; + return (x[j1]); + } + + /* can insert new points */ + if (n < MAXPOINTS) { + for (j1 = j = 1; j < n; j++) + if (x[j] - x[j - 1] > x[j1] - x[j1 - 1]) + j1 = j; + j = (x[j1] + x[j1 - 1]) / 2; + cj = ComputeScorePairs(j, pc, DMS); + for (j2 = n - 1; j2 >= j1; j2--) { + f[j2 + 1] = f[j2]; + x[j2 + 1] = x[j2]; + } + x[j1] = j; + f[j1] = cj; + n++; + goto restart; + } + + /* if a valley is found, split into two searches */ + for (j = 1; j < n - 1; j++) + if (f[j] < f[j - 1] && f[j] < f[j + 1]) { + for (j2 = j; j2 < n; j2++) { + x2[j2 - j] = x[j2]; + f2[j2 - j] = f[j2]; + } + j2 = FindMaxPamPairs(x2, f2, n - j, &c2, pc, DMS); + j1 = FindMaxPamPairs(x, f, j + 1, &c1, pc, DMS); + if (c1 > c2) { + *res = c1; + return (j1); + } else { + *res = c2; + return (j2); + } + } + + /* discard end farthest from maximum */ + for (j1 = j = 0; j < n; j++) + if (f[j] > f[j1]) + j1 = j; + if (x[j1] - x[0] > x[n - 1] - x[j1]) + for (j = 1; j < n; j++) { + f[j - 1] = f[j]; + x[j - 1] = x[j]; + } + n--; + goto restart; +} + +/* PamVariance by maximum likelihood, as ComputeVarianceML */ +static double VarianceML(const PairCounts *pc, int dmsnumb, const DayMatrix *DMS, + int DMSLen, const PamModel *m) { + double e[MATRIX_DIM], le[MATRIX_DIM], l2e[MATRIX_DIM]; + double FD1, FD1d, p, logL1, logL2, incr, x1, Mp, logMp, log2Mp, vw, rmpam; + const double *Va, *Wb; + int j, k, l; + + for (j = DMSLen; j > 3 && DMS[j - 1].FixedDel >= -1; j--) + ; + FD1 = log(10.0) * (DMS[1].FixedDel - DMS[j - 1].FixedDel) + / log(DMS[1].PamNumber / DMS[j - 1].PamNumber); + FD1d = pc->delF * FD1 / 10; + p = DMS[dmsnumb].PamNumber; + + for (;;) { + /* exp(p*logPAM1), logPAM1*exp(p*logPAM1) and logPAM1^2*exp(p*logPAM1) + are V diag(f(lam)) W with f = exp(p*lam), lam*exp(p*lam), lam^2*exp(p*lam) */ + for (l = 0; l < m->nlam; l++) { + e[l] = exp(p * m->lam[l]); + le[l] = m->lam[l] * e[l]; + l2e[l] = m->lam[l] * le[l]; + } + logL1 = FD1d / p; + logL2 = -FD1d / (p * p); + for (k = 0; k < pc->npairs; k++) { + Va = m->V + (pc->idx[k] / MATRIX_DIM) * MATRIX_DIM; + Wb = m->W + pc->idx[k] % MATRIX_DIM; + Mp = logMp = log2Mp = 0; + for (l = 0; l < m->nlam; l++) { + vw = Va[l] * Wb[l * MATRIX_DIM]; + Mp += vw * e[l]; + logMp += vw * le[l]; + log2Mp += vw * l2e[l]; + } + /* symbols not in the model (e.g. X), ignore entry */ + if (Mp == 0 && logMp == 0 && log2Mp == 0) + continue; + x1 = logMp / Mp; + logL1 += pc->cnt[k] * x1; + logL2 += pc->cnt[k] * (log2Mp / Mp - x1 * x1); + } + + if (logL2 >= 0) + break; + incr = -logL1 / logL2; + /* do not go beyond the limit of the DMS array */ + if (p + incr > DMS[DMSLen - 1].PamNumber) + incr = DMS[DMSLen - 1].PamNumber - p; + x1 = 0.04 * (p + 25); + if (fabs(incr) > x1) + incr = (incr < 0) ? -x1 : x1; + p += incr; + if (fabs(incr) < 1e-8 * p) + break; + } + if (logL2 >= 0 && p >= DMS[DMSLen - 1].PamNumber) { + /* at the end of the DMS, estimate the variance by the idealized + mutations model */ + p = DMS[DMSLen - 1].PamNumber; + rmpam = pow(1 - 0.20 / 19, -p); + logL2 = -pc->nmatch * 19 * log(1 - 0.20 / 19) * log(1 - 0.20 / 19) + / ((rmpam - 1) * (rmpam + 19)); + } + return -1 / logL2; +} + +/* PamVariance for very small pam, as ComputeSmallVariance */ +static double VarianceSmall(const PairCounts *pc, int dmsnumb, const DayMatrix *DMS, + const PamModel *m) { + double a1, a2, b0, c, t1, t2, t3, t4, t5, var; + const double *L1 = m->L1, *L2 = m->L2, *L3 = m->L3; + int ij, k, n = MATRIX_DIM; + + a1 = a2 = 0; + b0 = pc->delF * 0.7434; + for (k = 0; k < pc->npairs; k++) { + ij = pc->idx[k]; + c = pc->cnt[k]; + if (L1[ij] == 0) + continue; + if (ij % (n + 1) == 0) { + a1 += c * L1[ij]; + a2 += c * 0.5 * (L2[ij] - L1[ij] * L1[ij]); + } else { + b0 += c; + a1 += c * 0.5 * L2[ij] / L1[ij]; + a2 += c * (L3[ij] / 6 - L2[ij] * L2[ij] / 8 / L1[ij]) / L1[ij]; + } + } + + t5 = b0 * b0; + t3 = a1 * a1; + t2 = 2.0 * t3 + a2 * t5; + t4 = b0 + 1.0; + t1 = (3.0 * b0 + 2.0) * a2 + t2; + var = (4.0 * t3 * t3 + + ((32.0 + 4.0 * t5 + 24.0 * b0) * t3 + + (28.0 * b0 + 12.0 + (23.0 + 8.0 * b0 + t5) * t5) * a2) + * a2) / t3 * t4 / (t1 * t1); + if (a1 == 0 || t1 == 0) + var = DMS[dmsnumb].PamNumber * DMS[dmsnumb].PamNumber; + return var; +} + +void EstimatePamEigen(const char* o1, const char* o2, int len, DayMatrix* DMS, + int DMSLen, const PamModel* model, double* result) { + PairCounts pc; + double f[3 * MAXPOINTS + 1], sim; + int k, x[3 * MAXPOINTS + 1], dmsnumb; + + /* Increase the length to be compatible with the core */ + DMSLen += 1; + + CountPairs(o1, o2, len, &pc); + for (k = 0; k <= 3 * MAXPOINTS; k++) { + x[k] = k * (DMSLen - 2) / (3 * MAXPOINTS) + 1; + f[k] = ComputeScorePairs(x[k], &pc, DMS); + } + dmsnumb = FindMaxPamPairs(x, f, 3 * MAXPOINTS + 1, &sim, &pc, DMS); + + result[0] = ComputeScorePairs(dmsnumb, &pc, DMS); + result[1] = DMS[dmsnumb].PamNumber; + if (dmsnumb <= 1) + result[2] = VarianceSmall(&pc, dmsnumb, DMS, model); + else + result[2] = VarianceML(&pc, dmsnumb, DMS, DMSLen, model); +} diff --git a/c_source/EstimatePamEigen.h b/c_source/EstimatePamEigen.h new file mode 100644 index 0000000..8c187b0 --- /dev/null +++ b/c_source/EstimatePamEigen.h @@ -0,0 +1,26 @@ +#ifndef ESTIMATEPAMEIGEN_H_ +#define ESTIMATEPAMEIGEN_H_ + +#include "EstimatePam.h" + +/* logPAM1 and the data derived from it that EstimatePamEigen needs. It only + depends on logPAM1, so it is computed once (see MutipleAlEnv in pyopa.pyx). + logPAM1 = V diag(lam) W with W = V^-1, restricted to the nlam letters for + which logPAM1 has non-zero rows/columns; V and W are zero elsewhere. */ +typedef struct PamModel { + int nlam; + double lam[MATRIX_DIM]; + double V[MATRIX_DIM * MATRIX_DIM]; /* V[a * MATRIX_DIM + k] */ + double W[MATRIX_DIM * MATRIX_DIM]; /* W[k * MATRIX_DIM + b] */ + double L1[MATRIX_DIM * MATRIX_DIM]; /* logPAM1 */ + double L2[MATRIX_DIM * MATRIX_DIM]; /* logPAM1^2 */ + double L3[MATRIX_DIM * MATRIX_DIM]; /* logPAM1^3 */ +} PamModel; + +/* Same results as EstimatePam, but evaluates exp(pam * logPAM1) and its + derivatives from the eigendecomposition, only for the pairs of symbols + occurring in the alignment. */ +void EstimatePamEigen(const char* o1, const char* o2, int len, DayMatrix* DMS, + int DMSLen, const PamModel* model, double* result); + +#endif /* ESTIMATEPAMEIGEN_H_ */ diff --git a/pyopa/backend/cpyopa.pxd b/pyopa/backend/cpyopa.pxd index 4ccad55..e8971ad 100644 --- a/pyopa/backend/cpyopa.pxd +++ b/pyopa/backend/cpyopa.pxd @@ -37,6 +37,20 @@ cdef extern from "EstimatePam.h": void CreateOrigDayMatrix(double* log_pam1, double pam, double* new_matrix) +cdef extern from "EstimatePamEigen.h": + ctypedef struct PamModel: + int nlam + double lam[26] + double V[676] + double W[676] + double L1[676] + double L2[676] + double L3[676] + + void EstimatePamEigen(const char* o1, const char* o2, int len, DayMatrix* DMS, int DMSLen, + const PamModel* model, double* result) nogil + + cdef extern from "DynProgr_sse_double.h": ctypedef struct ProfileDouble: pass diff --git a/pyopa/backend/pyopa.pyx b/pyopa/backend/pyopa.pyx index 614bac2..8d3ea8c 100644 --- a/pyopa/backend/pyopa.pyx +++ b/pyopa/backend/pyopa.pyx @@ -734,10 +734,60 @@ cdef class MutipleAlEnv: # number of estimate_many_pam calls currently using _c_dayMatrices without the GIL cdef int _active_estimates cdef log_pam1 + cdef cpyopa.PamModel _pam_model + # whether estimate_pam can use the eigendecomposition in _pam_model + cdef bint _use_eigen def __init__(self, envs, log_pam1): self.create_day_matrices(envs) self.log_pam1 = log_pam1 + self._use_eigen = self._create_pam_model(log_pam1.float64_matrix) + + def _create_pam_model(self, Q): + """ + Computes the eigendecomposition logPAM1 = V diag(lam) W, W = V^-1, used by + EstimatePamEigen. logPAM1 is the rate matrix of a reversible process, so it is + symmetrized with its stationary distribution and decomposed with eigh, restricted + to the symbols for which it has non-zero rows/columns. + + :return: False if logPAM1 is not reversible, in which case estimate_pam falls back + to the slower EstimatePam + """ + Q = np.asarray(Q, dtype=np.float64) + active = np.flatnonzero(np.abs(Q).sum(axis=0) + np.abs(Q).sum(axis=1) > 0) + if len(active) == 0: + return False + A = Q[np.ix_(active, active)] + # stationary distribution pi: A pi = 0 + w, v = np.linalg.eig(A) + pi = np.real(v[:, np.argmin(np.abs(w))]) + pi /= pi.sum() + if not np.all(pi > 0): + return False + s = np.sqrt(pi) + S = A * s[np.newaxis, :] / s[:, np.newaxis] + # S is symmetric iff logPAM1 is reversible (detailed balance) + if not np.allclose(S, S.T, rtol=1e-8, atol=1e-12 * np.abs(S).max()): + return False + lam, U = np.linalg.eigh((S + S.T) / 2) + + V = np.zeros((26, 26)) + W = np.zeros((26, 26)) + V[np.ix_(active, np.arange(len(active)))] = U * s[:, np.newaxis] + W[np.ix_(np.arange(len(active)), active)] = U.T / s[np.newaxis, :] + Q2 = Q @ Q + Q3 = Q @ Q2 + + self._pam_model.nlam = len(active) + for i in range(26): + self._pam_model.lam[i] = lam[i] if i < len(active) else 0.0 + for j in range(26): + self._pam_model.V[i * 26 + j] = V[i, j] + self._pam_model.W[i * 26 + j] = W[i, j] + self._pam_model.L1[i * 26 + j] = Q[i, j] + self._pam_model.L2[i * 26 + j] = Q2[i, j] + self._pam_model.L3[i * 26 + j] = Q3[i, j] + return True def __cinit__(self): self._c_dayMatrices = NULL @@ -765,6 +815,31 @@ cdef class MutipleAlEnv: pam_dist.data, matrices.data, self.dms_len) def estimate_pam(self, s1, s2): + """ + Estimates the PAM distance of two aligned sequences. + + :return: list of (score, pam distance, variance) + """ + if not self._use_eigen: + return self._estimate_pam_reference(s1, s2) + + if len(s1) != len(s2): + raise ValueError('The length of s1 and s2 must be equal!') + + if self._c_dayMatrices is NULL: + raise RuntimeError('The day matrices have been freed') + + cdef double result[3] + cpyopa.EstimatePamEigen(s1.s_norm.encode('utf-8'), s2.s_norm.encode('utf-8'), len(s1), + self._c_dayMatrices, self.dms_len, &self._pam_model, result) + return [result[0], result[1], result[2]] + + def _estimate_pam_reference(self, s1, s2): + """ + estimate_pam with the original EstimatePam, which computes matrix exponentials + of logPAM1 instead of using its eigendecomposition. Used as a fallback if logPAM1 + is not reversible, and as a reference in tests. + """ if len(s1) != len(s2): raise ValueError('The length of s1 and s2 must be equal!') @@ -776,12 +851,7 @@ cdef class MutipleAlEnv: cpyopa.EstimatePam(s1.s_norm.encode('utf-8'), s2.s_norm.encode('utf-8'), len(s1), self._c_dayMatrices, self.dms_len, logpam.data, result) - ret = [0.0] * 3 - ret[0] = result[0] - ret[1] = result[1] - ret[2] = result[2] - - return ret + return [result[0], result[1], result[2]] @cython.cdivision(True) def estimate_many_pam(self, seqs, num_threads=0): @@ -816,6 +886,8 @@ cdef class MutipleAlEnv: cdef cpyopa.DayMatrix* c_dayMatrices = self._c_dayMatrices cdef int dms_len = self.dms_len + cdef cpyopa.PamModel* pam_model = &self._pam_model + cdef bint use_eigen = self._use_eigen cdef int nthreads = num_threads seqs_bytes = [s.s_norm.encode('utf-8') for s in seqs] @@ -838,15 +910,26 @@ cdef class MutipleAlEnv: for seq_j in range(seq_i + 1, n): # condensed index of pair (seq_i, seq_j) k = seq_i * n - seq_i * (seq_i + 1) / 2 + seq_j - seq_i - 1 - cpyopa.EstimatePam( - c_seqs[seq_i], - c_seqs[seq_j], - seq_len, - c_dayMatrices, - dms_len, - logpam_ptr, - &result_ptr[3 * k] - ) + if use_eigen: + cpyopa.EstimatePamEigen( + c_seqs[seq_i], + c_seqs[seq_j], + seq_len, + c_dayMatrices, + dms_len, + pam_model, + &result_ptr[3 * k] + ) + else: + cpyopa.EstimatePam( + c_seqs[seq_i], + c_seqs[seq_j], + seq_len, + c_dayMatrices, + dms_len, + logpam_ptr, + &result_ptr[3 * k] + ) finally: self._active_estimates -= 1 free(c_seqs) diff --git a/setup.py b/setup.py index 9f95b36..64fad71 100644 --- a/setup.py +++ b/setup.py @@ -15,7 +15,7 @@ c_dir = 'c_source/' c_sources = ['Python_extension.c', 'DynProgr_sse_short.c', 'DynProgr_sse_byte.c', - 'DynProgr_scalar.c', 'EstimatePam.c', 'Page_size.c', 'DynProgr_sse_double.c'] + 'DynProgr_scalar.c', 'EstimatePam.c', 'EstimatePamEigen.c', 'Page_size.c', 'DynProgr_sse_double.c'] def openmp_flags(compiler_type): diff --git a/test/test_estimate_pam.py b/test/test_estimate_pam.py new file mode 100644 index 0000000..ed3a121 --- /dev/null +++ b/test/test_estimate_pam.py @@ -0,0 +1,91 @@ +import random +import types +import unittest + +import numpy as np +import pyopa + +AMINO_ACIDS = 'ACDEFGHIKLMNPQRSTVWY' + + +def random_alignment(rng, length, subst_rate, gap_rate): + """returns two aligned rows of equal length with random substitutions and gaps""" + a = ''.join(rng.choice(AMINO_ACIDS) for _ in range(length)) + b = ''.join('_' if rng.random() < gap_rate else + (rng.choice(AMINO_ACIDS) if rng.random() < subst_rate else c) for c in a) + a = ''.join('_' if rng.random() < gap_rate / 2 else c for c in a) + return pyopa.Sequence(a), pyopa.Sequence(b) + + +class EstimatePamTest(unittest.TestCase): + """estimate_pam (eigendecomposition of logPAM1) against the original EstimatePam""" + + def setUp(self): + self.defaults = pyopa.load_default_environments() + self.dms = pyopa.MutipleAlEnv(self.defaults['environments'], self.defaults['log_pam1']) + + def assertMatchesReference(self, dms, s1, s2): + res = dms.estimate_pam(s1, s2) + ref = dms._estimate_pam_reference(s1, s2) + msg = '%s / %s: %s != %s' % (s1.convert_readable(), s2.convert_readable(), res, ref) + # the score and the pam distance come from the same search over the DMS + self.assertEqual(res[:2], ref[:2], msg) + # the variance is computed with a different, but equivalent, formula + np.testing.assert_allclose(res[2], ref[2], rtol=1e-8, err_msg=msg) + + def test_random_alignments(self): + rng = random.Random(7) + for _ in range(300): + length = rng.choice([10, 50, 200, 1000]) + s1, s2 = random_alignment(rng, length, subst_rate=rng.random() ** 2, + gap_rate=rng.random() * 0.1) + self.assertMatchesReference(self.dms, s1, s2) + + def test_edge_cases(self): + cases = [ + ('ACDEFGHIKLMNPQRSTVWY', 'ACDEFGHIKLMNPQRSTVWY'), # identical, smallest pam + ('A', 'A'), + ('A', 'C'), + ('AAAAAAAAAA', 'CCCCCCCCCC'), # saturated, largest pam + ('ACDE__GHIK', 'AC__FGHIKL'), + ('ACDEXXGHIK', 'ACDEFGXXIK'), # X is not counted + ('ABZDEFGHIK', 'ABZDEFGHIK'), # ambiguous symbols + ('___', 'ACD'), + ('A_C', 'A_C'), + ] + for a, b in cases: + self.assertMatchesReference(self.dms, pyopa.Sequence(a), pyopa.Sequence(b)) + + def test_estimate_many_pam(self): + rng = random.Random(3) + ancestor = ''.join(rng.choice(AMINO_ACIDS) for _ in range(300)) + seqs = [pyopa.Sequence(ancestor)] + [pyopa.Sequence(''.join( + '_' if rng.random() < 0.02 else (rng.choice(AMINO_ACIDS) if rng.random() < 0.1 * i else c) + for c in ancestor)) for i in range(1, 8)] + score, dist, var = self.dms.estimate_many_pam(seqs, num_threads=2) + k = 0 + for i in range(len(seqs)): + for j in range(i + 1, len(seqs)): + self.assertEqual([score[k], dist[k], var[k]], self.dms.estimate_pam(seqs[i], seqs[j])) + k += 1 + + def test_non_reversible_log_pam1_falls_back(self): + # breaking detailed balance of logPAM1 (columns still sum to 0) + # makes estimate_pam use the original EstimatePam + q = self.defaults['log_pam1'].float64_matrix.copy() + a, b = (ord(c) - ord('A') for c in 'AC') + q[b, b] -= 0.5 * q[a, b] + q[a, b] *= 1.5 + # MutipleAlEnv only uses the float64_matrix of log_pam1 + dms = pyopa.MutipleAlEnv(self.defaults['environments'], types.SimpleNamespace(float64_matrix=q)) + + rng = random.Random(5) + seqs = [random_alignment(rng, 200, subst_rate=0.3, gap_rate=0.05) for _ in range(20)] + for s1, s2 in seqs: + self.assertEqual(dms.estimate_pam(s1, s2), dms._estimate_pam_reference(s1, s2)) + score, dist, var = dms.estimate_many_pam(list(seqs[0])) + self.assertEqual([score[0], dist[0], var[0]], dms._estimate_pam_reference(*seqs[0])) + + +if __name__ == '__main__': + unittest.main() From 1e4bb800d45b82915d263eae9115072075ba749a Mon Sep 17 00:00:00 2001 From: Adrian Altenhoff Date: Thu, 29 Oct 2026 21:53:01 +0100 Subject: [PATCH 2/2] Add the eigendecomposition of logPAM1 to the changelog --- CHANGELOG.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index df049af..42435e4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,13 @@ changes are released with the following version. - Requires Python 3.10 or newer; no more wheels for Python 3.6 to 3.9. Package metadata moved from `setup.py` to `pyproject.toml` ([#27](https://github.com/DessimozLab/pyopa/pull/27)). +- `MutipleAlEnv.estimate_pam` and `estimate_many_pam` are about 30 times + faster: they evaluate `exp(pam * logPAM1)` from an eigendecomposition + of logPAM1, only for the pairs of symbols in the alignment. Scores and + PAM distances are unchanged, variances agree to about 1e-9 relative. + If a custom logPAM1 is not reversible, the original implementation is + used. logPAM1 is now read only when the `MutipleAlEnv` is created + ([#30](https://github.com/DessimozLab/pyopa/pull/30)). - Wheels are built against numpy 2.x and work with numpy 1.x and 2.x ([#27](https://github.com/DessimozLab/pyopa/pull/27)). - The extension is built with OpenMP where the compiler supports it, and