r"""Decompose flux network into dominant reaction paths.

    Parameters
    ----------
    F : (M, M) scipy.sparse matrix
        The flux network (matrix of netflux values)
    A : array_like
        The set of starting states
    B : array_like
        The set of end states
    fraction : float, optional
        Fraction of total flux to assemble in pathway decomposition
    maxiter : int, optional
        Maximum number of pathways for decomposition
    tol : float, optional
        Floating point tolerance. The iteration is terminated once the
        relative capacity of all discovered path matches the desired
        fraction within floating point tolerance

    Returns
    -------
    paths : list
        List of dominant reaction pathways
    capacities: list
        List of capacities corresponding to each reactions pathway in paths

    References
    ----------
    .. [1] P. Metzner, C. Schuette and E. Vanden-Eijnden.
        Transition Path Theory for Markov Jump Processes.
        Multiscale Model Simul 7: 1192-1219 (2009)

def pathways(F, A, B, fraction=1.0, maxiter=1000, tol=1e-14):
    r"""Decompose flux network into dominant reaction paths.

    Parameters
    ----------
    F : (M, M) scipy.sparse matrix
        The flux network (matrix of netflux values)
    A : array_like
        The set of starting states
    B : array_like
        The set of end states
    fraction : float, optional
        Fraction of total flux to assemble in pathway decomposition
    maxiter : int, optional
        Maximum number of pathways for decomposition
    tol : float, optional
        Floating point tolerance. The iteration is terminated once the
        relative capacity of all discovered path matches the desired
        fraction within floating point tolerance

    Returns
    -------
    paths : list
        List of dominant reaction pathways
    capacities: list
        List of capacities corresponding to each reactions pathway in paths

    References
    ----------
    .. [1] P. Metzner, C. Schuette and E. Vanden-Eijnden.
        Transition Path Theory for Markov Jump Processes.
        Multiscale Model Simul 7: 1192-1219 (2009)

    """
    F, a, b = add_endstates(F, A, B)
    A = [a, ]
    B = [b, ]

    """Total flux"""
    TF = F.tocsr()[A, :].sum()

    """Total capacity fo all previously found reaction paths"""
    CF = 0.0
    niter = 0

    """List of dominant reaction pathways"""
    paths = []
    """List of corresponding capacities"""
    capacities = []

    while True:
        """Find dominant pathway of flux-network"""
        try:
            path = pathway(F, A, B)
        except PathwayError:
            break
        """Compute capacity of current pathway"""
        c = capacity(F, path)
        """Remove artifical end-states"""
        path = path[1:-1]
        """Append to lists"""
        paths.append(np.array(path))
        capacities.append(c)
        """Update capacity of all previously found paths"""
        CF += c
        """Remove capacity along given path from flux-network"""
        F = remove_path(F, path)
        niter += 1
        """Current flux numerically equals fraction * total flux or is
        greater equal than fraction * total flux"""
        if (abs(CF/TF - fraction) <= tol) or (CF/TF >= fraction):
            break
        if niter > maxiter:
            warnings.warn("Maximum number of iterations reached", RuntimeWarning)
            break
    return paths, capacities