def findpaths(CIJ, qmax, sources, savepths=False):
    '''
    Paths are sequences of linked nodes, that never visit a single node
    more than once. This function finds all paths that start at a set of
    source nodes, up to a specified length. Warning: very memory-intensive.

    Parameters
    ----------
    CIJ : NxN np.ndarray
        binary directed/undirected connection matrix
    qmax : int
        maximal path length
    sources : Nx1 np.ndarray
        source units from which paths are grown
    savepths : bool
        True if all paths are to be collected and returned. This functionality
        is currently not enabled.

    Returns
    -------
    Pq : NxNxQ np.ndarray
        Path matrix with P[i,j,jq] = number of paths from i to j with length q
    tpath : int
        total number of paths found
    plq : Qx1 np.ndarray
        path length distribution as a function of q
    qstop : int
        path length at which findpaths is stopped
    allpths : None
        a matrix containing all paths up to qmax. This function is extremely
        complicated and reimplementing it in bctpy is not straightforward.
    util : NxQ np.ndarray
        node use index

    Notes
    -----
    Note that Pq(:,:,N) can only carry entries on the diagonal, as all
    "legal" paths of length N-1 must terminate.  Cycles of length N are
    possible, with all vertices visited exactly once (except for source and
    target). 'qmax = N' can wreak havoc (due to memory problems).

    Note: Weights are discarded.
    Note: I am certain that this algorithm is rather inefficient -
    suggestions for improvements are welcome.

    '''
    CIJ = binarize(CIJ, copy=True)  # ensure CIJ is binary
    n = len(CIJ)
    k = np.sum(CIJ)
    pths = []
    Pq = np.zeros((n, n, qmax))
    util = np.zeros((n, qmax))

    # this code is for pathlength=1
    # paths are seeded from sources
    q = 1
    for j in range(n):
        for i in range(len(sources)):
            i_s = sources[i]
            if CIJ[i_s, j] == 1:
                pths.append([i_s, j])
    pths = np.array(pths)

    # calculate the use index per vertex (for paths of length 1)
    util[:, q], _ = np.histogram(pths, bins=n)
    # now enter the found paths of length 1 into the pathmatrix Pq
    for nrp in range(np.size(pths, axis=0)):
        Pq[pths[nrp, 0], pths[nrp, q], q - 1] += 1

    # begin saving allpths
    if savepths:
        allpths = pths.copy()
    else:
        allpths = []

    npthscnt = k

    # big loop for all other pathlengths q
    for q in range(2, qmax + 1):
        # to keep track of time...
        print((
            'current pathlength (q=i, number of paths so far (up to q-1)=i' % (q, np.sum(Pq))))

        # old paths are now in 'pths'
        # new paths are about to be collected in 'npths'
        # estimate needed allocation for new paths
        len_npths = np.min((np.ceil(1.1 * npthscnt * k / n), 100000000))
        npths = np.zeros((q + 1, len_npths))

        # find the unique set of endpoints of 'pths'
        endp = np.unique(pths[:, q - 1])
        npthscnt = 0

        for i in endp:  # set of endpoints of previous paths
            # in 'pb' collect all previous paths with 'i' as their endpoint
            pb, = np.where(pths[:, q - 1] == i)
            # find the outgoing connections from i (breadth-first)
            nendp, = np.where(CIJ[i, :] == 1)
            # if i is not a dead end
            if nendp.size:
                for j in nendp:  # endpoints of next edge
                    # find new paths -- only legal ones, no vertex twice
                    # visited
                    pb_temp = pb[np.sum(j == pths[pb, 1:q], axis=1) == 0]

                    # add new paths to 'npths'
                    pbx = pths[pb_temp - 1, :]
                    npx = np.ones((len(pb_temp), 1)) * j
                    npths[:, npthscnt:npthscnt + len(pb_temp)] = np.append(
                        pbx, npx, axis=1).T
                    npthscnt += len(pb_temp)
                    # count new paths and add the number to P
                    Pq[:n, j, q -
                        1] += np.histogram(pths[pb_temp - 1, 0], bins=n)[0]

        # note: 'npths' now contains a list of all the paths of length q
        if len_npths > npthscnt:
            npths = npths[:, :npthscnt]

        # append the matrix of all paths
        # FIXME
        if savepths:
            raise NotImplementedError("Sorry allpaths is not yet implemented")

        # calculate the use index per vertex (correct for cycles, count
        # source/target only once)
        util[:, q - 1] += (np.histogram(npths[:, :npthscnt], bins=n)[0] -
                           np.diag(Pq[:, :, q - 1]))

        # elininate cycles from "making it" to the next level, so that "pths"
        # contains all the paths that have a chance of being continued
        if npths.size:
            pths = np.squeeze(npths[:, np.where(npths[0, :] != npths[q, :])]).T
        else:
            pths = []

        # if there are no 'pths' paths left, end the search
        if not pths.size:
            qstop = q
            tpath = np.sum(Pq)
            plq = np.sum(np.sum(Pq, axis=0), axis=0)
            return

    qstop = q
    tpath = np.sum(Pq)  # total number of paths
    plq = np.sum(np.sum(Pq, axis=0), axis=0)  # path length distribution

    return Pq, tpath, plq, qstop, allpths, util
Paths are sequences of linked nodes, that never visit a single node
    more than once. This function finds all paths that start at a set of
    source nodes, up to a specified length. Warning: very memory-intensive.

    Parameters
    ----------
    CIJ : NxN np.ndarray
        binary directed/undirected connection matrix
    qmax : int
        maximal path length
    sources : Nx1 np.ndarray
        source units from which paths are grown
    savepths : bool
        True if all paths are to be collected and returned. This functionality
        is currently not enabled.

    Returns
    -------
    Pq : NxNxQ np.ndarray
        Path matrix with P[i,j,jq] = number of paths from i to j with length q
    tpath : int
        total number of paths found
    plq : Qx1 np.ndarray
        path length distribution as a function of q
    qstop : int
        path length at which findpaths is stopped
    allpths : None
        a matrix containing all paths up to qmax. This function is extremely
        complicated and reimplementing it in bctpy is not straightforward.
    util : NxQ np.ndarray
        node use index

    Notes
    -----
    Note that Pq(:,:,N) can only carry entries on the diagonal, as all
    "legal" paths of length N-1 must terminate.  Cycles of length N are
    possible, with all vertices visited exactly once (except for source and
    target). 'qmax = N' can wreak havoc (due to memory problems).

    Note: Weights are discarded.
    Note: I am certain that this algorithm is rather inefficient -
    suggestions for improvements are welcome.