include("lattices.jl")

mutable struct ElectronTinder
    """ Database to match lonely up spins with lonely down spins by site"""
    lonely_up::Array{Int}
    lonely_down::Array{Int}
end

function addUser!(t::ElectronTinder,site::Int,spin::Int)
    spin==1 ? append!(t.lonely_up,site) : append!(t.lonely_down,site)
end


function SwipeRight(database::ElectronTinder)
    """ Creates list of matches (moves) generated by
    electron tinder's database """
    @assert length(database.lonely_up) == length(database.lonely_down)
    lu = database.lonely_up
    ld = database.lonely_down
    nmatches = length(lu)
    matches = Array{ExchangeMove}(undef,nmatches)
    for i in 1:nmatches
        matches[i] = ExchangeMove((lu[i],ld[i]))
    end
    return matches
end

using Statistics: mean
using LightGraphs: neighborhood
function Swap(chain::DoubleGutzwillerChain,regionsize::Int)
    """ regionsize ≈ r, the radius of the region selected """

    states = get_State(chain) # (state1, state2)
    lg = get_Model(chain)["lattice"].graph
    ns = prod(get_Model(chain)["dims"])

    type = typeof(states[1].det_A_up)
    swap_per_site = zeros(type,ns)

    # redundancy for mixed BC lattice
    lg0 = SimpleGraph(adjacency_matrix(lg))

    for v = 1:ns
        sites = unique(neighborhood(lg0,v,regionsize))
        swap_per_site[v] = compute_swapregion(chain,sites)
    end

    return mean(swap_per_site)
end

function compute_swapregion(chain::DoubleGutzwillerChain,
                            sites::Array)
    """ This function is somewhat nontrivial, the following is the
    layout of how it works. First, a list of exchange moves is created
    for each replica that maps the current configuration to the new
    one after swapping. This happens in a two stages: first the spin
    spin configurations inside the regions are compared to see where they
    disagree. Every site where they disagree, an entry is added the
    appropriate spin list in Electron Tinder. Finally, the completed
    ElectronTinder database is fed into SwipeRight to generate a list of moves.

    After this, the list of moves are submitted to preexisting functions
    to compute the determinant ratio produced by that sequence of moves.

    Finally what is returned is the expectation value of the swap operator
    acting on this region which is,
    det(swapped copy 1)*det(swapped copy 2) / (det(og copy 1) * det(og copy 2))

    Later, I may modify this to be more general/support more than 2 copies,
    but for now it seems like thats more trouble than its worth.
    """
    s0 = get_State(chain)
    (spinc1, spinc2) = (s0[1].spin_config.sc,
                        s0[2].spin_config.sc)

    if sum(spinc1[sites])!=sum(spinc2[sites]) # regions must have same quantum #
        return 0
    end
    if spinc1[sites] == spinc2[sites]
        return 1
    end


    Tinder = (ElectronTinder([],[]),ElectronTinder([],[]))
    for site in sites
        if spinc1[site] != spinc2[site]
            addUser!(Tinder[1],site,spinc1[site])
            addUser!(Tinder[2],site,spinc2[site])
        end
    end

    (MoveList1, MoveList2) = SwipeRight.(Tinder)
    @assert size(MoveList1) == size(MoveList2)



    # now because update_chain! modifies the state, make copies of
    # the two states to update
    free_chains = deepcopy(get_Replicas(chain))
    free_states = get_State.(free_chains)

    @assert length(MoveList1) == length(MoveList2) "compute swapregion: move lists
    for two copies are of different length"

    # try to do "fast" determinant update with all moves simultaneously
    n_moves = length(MoveList1)
    num_particles = length(free_states[1].r_up) # not very pretty

    type = typeof(s0[1].det_A_up) # complexf64 or f64

    U_up_copy1 = zeros(type,(num_particles,n_moves))
    U_down_copy1 = zeros(type,(num_particles,n_moves))
    V_up_copy1 = zeros(type,(num_particles,n_moves))
    V_down_copy1 = zeros(type,(num_particles,n_moves))

    U_up_copy2 = zeros(type,(num_particles,n_moves))
    U_down_copy2 = zeros(type,(num_particles,n_moves))
    V_up_copy2 = zeros(type,(num_particles,n_moves))
    V_down_copy2 = zeros(type,(num_particles,n_moves))

    for i in 1:n_moves
        # update_Rs!(free_states[1],MoveList1[i])
        # update_Rs!(free_states[2],MoveList2[i])

        U_up_copy1[:,i], V_up_copy1[:,i], U_down_copy1[:,i], V_down_copy1[:,i] =
        get_update_vectors(free_states[1],MoveList1[i])

        U_up_copy2[:,i], V_up_copy2[:,i], U_down_copy2[:,i], V_down_copy2[:,i] =
        get_update_vectors(free_states[2],MoveList2[i])

    end

    # its important to do the determinant ratio updates all at once
    # brute force reevaluation of determinants code left in for debugging

    #
    # filled_states = get_Wavefunctions(free_chains[1])
    #
    # det1up = get_Determinant(free_states[1].r_up,filled_states)
    # det1down = get_Determinant(free_states[1].r_down,filled_states)
    # det2up = get_Determinant(free_states[2].r_up,filled_states)
    # det2down = get_Determinant(free_states[2].r_down,filled_states)


    det_ratio_1u = det_ratio_factor(free_states[1].A_inv_up,U_up_copy1,V_up_copy1)
    det_ratio_1d = det_ratio_factor(free_states[1].A_inv_down,U_down_copy1,V_down_copy1)
    det_ratio_2u = det_ratio_factor(free_states[2].A_inv_up,U_up_copy2,V_up_copy2)
    det_ratio_2d = det_ratio_factor(free_states[2].A_inv_down,U_down_copy2,V_down_copy2)

    fast_ratio = det_ratio_1u * det_ratio_1d * det_ratio_2u * det_ratio_2d

    # ss1 = get_test_state(free_chains[1],free_states[1].r_up, free_states[1].r_down)
    # ss2 = get_test_state(free_chains[2],free_states[2].r_up, free_states[2].r_down)

    # sdp = det1up * det1down * det2up * det2down
    # odp = s0[1].det_A_up * s0[1].det_A_down * s0[2].det_A_up * s0[2].det_A_down


    # print("sdp: ",sdp," odp: ",odp,"\n")


    return fast_ratio
end
