docs: def get_pij_matrix(t, diag, A, A_inv):
      """
      Calculates the probability matrix of substitutions i->j over time t,
      given the normalised generator diagonalisation.
  
  
      :param t: time
      :type t: float
      :return: probability matrix
      :rtype: numpy.ndarray
      """
      return A.dot(np.diag(np.exp(diag * t))).dot(A_inv)