def coarsegrain(F, sets):
    r"""Coarse-grains the flux to the given sets

    $fc_{i,j} = \sum_{i \in I,j \in J} f_{i,j}$
    Note that if you coarse-grain a net flux, it does not necessarily have a net
    flux property anymore. If want to make sure you get a netflux,
    use to_netflux(coarsegrain(F,sets)).

    Parameters
    ----------
    F : (n, n) ndarray
        Matrix of flux values between pairs of states.
    sets : list of array-like of ints
        The sets of states onto which the flux is coarse-grained.

    """
    nnew = len(sets)
    Fc = np.zeros((nnew, nnew))
    for i in range(0, nnew - 1):
        for j in range(i + 1, nnew):
            I = list(sets[i])
            J = list(sets[j])
            Fc[i, j] = np.sum(F[I, :][:, J])
            Fc[j, i] = np.sum(F[J, :][:, I])
    return Fc