docs: def flux_matrix(T, pi, qminus, qplus, netflux=True):
      r"""Compute the TPT flux network for the reaction A-->B.
  
      Parameters
      ----------
      T : (M, M) ndarray
          transition matrix
      pi : (M,) ndarray
          Stationary distribution corresponding to T
      qminus : (M,) ndarray
          Backward comittor
      qplus : (M,) ndarray
          Forward committor
      netflux : boolean
          True: net flux matrix will be computed
          False: gross flux matrix will be computed
  
      Returns
      -------
      flux : (M, M) ndarray
          Matrix of flux values between pairs of states.
  
      Notes
      -----
      Computation of the flux network relies on transition path theory
      (TPT). The central object used in transition path theory is the
      forward and backward comittor function.
  
      See also
      --------
      committor.forward_committor, committor.backward_committor
  
  
      """
      ind = np.diag_indices(T.shape[0])
      flux = pi[:, np.newaxis] * qminus[:, np.newaxis] * T * qplus[np.newaxis, :]
      """Remove self fluxes f_ii"""
      flux[ind] = 0.0
      """Return net or gross flux"""
      if netflux:
          return to_netflux(flux)
      else:
          return flux