def flux_matrix(T, pi, qminus, qplus, netflux=True):
    r"""Compute the flux.

    Parameters
    ----------
    T : (M, M) scipy.sparse matrix
        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) scipy.sparse matrix
        Matrix of flux values between pairs of states.

    """
    D1 = diags((pi * qminus,), (0,))
    D2 = diags((qplus,), (0,))

    flux = D1.dot(T.dot(D2))

    """Remove self-fluxes"""
    flux = flux - diags(flux.diagonal(), 0)

    """Return net or gross flux"""
    if netflux:
        return to_netflux(flux)
    else:
        return flux