r"""Adds artifical end states replacing source and sink sets.

    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
    -------
    F_new : (M+2, M+2) scipy.sparse matrix
        The artifical flux network with extra end states
    a_new : int
        The new single source a_new = M
    b_new : int
        The new single sink b_new = M+1

def add_endstates(F, A, B):
    r"""Adds artifical end states replacing source and sink sets.

    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
    -------
    F_new : (M+2, M+2) scipy.sparse matrix
        The artifical flux network with extra end states
    a_new : int
        The new single source a_new = M
    b_new : int
        The new single sink b_new = M+1

    """

    """Outgoing currents from A"""
    F = F.tocsr()
    outA = (F[A, :].sum(axis=1)).getA()[:, 0]

    """Incoming currents into B"""
    F = F.tocsc()
    inB = (F[:, B].sum(axis=0)).getA()[0, :]

    F = F.tocoo()
    M = F.shape[0]

    data_old = F.data
    row_old = F.row
    col_old = F.col

    """Add currents from new A=[n,] to all states in A"""
    row1 = np.zeros(outA.shape[0], dtype=np.int)
    row1[:] = M
    col1 = np.array(A)
    data1 = outA

    """Add currents from old B to new B=[n+1,]"""
    row2 = np.array(B)
    col2 = np.zeros(inB.shape[0], dtype=np.int)
    col2[:] = M + 1
    data2 = inB

    """Stack data, row and col arrays"""
    data = np.hstack((data_old, data1, data2))
    row = np.hstack((row_old, row1, row2))
    col = np.hstack((col_old, col1, col2))

    """New netflux matrix"""
    F_new = coo_matrix((data, (row, col)), shape=(M + 2, M + 2))

    return F_new, M, M + 1