r"""Compute the dominant reaction-pathway.

    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

    Returns
    -------
    w : list
        The dominant reaction-pathway

def pathway(F, A, B):
    r"""Compute the dominant reaction-pathway.

    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

    Returns
    -------
    w : list
        The dominant reaction-pathway

    """
    if F.nnz == 0:
        raise PathwayError('no more pathways left: Flux matrix does not contain any positive entries')
    b1, b2, F = find_bottleneck(F, A, B)
    if np.any(A == b1):
        wL = [b1, ]
    elif np.any(B == b1):
        raise PathwayError(("Roles of vertices b1 and b2 are switched."
                            "This should never happen for a correct flux network"
                            "obtained from a reversible transition meatrix."))
    else:
        wL = pathway(F, A, [b1, ])
    if np.any(B == b2):
        wR = [b2, ]
    elif np.any(A == b2):
        raise PathwayError(("Roles of vertices b1 and b2 are switched."
                            "This should never happen for a correct flux network"
                            "obtained from a reversible transition meatrix."))
    else:
        wR = pathway(F, [b2, ], B)
    return wL + wR