Calculate the probability of substitution i->j over time t, given the mutation rate mu:
    For F81 (and JC which is a simpler version of it)
    Pij(t) = \pi_j (1 - exp(-mu t)) + exp(-mu t), if i == j, \pi_j (1 - exp(-mu t)), otherwise
    [Gascuel "Mathematics of Evolution and Phylogeny" 2005].

    :param frequencies: array of state frequencies \pi_i
    :type frequencies: numpy.array
    :param mu: mutation rate: \mu = 1 / (1 - sum_i \pi_i^2)
    :type mu: float
    :param t: time t
    :type t: float
    :param sf: scaling factor by which t should be multiplied.
    :type sf: float
    :return: probability matrix
    :rtype: numpy.ndarray

def get_f81_pij(t, frequencies, mu):
    """
    Calculate the probability of substitution i->j over time t, given the mutation rate mu:
    For F81 (and JC which is a simpler version of it)
    Pij(t) = \pi_j (1 - exp(-mu t)) + exp(-mu t), if i == j, \pi_j (1 - exp(-mu t)), otherwise
    [Gascuel "Mathematics of Evolution and Phylogeny" 2005].

    :param frequencies: array of state frequencies \pi_i
    :type frequencies: numpy.array
    :param mu: mutation rate: \mu = 1 / (1 - sum_i \pi_i^2)
    :type mu: float
    :param t: time t
    :type t: float
    :param sf: scaling factor by which t should be multiplied.
    :type sf: float
    :return: probability matrix
    :rtype: numpy.ndarray
    """

    # if mu == inf (e.g. just one state) and t == 0, we should prioritise mu
    exp_mu_t = 0. if (mu == np.inf) else np.exp(-mu * t)
    return (1 - exp_mu_t) * frequencies + np.eye(len(frequencies)) * exp_mu_t