code: def flux(self, a, b):
          r
          # if b<=a:
          # raise ValueError("State index b has to be strictly larger than state index a")
          qminus = self.committor_backward(a, b)
          qplus = self.committor_forward(a, b)
          P = self.transition_matrix()
          pi = self.stationary_distribution()
  
          flux = pi[:, np.newaxis] * qminus[:, np.newaxis] * P * qplus[np.newaxis, :]
          ind = np.diag_indices(P.shape[0])
          flux[ind] = 0.0
          return flux
docs: r"""The flux network for the reaction from
          A=[0,...,a] => B=[b,...,M].
  
          Parameters
          ----------
          a : int
              State index
          b : int
              State index
  
          Returns
          -------
          flux : (M, M) ndarray
              Matrix of flux values between pairs of states.