def reorder_mod(A, ci):
    '''
    This function reorders the connectivity matrix by modular structure and
    may hence be useful in visualization of modular structure.

    Parameters
    ----------
    A : NxN np.ndarray
        binary/weighted connectivity matrix
    ci : Nx1 np.ndarray
        module affiliation vector

    Returns
    -------
    On : Nx1 np.ndarray
        new node order
    Ar : NxN np.ndarray
        reordered connectivity matrix
    '''
    # TODO update function with 2015 changes

    from scipy import stats
    _, max_module_size = stats.mode(ci)
    u, ci = np.unique(ci, return_inverse=True)  # make consecutive
    n = np.size(ci)  # number of nodes
    m = np.size(u)  # number of modules

    nm = np.zeros((m,))  # number of nodes in modules
    knm = np.zeros((n, m))  # degree to other modules

    for i in range(m):
        nm[i] = np.size(np.where(ci == i))
        knm[:, i] = np.sum(A[:, ci == i], axis=1)

    am = np.zeros((m, m))  # relative intermodular connectivity
    for i in range(m):
        am[i, :] = np.sum(knm[ci == i, :], axis=0)
    am /= np.outer(nm, nm)

    # 1. Arrange densely connected modules together
    # symmetrized intermodular connectivity
    i, j = np.where(np.tril(am, -1) + 1)
    s = (np.tril(am, -1) + 1)[i, j]
    ord = np.argsort(s)[::-1]  # sort by high relative connectivity
    i = i[ord]
    j = j[ord]
    i += 1
    j += 1  # fix off by 1 error so np.where doesnt
    om = np.array((i[0], j[0]))  # catch module 0
    i[0] = 0
    j[0] = 0
    while len(om) < m:  # while not all modules ordered
        ui, = np.where(np.logical_and(
            i, np.logical_or(j == om[0], j == om[-1])))
        uj, = np.where(np.logical_and(
            j, np.logical_or(i == om[0], i == om[-1])))

        if np.size(ui):
            ui = ui[0]
        if np.size(uj):
            uj = uj[0]

        if ui == uj:
            i[ui] = 0
            j[uj] = 0
            continue

        if not np.size(ui):
            ui = np.inf
        if not np.size(uj):
            uj = np.inf
        if ui < uj:
            old = j[ui]
            new = i[ui]
        if uj < ui:
            old = i[uj]
            new = j[uj]
        if old == om[0]:
            om = np.append((new,), om)
        if old == om[-1]:
            om = np.append(om, (new,))

        i[i == old] = 0
        j[j == old] = 0

    print(om)

    # 2. Reorder nodes within modules
    on = np.zeros((n,), dtype=int)
    for y, x in enumerate(om):
        ind, = np.where(ci == x - 1)  # indices
        pos, = np.where(om == x)  # position
        # NOT DONE! OE NOES

        mod_imp = np.array((om, np.sign(np.arange(m) - pos),
                            np.abs(np.arange(m) - pos), am[x - 1, om - 1])).T
        print(np.shape((mod_imp[:, 3][::-1], mod_imp[:, 2])))
        ix = np.lexsort((mod_imp[:, 3][::-1], mod_imp[:, 2]))
        mod_imp = mod_imp[ix]
        # at this point mod_imp agrees with the matlab version
        signs = mod_imp[:, 1]
        mod_imp = np.abs(mod_imp[:, 0] * mod_imp[:, 1])
        mod_imp = np.append(mod_imp[1:], x)
        mod_imp = np.array(mod_imp - 1, dtype=int)
        print(mod_imp, signs)
        # at this point mod_imp is the absolute value of that in the matlab
        # version.  this limitation comes from sortrows ability to deal with
        # negative indices, which we would have to do manually.

        # instead, i punt on its importance; i only bother to order by the
        # principal dimension.  some within-module orderings
        # may potentially be a little bit out of order.

        # ksmi=knm[ind,:].T[mod_imp[::-1]]
        # reverse mod_imp to sort by the first column first and so on
        # print ksmi
        # for i,sin in enumerate(signs):
        #	if sin==-1:
        #		ksmi[i,:]=ksmi[i,:][::-1]
        # print ksmi
        # print np.shape(ksmi)

        # ^ this is unworkable and wrong, lexsort alone cannot handle the
        # negative indices problem of sortrows.  you would pretty much need
        # to rewrite sortrows to do lexsort plus negative indices; the algorithm
        # cant be further simplified.

        ord = np.lexsort(knm[np.ix_(ind, mod_imp[::-1])])
        # ord=np.lexsort(knm[ind,:].T[mod_imp[::-1]])
        if signs[mod_imp[0]] == -1:
            ord = ord[::-1]
            # reverse just the principal level and punt on the other levels.
            # this will basically be fine for most purposes and probably won't
            # ever show a difference for weighted graphs.
        on[ind[ord]] = y * int(max_module_size) + \
            np.arange(nm[x - 1], dtype=int)

    on = np.argsort(on)
    ar = A[np.ix_(on, on)]

    return on, ar