docs: def flux_consumers(F, rtol=1e-05, atol=1e-12):
      r"""Return indexes of states that are net flux producers.
  
      Parameters
      ----------
      F : (n, n) ndarray
          Matrix of flux values between pairs of states.
      rtol : float
          relative tolerance. fulfilled if max(outflux-influx, 0) / max(outflux,influx) < rtol
      atol : float
          absolute tolerance. fulfilled if max(outflux-influx, 0) < atol
  
      Returns
      -------
      producers : (n) ndarray of int
          indexes of states that are net flux producers. May include "dirty" producers, i.e.
          states that have influx but still produce more outflux and thereby violate flux
          conservation.
  
      """
      # can be used with sparse or dense
      n = np.shape(F)[0]
      influxes = np.array(np.sum(F, axis=0)).flatten()  # all that flows in
      outfluxes = np.array(np.sum(F, axis=1)).flatten()  # all that flows out
      # net in flux absolute
      con_abs = np.maximum(influxes - outfluxes, np.zeros(n))
      # net in flux relative
      con_rel = con_abs / (np.maximum(outfluxes, influxes))
      # return all indexes that are produces in terms of absolute and relative tolerance
      return list(np.where((con_abs > atol) * (con_rel > rtol))[0])