def simulated_annealing(coordinates, route=None, epochs=1, T=20.0, eta=0.995):
    """Simulated Annealing for TSP

    T(n) = T * eta**n

    Args:
        coordinates (np.ndarray): initial route consisting of coordinates
        route (array): inital route
        epochs (int): number of loops through all cities.
        T (float): Initial temperature
        eta (float): Exponential cooling rate

    Returns:
        route (np.ndarray): optimized route
        distances (np.ndarray): travel distance after each epoch
    """
    np.random.seed(7)

    N = coordinates.shape[0]

    if route is None:
        route = np.random.permutation(N)

    distances = []
    distances.append(distance(coordinates, route))

    T_n = [T * eta**i for i in range(epochs)]

    for e in range(epochs):
        if e==0:
            prev_distance = distance(coordinates, route)

        for i, idx in enumerate(np.random.permutation(N)):
            route[idx], route[(idx+1) % N] = route[(idx+1) % N], route[idx]
            curr_distance = distance(coordinates, route)

            dE = curr_distance - prev_distance
            if dE > 0:
                p = np.exp(-dE/T_n[e])

                if np.random.binomial(1,p,1) == 0: # undo changes
                    route[idx], route[(idx+1) % N] = route[(idx+1) % N], route[idx]
                    # prev_distance remains unchaged
                else: # accept
                    prev_distance = curr_distance
            else: # accept
                prev_distance = curr_distance


        distances.append(distance(coordinates, route))

    print(distance)
    return np.asarray(route), np.asarray(distances)