File: test_matrix_exponential_integration.py

package info (click to toggle)
python-cogent 2023.2.12a1%2Bdfsg-2%2Bdeb12u1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 12,416 kB
  • sloc: python: 89,165; makefile: 117; sh: 16
file content (150 lines) | stat: -rw-r--r-- 5,027 bytes parent folder | download
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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
from unittest import TestCase, main

import numpy as np

from numpy import array, diag, dot, exp

import cogent3.maths.matrix_exponentiation as cmme

from cogent3.maths import matrix_exponential_integration as expm


__author__ = "Ben Kaehler"
__copyright__ = "Copyright 2007-2014, The Cogent Project"
__credits__ = ["Ben Kaehler", "Ananias Iliadis", "Gavin Huttley"]
__license__ = "BSD-3"
__version__ = "2023.2.12a1"
__maintainer__ = "Ben Kaehler"
__email__ = "benjamin.kaehler@anu.edu.au"
__status__ = "Production"

from numpy.testing import assert_allclose


class TestIntegratingExponentiator(TestCase):
    def setUp(self) -> None:
        self.result = 0.7295333
        q = array([[0.5, 0.2, 0.1, 0.2]] * 4)
        for i in range(4):
            q[i, i] = 0.0
            q[i, i] = -sum(
                q[
                    i,
                ]
            )
        self.q = q
        self.p0 = array([0.2, 0.3, 0.3, 0.2])

    def test_van_loan_integrating_exponentiator(self):
        """VanLoanIntegratingExponentiator should reproduce Felsenstein
        analytic result, should throw if we pass it a defected matrix and ask
        it to use CheckedExponentiator, will work with a defective matrix (that
        we can integrate by hand) if we use the default RobustExponentiator,
        and should work for different choices of R and exponentiatior."""
        # Result from Von Bing's R code
        I = expm.VanLoanIntegratingExponentiator(self.q, -diag(self.q))(1.0)
        assert_allclose(dot(self.p0, I), self.result)

        self.assertRaises(
            ArithmeticError,
            expm.VanLoanIntegratingExponentiator,
            self.q,
            -diag(self.q),
            cmme.CheckedExponentiator,
        )

        Q = array([[1.0, 1.0], [0.0, 1.0]])

        def integral(t):
            return array(
                [[exp(t) - 1.0, exp(t) * (t - 1.0) + 1.0], [0.0, exp(t) - 1.0]]
            )

        assert_allclose(expm.VanLoanIntegratingExponentiator(Q)(1.0), integral(1.0))
        assert_allclose(expm.VanLoanIntegratingExponentiator(Q)(2.0), integral(2.0))

        R = array([[1.0], [1.0]])
        assert_allclose(
            expm.VanLoanIntegratingExponentiator(Q, R, cmme.TaylorExponentiator)(1.0),
            dot(integral(1.0), R),
        )

    def test_von_bing_integrating_exponentiator(self):
        """VonBingIntegratingExponentiator should reproduce Felsenstein
        analytic result, should throw if we pass it a defective matrix, and
        should match results obtained from VanLoanIntegratingExponentiator for
        a diagonisable matrix."""
        # Result from Von Bing's R code.
        I = expm.VonBingIntegratingExponentiator(self.q)(1.0)
        assert_allclose(dot(dot(self.p0, I), -diag(self.q)), self.result)

        self.assertRaises(
            ArithmeticError,
            expm.VonBingIntegratingExponentiator,
            array([[1.0, 1.0], [0.0, 1.0]]),
        )

        p = array(
            [
                [0.86758487, 0.05575623, 0.0196798, 0.0569791],
                [0.01827347, 0.93312148, 0.02109664, 0.02750842],
                [0.04782582, 0.1375742, 0.80046869, 0.01413129],
                [0.23022035, 0.22306947, 0.06995306, 0.47675713],
            ]
        )

        assert_allclose(
            expm.VonBingIntegratingExponentiator(p)(1.0),
            expm.VanLoanIntegratingExponentiator(
                p, exponentiator=cmme.FastExponentiator
            )(1.0),
        )
        assert_allclose(
            expm.VonBingIntegratingExponentiator(p)(2.0),
            expm.VanLoanIntegratingExponentiator(
                p, exponentiator=cmme.FastExponentiator
            )(2.0),
        )

    def test_repr(self):
        """repr() works for the integrating exponentiators"""
        for klass in (
            expm.VanLoanIntegratingExponentiator,
            expm.VonBingIntegratingExponentiator,
        ):
            i = klass(self.q)
            g = repr(i)
            self.assertIsInstance(g, str)

    def test_calc_number_subs(self):
        """correctly compute ENS"""
        mprobs = diag([0.1, 0.2, 0.3, 0.4])
        moprobs = array([0.1, 0.2, 0.3, 0.4])

        def get_calibrated_Q(R):
            Q = dot(R, mprobs)
            diag_add = diag(np.sum(Q, axis=1))
            to_divide = np.dot(moprobs, np.sum(Q, axis=1))
            Q -= diag_add
            Q /= to_divide
            return Q

        R = array([[0, 2, 1, 1], [2, 0, 1, 1], [1, 1, 0, 2], [1, 1, 2, 0]], dtype=float)

        Q = get_calibrated_Q(R)
        length = 0.1
        got = expm.expected_number_subs(moprobs, Q, length)
        assert_allclose(got, length)
        # case 2, length != ENS

        A = array(
            [[0, 1, 1, 1], [2, 0, 1, 1], [1, 1, 0, 40], [1, 1, 1, 0]], dtype=float
        )
        Q = get_calibrated_Q(A)
        length = 0.2
        got = expm.expected_number_subs(moprobs, Q, length)
        self.assertNotAlmostEqual(got, length)


if __name__ == "__main__":
    main()