code: def get_f81_pij(t, frequencies, mu):
      
  
      # 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
docs: 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