def align_matrices(m1, m2, dfun='sqrdiff', verbose=False, H=1e6, Texp=1,
                   T0=1e-3, Hbrk=10):
    '''
    This function aligns two matrices relative to one another by reordering
    the nodes in M2.  The function uses a version of simulated annealing.

    Parameters
    ----------
    M1 : NxN np.ndarray
        first connection matrix
    M2 : NxN np.ndarray
        second connection matrix
    dfun : str
        distance metric to use for matching
            'absdiff' : absolute difference
            'sqrdiff' : squared difference (default)
            'cosang' : cosine of vector angle
    verbose : bool
        print out cost at each iteration. Default False.
    H : int
        annealing parameter, default value 1e6
    Texp : int
        annealing parameter, default value 1. Coefficient of H s.t.
        Texp0=1-Texp/H
    T0 : float
        annealing parameter, default value 1e-3
    Hbrk : int
        annealing parameter, default value = 10. Coefficient of H s.t.
        Hbrk0 = H/Hkbr

    Returns
    -------
    Mreordered : NxN np.ndarray
        reordered connection matrix M2
    Mindices : Nx1 np.ndarray
        reordered indices
    cost : float
        objective function distance between M1 and Mreordered

    Notes
    -----
    Connection matrices can be weighted or binary, directed or undirected.
    They must have the same number of nodes.  M1 can be entered in any
    node ordering.

    Note that in general, the outcome will depend on the initial condition
    (the setting of the random number seed).  Also, there is no good way to
    determine optimal annealing parameters in advance - these parameters
    will need to be adjusted "by hand" (particularly H, Texp, T0, and Hbrk).
    For large and/or dense matrices, it is highly recommended to perform
    exploratory runs varying the settings of 'H' and 'Texp' and then select
    the best values.

    Based on extensive testing, it appears that T0 and Hbrk can remain
    unchanged in most cases.  Texp may be varied from 1-1/H to 1-10/H, for
    example.  H is the most important parameter - set to larger values as
    the problem size increases.  Good solutions can be obtained for
    matrices up to about 100 nodes.  It is advisable to run this function
    multiple times and select the solution(s) with the lowest 'cost'.

    If the two matrices are related it may be very helpful to pre-align them
    by reordering along their largest eigenvectors:
       [v,~] = eig(M1); v1 = abs(v(:,end)); [a1,b1] = sort(v1);
       [v,~] = eig(M2); v2 = abs(v(:,end)); [a2,b2] = sort(v2);
       [a,b,c] = overlapMAT2(M1(b1,b1),M2(b2,b2),'dfun',1);

    Setting 'Texp' to zero cancels annealing and uses a greedy algorithm
    instead.
    '''
    n = len(m1)
    if n < 2:
        raise BCTParamError("align_matrix will infinite loop on a singleton "
                            "or null matrix.")

    # define maxcost (greatest possible difference) and lowcost
    if dfun in ('absdiff', 'absdff'):
        maxcost = np.sum(np.abs(np.sort(m1.flat) - np.sort(m2.flat)[::-1]))
        lowcost = np.sum(np.abs(m1 - m2)) / maxcost
    elif dfun in ('sqrdiff', 'sqrdff'):
        maxcost = np.sum((np.sort(m1.flat) - np.sort(m2.flat)[::-1])**2)
        lowcost = np.sum((m1 - m2)**2) / maxcost
    elif dfun == 'cosang':
        maxcost = np.pi / 2
        lowcost = np.arccos(np.dot(m1.flat, m2.flat) /
                            np.sqrt(np.dot(m1.flat, m1.flat) * np.dot(m2.flat, m2.flat))) / maxcost
    else:
        raise BCTParamError('dfun must be absdiff or sqrdiff or cosang')

    mincost = lowcost
    anew = np.arange(n)
    amin = np.arange(n)
    h = 0
    hcnt = 0

    # adjust annealing parameters from user provided coefficients
    # H determines the maximal number of steps (user-provided)
    # Texp determines the steepness of the temperature gradient
    Texp = 1 - Texp / H
    # T0 sets the initial temperature and scales the energy term (user provided)
    # Hbrk sets a break point for the stimulation
    Hbrk = H / Hbrk

    while h < H:
        h += 1
        hcnt += 1
        # terminate if no new mincost has been found for some time
        if hcnt > Hbrk:
            break
        # current temperature
        T = T0 * (Texp**h)

        # choose two positions at random and flip them
        atmp = anew.copy()
        r1, r2 = rng.randint(n, size=(2,))
        while r1 == r2:
            r2 = rng.randint(n)
        atmp[r1] = anew[r2]
        atmp[r2] = anew[r1]
        m2atmp = m2[np.ix_(atmp, atmp)]
        if dfun in ('absdiff', 'absdff'):
            costnew = np.sum(np.abs(m1 - m2atmp)) / maxcost
        elif dfun in ('sqrdiff', 'sqrdff'):
            costnew = np.sum((m1 - m2atmp)**2) / maxcost
        elif dfun == 'cosang':
            costnew = np.arccos(np.dot(m1.flat, m2atmp.flat) / np.sqrt(
                np.dot(m1.flat, m1.flat) * np.dot(m2.flat, m2.flat))) / maxcost

        if costnew < lowcost or rng.random_sample() < np.exp(-(costnew - lowcost) / T):
            anew = atmp
            lowcost = costnew
            # is this the absolute best?
            if lowcost < mincost:
                amin = anew
                mincost = lowcost
                if verbose:
                    print('step %i ... current lowest cost = %f' % (h, mincost))
                hcnt = 0
            # if the cost is 0 we're done
            if mincost == 0:
                break
    if verbose:
        print('step %i ... final lowest cost = %f' % (h, mincost))

    M_reordered = m2[np.ix_(amin, amin)]
    M_indices = amin
    cost = mincost
    return M_reordered, M_indices, cost