def makefractalCIJ(mx_lvl, E, sz_cl, seed=None):
    '''
    This function generates a directed network with a hierarchical modular
    organization. All modules are fully connected and connection density
    decays as 1/(E^n), with n = index of hierarchical level.

    Parameters
    ----------
    mx_lvl : int
        number of hierarchical levels, N = 2^mx_lvl
    E : int
        connection density fall off per level
    sz_cl : int
        size of clusters (must be power of 2)
    seed : hashable, optional
        If None (default), use the np.random's global random state to generate random numbers.
        Otherwise, use a new np.random.RandomState instance seeded with the given value.

    Returns
    -------
    CIJ : NxN np.ndarray
        connection matrix
    K : int
        number of connections present in output CIJ
    '''
    rng = get_rng(seed)
    # make a stupid little template
    t = np.ones((2, 2)) * 2

    # compute N and cluster size
    n = 2**mx_lvl
    sz_cl -= 1

    for lvl in range(1, mx_lvl):
        s = 2**(lvl + 1)
        CIJ = np.ones((s, s))
        grp1 = range(int(s / 2))
        grp2 = range(int(s / 2), s)
        ix1 = np.add.outer(np.array(grp1) * s, grp1).flatten()
        ix2 = np.add.outer(np.array(grp2) * s, grp2).flatten()
        CIJ.flat[ix1] = t  # numpy indexing is teh sucks :(
        CIJ.flat[ix2] = t
        CIJ += 1
        t = CIJ.copy()

    CIJ -= (np.ones((s, s)) + mx_lvl * np.eye(s))

    # assign connection probabilities
    ee = mx_lvl - CIJ - sz_cl
    ee = (ee > 0) * ee
    prob = (1 / E**ee) * (np.ones((s, s)) - np.eye(s))
    CIJ = (prob > rng.random_sample((n, n)))

    # count connections
    k = np.sum(CIJ)

    return np.array(CIJ, dtype=int), k