def null_model_dir_sign(W, bin_swaps=5, wei_freq=.1, seed=None):
    '''
    This function randomizes an directed network with positive and
    negative weights, while preserving the degree and strength
    distributions. This function calls randmio_dir.m

    Parameters
    ----------
    W : NxN np.ndarray
        directed weighted connection matrix
    bin_swaps : int
        average number of swaps in each edge binary randomization. Default
        value is 5. 0 swaps implies no binary randomization.
    wei_freq : float
        frequency of weight sorting in weighted randomization. 0<=wei_freq<1.
        wei_freq == 1 implies that weights are sorted at each step.
        wei_freq == 0.1 implies that weights sorted each 10th step (faster,
            default value)
        wei_freq == 0 implies no sorting of weights (not recommended)
    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
    -------
    W0 : NxN np.ndarray
        randomized weighted connection matrix
    R : 4-tuple of floats
        Correlation coefficients between strength sequences of input and
        output connection matrices, rpos_in, rpos_out, rneg_in, rneg_out

    Notes
    -----
    The value of bin_swaps is ignored when binary topology is fully
       connected (e.g. when the network has no negative weights).
    Randomization may be better (and execution time will be slower) for
       higher values of bin_swaps and wei_freq. Higher values of bin_swaps may
       enable a more random binary organization, and higher values of wei_freq
       may enable a more accurate conservation of strength sequences.
    R are the correlation coefficients between positive and negative
       in-strength and out-strength sequences of input and output connection
       matrices and are used to evaluate the accuracy with which strengths
       were preserved. Note that correlation coefficients may be a rough
       measure of strength-sequence accuracy and one could implement more
       formal tests (such as the Kolmogorov-Smirnov test) if desired.
    '''
    rng = get_rng(seed)
    W = W.copy()
    n = len(W)
    np.fill_diagonal(W, 0)  # clear diagonal
    Ap = (W > 0)  # positive adjmat

    if np.size(np.where(Ap.flat)) < (n * (n - 1)):
        W_r = randmio_und_signed(W, bin_swaps, seed=rng)
        Ap_r = W_r > 0
        An_r = W_r < 0
    else:
        Ap_r = Ap
        An_r = An

    W0 = np.zeros((n, n))
    for s in (1, -1):
        if s == 1:
            Acur = Ap
            A_rcur = Ap_r
        else:
            Acur = An
            A_rcur = An_r

        Si = np.sum(W * Acur, axis=0)  # positive in-strength
        So = np.sum(W * Acur, axis=1)  # positive out-strength
        Wv = np.sort(W[Acur].flat)  # sorted weights vector
        i, j = np.where(A_rcur)
        Lij, = np.where(A_rcur.flat)  # weights indices

        P = np.outer(So, Si)

        if wei_freq == 0:  # get indices of Lij that sort P
            Oind = np.argsort(P.flat[Lij])  # assign corresponding sorted
            W0.flat[Lij[Oind]] = s * Wv  # weight at this index
        else:
            wsize = np.size(Wv)
            wei_period = np.round(1 / wei_freq)  # convert frequency to period
            lq = np.arange(wsize, 0, -wei_period, dtype=int)
            for m in lq:  # iteratively explore at this period
                # get indices of Lij that sort P
                Oind = np.argsort(P.flat[Lij])
                R = rng.permutation(m)[:np.min((m, wei_period))]
                for q, r in enumerate(R):
                    # choose random index of sorted expected weight
                    o = Oind[r]
                    W0.flat[Lij[o]] = s * Wv[r]  # assign corresponding weight

                    # readjust expected weighted probability for i[o],j[o]
                    f = 1 - Wv[r] / So[i[o]]
                    P[i[o], :] *= f
                    f = 1 - Wv[r] / So[j[o]]
                    P[j[o], :] *= f

                    # readjust in-strength of i[o]
                    So[i[o]] -= Wv[r]
                    # readjust out-strength of j[o]
                    Si[j[o]] -= Wv[r]

                O = Oind[R]
                # remove current indices from further consideration
                Lij = np.delete(Lij, O)
                i = np.delete(i, O)
                j = np.delete(j, O)
                Wv = np.delete(Wv, O)

    rpos_in = np.corrcoef(np.sum(W * (W > 0), axis=0),
                          np.sum(W0 * (W0 > 0), axis=0))
    rpos_ou = np.corrcoef(np.sum(W * (W > 0), axis=1),
                          np.sum(W0 * (W0 > 0), axis=1))
    rneg_in = np.corrcoef(np.sum(-W * (W < 0), axis=0),
                          np.sum(-W0 * (W0 < 0), axis=0))
    rneg_ou = np.corrcoef(np.sum(-W * (W < 0), axis=1),
                          np.sum(-W0 * (W0 < 0), axis=1))
    return W0, (rpos_in[0, 1], rpos_ou[0, 1], rneg_in[0, 1], rneg_ou[0, 1])