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 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()