def flux_producers(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.

    """
    n = F.shape[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 out flux absolute
    prod_abs = np.maximum(outfluxes - influxes, np.zeros(n))
    # net out flux relative
    prod_rel = prod_abs / (np.maximum(outfluxes, influxes))
    # return all indexes that are produces in terms of absolute and relative tolerance
    return list(np.where((prod_abs > atol) * (prod_rel > rtol))[0])