Skip to content

Substitution Models

Rate matrix eigendecomposition and transition probability computation.

EigenDecomp dataclass

Cached eigendecomposition of a reversible CTMC rate matrix.

Attributes:

Name Type Description
n_states int

Number of character states (4 for DNA).
freqs ndarray

float64, shape (n_states,). Stationary frequencies.
sqrt_pi ndarray

float64, shape (n_states,). Element-wise sqrt(freqs).
inv_sqrt_pi ndarray

float64, shape (n_states,). Element-wise 1 / sqrt(freqs).
eigvals ndarray

float64, shape (n_states,). Eigenvalues of the symmetrised matrix S, in ascending order.
U ndarray

float64, shape (n_states, n_states). Orthonormal eigenvectors of S (columns); used in the P(t) computation.
Source code in src/hifuku/substitution.py 
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
@dataclass
class EigenDecomp:
    """
    Cached eigendecomposition of a reversible CTMC rate matrix.

    Attributes
    ----------
    n_states : int
        Number of character states (4 for DNA).
    freqs : np.ndarray
        float64, shape (n_states,). Stationary frequencies.
    sqrt_pi : np.ndarray
        float64, shape (n_states,). Element-wise ``sqrt(freqs)``.
    inv_sqrt_pi : np.ndarray
        float64, shape (n_states,). Element-wise ``1 / sqrt(freqs)``.
    eigvals : np.ndarray
        float64, shape (n_states,). Eigenvalues of the symmetrised matrix S,
        in ascending order.
    U : np.ndarray
        float64, shape (n_states, n_states). Orthonormal eigenvectors of S
        (columns); used in the P(t) computation.
    """

    n_states: int           # 4 for DNA, 20 for protein
    freqs: np.ndarray       # shape (n_states,)
    sqrt_pi: np.ndarray     # shape (n_states,)
    inv_sqrt_pi: np.ndarray # shape (n_states,)
    eigvals: np.ndarray     # shape (n_states,)
    U: np.ndarray           # shape (n_states, n_states)

    def p_t(self, t: float) -> np.ndarray:
        """
        Return the transition probability matrix P(t) as float64 (n_states, n_states).

        P(t)[i,j] = probability of transitioning from state i to state j in time t.
        Rows sum to 1.  P(0) = identity.

        Notes
        -----
        - Felsenstein (1981) : transition matrix via matrix exponentiation.
        - Fourment et al. (2025) : .p_t pattern.
        """
        exp_lt = np.exp(self.eigvals * t)                   # (n_states,)
        M = (self.U * exp_lt) @ self.U.T                    # (n_states, n_states)
        # Tavaré (1986) : P = D_inv_sqrt M D_sqrt
        return self.inv_sqrt_pi[:, None] * M * self.sqrt_pi[None, :]

p_t 

p_t(t: float) -> np.ndarray

Return the transition probability matrix P(t) as float64 (n_states, n_states).

P(t)[i,j] = probability of transitioning from state i to state j in time t. Rows sum to 1. P(0) = identity.

Notes
Source code in src/hifuku/substitution.py 
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def p_t(self, t: float) -> np.ndarray:
    """
    Return the transition probability matrix P(t) as float64 (n_states, n_states).

    P(t)[i,j] = probability of transitioning from state i to state j in time t.
    Rows sum to 1.  P(0) = identity.

    Notes
    -----
    - Felsenstein (1981) : transition matrix via matrix exponentiation.
    - Fourment et al. (2025) : .p_t pattern.
    """
    exp_lt = np.exp(self.eigvals * t)                   # (n_states,)
    M = (self.U * exp_lt) @ self.U.T                    # (n_states, n_states)
    # Tavaré (1986) : P = D_inv_sqrt M D_sqrt
    return self.inv_sqrt_pi[:, None] * M * self.sqrt_pi[None, :]

jc69 

jc69() -> EigenDecomp

JC69 rate matrix eigendecomposition.

Equal stationary frequencies (1/4 each) and equal exchangeabilities. Normalised so that the mean substitution rate equals 1.

Returns:

Type Description
EigenDecomp

Ready to pass to log_likelihood or log_likelihood_gpu.
Notes
Source code in src/hifuku/substitution.py 
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
def jc69() -> EigenDecomp:
    """
    JC69 rate matrix eigendecomposition.

    Equal stationary frequencies (1/4 each) and equal exchangeabilities.
    Normalised so that the mean substitution rate equals 1.

    Returns
    -------
    EigenDecomp
        Ready to pass to ``log_likelihood`` or ``log_likelihood_gpu``.

    Notes
    -----
    - Felsenstein (1981) : JC69 as special case of GTR.
    """
    freqs = np.full(_N_DNA, 0.25)
    # Normalised JC69: Q_ii = -1, Q_ij = 1/3 for i ≠ j
    # Mean rate = -sum_i pi_i Q_ii = -4*(1/4)*(-1) = 1  ✓
    Q = np.full((_N_DNA, _N_DNA), 1.0 / 3.0)
    np.fill_diagonal(Q, -1.0)
    return _build_eigen(Q, freqs)

gtr 

gtr(rates: ndarray, freqs: ndarray) -> EigenDecomp

GTR rate matrix eigendecomposition.

Parameters:

Name Type Description Default
rates float64 (6,) - exchangeability parameters (AC, AG, AT, CG, CT, GT). 
These are symmetric: r_{ij} = r_{ji}.
required
freqs float64 (4,) - stationary base frequencies (A, C, G, T). 
Need not sum exactly to 1; normalized internally.
required
Notes
Source code in src/hifuku/substitution.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
def gtr(rates: np.ndarray, freqs: np.ndarray) -> EigenDecomp:
    """
    GTR rate matrix eigendecomposition.

    Parameters
    ----------
    rates : float64 (6,) - exchangeability parameters (AC, AG, AT, CG, CT, GT).
            These are symmetric: r_{ij} = r_{ji}.
    freqs : float64 (4,) - stationary base frequencies (A, C, G, T).
            Need not sum exactly to 1; normalized internally.

    Notes
    -----
    - Tavaré (1986) : GTR rate matrix parameterization.
    - Fourment et al. (2025) : eigendecomposition via symmetrisation.
    """
    rates = np.asarray(rates, dtype=np.float64)
    freqs = np.asarray(freqs, dtype=np.float64)

    if rates.shape != (6,):
        raise ValueError(f"rates must have shape (6,), got {rates.shape}")
    if freqs.shape != (4,):
        raise ValueError(f"freqs must have shape (4,), got {freqs.shape}")
    if np.any(rates <= 0):
        raise ValueError("All exchangeability rates must be positive")
    if np.any(freqs <= 0):
        raise ValueError("All stationary frequencies must be positive")

    freqs = freqs / freqs.sum()  # normalize

    Q = np.zeros((_N_DNA, _N_DNA), dtype=np.float64)
    for i in range(_N_DNA):
        for j in range(_N_DNA):
            if i != j:
                Q[i, j] = rates[_GTR_RATE_IDX[(i, j)]] * freqs[j]
    np.fill_diagonal(Q, -Q.sum(axis=1))

    # Normalise: mean substitution rate = 1
    # Fourment et al. (2025) : rate normalization
    mean_rate = -np.dot(freqs, np.diag(Q))
    Q /= mean_rate

    return _build_eigen(Q, freqs)