-
Notifications
You must be signed in to change notification settings - Fork 6
/
Copy pathPijTest.py
executable file
·60 lines (51 loc) · 2.68 KB
/
PijTest.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import unittest
import numpy as np
from pastml.models.CustomRatesModel import CustomRatesModel
from pastml.models.F81Model import F81Model
from pastml.models.HKYModel import HKYModel
from pastml.models.JCModel import JCModel
from pastml.models.HKYModel import A, G, T, C, HKY_STATES
class PijTest(unittest.TestCase):
def test_pij_f81(self):
for _ in range(10):
t = 10 * np.random.rand(1)[0]
freqs = np.random.rand(10)
freqs /= freqs.sum()
model = CustomRatesModel(sf=1, states=np.array(list('ABCDEFGHIJ')), forest_stats=None, frequencies=freqs,
rate_matrix=np.ones(shape=(10, 10), dtype=np.float64) - np.eye(10))
p_ij = model.get_Pij_t(t)
model = F81Model(sf=1, forest_stats=None, frequencies=freqs, states=np.array(list('ABCDEFGHIJ')))
p_ij_f81 = model.get_Pij_t(t)
self.assertTrue(np.allclose(p_ij, p_ij_f81),
msg='F81 P_ij calculation failed for time {} and frequencies {}'.format(t, freqs))
def test_pij_jc(self):
freqs = np.ones(10)
freqs /= freqs.sum()
for _ in range(10):
t = 10 * np.random.rand(1)[0]
model = CustomRatesModel(sf=1, states=np.array(list('ABCDEFGHIJ')), forest_stats=None, frequencies=freqs,
rate_matrix=np.ones(shape=(10, 10), dtype=np.float64) - np.eye(10))
p_ij = model.get_Pij_t(t)
model = JCModel(sf=1, forest_stats=None, states=np.array(list('ABCDEFGHIJ')))
p_ij_jc = model.get_Pij_t(t)
self.assertTrue(np.allclose(p_ij, p_ij_jc),
msg='JC P_ij calculation failed for time {} and frequencies {}'.format(t, freqs))
def test_pij_hky(self):
n = 4
for _ in range(10):
t = 10 * np.random.rand(1)[0]
kappa = 20 * np.random.rand(1)[0]
freqs = np.random.rand(n)
freqs /= freqs.sum()
rate_matrix = np.ones(shape=(n, n), dtype=np.float64) - np.eye(n)
rate_matrix[A, G] = kappa
rate_matrix[C, T] = kappa
rate_matrix[G, A] = kappa
rate_matrix[T, C] = kappa
model = CustomRatesModel(sf=1, states=HKY_STATES, forest_stats=None, frequencies=freqs,
rate_matrix=rate_matrix)
p_ij = model.get_Pij_t(t)
model = HKYModel(sf=1, forest_stats=None, kappa=kappa, frequencies=freqs)
p_ij_hky = model.get_Pij_t(t)
self.assertTrue(np.allclose(p_ij, p_ij_hky),
msg='HKY P_ij calculation failed for time {} and frequencies {}'.format(t, freqs))