export move_particle!, equilibrate_diagonal!

"""
    move_particle!(m::Model, e::Ensemble, c::Configuration) :: Tuple{Float64, Step}
Update that removes occupation in one kinkfree Orbital and adds it to another.
"""
function move_particle!(m::Model, e::Ensemble, c::Configuration) :: Tuple{Float64, Step}
    free_orbitals = filter(x -> isunaffected(c.kinks, x), c.occupations)
    if isempty(free_orbitals)
        return 1.0, Step()
    else
        x = rand( free_orbitals )
    end
    oe = filter(x -> isunaffected(c.kinks,x), setdiff!(sphere_with_same_spin(x, dk = ex_radius), c.occupations))

    #if there are no empty non interacting orbitals in neighbourhood make no change
    if isempty(oe)
        return 1.0, Step()
    end
    y = rand(oe)
    @assert x != y "same Configuration proposed."

    delta_di = ΔWdiag_element(m, e, c, y, x, ImgTime(0), ImgTime(1))
    @assert delta_di != Inf
    # weight change
    dw = exp(-(e.β*(energy(m,y)-energy(m,x)) + e.β*delta_di))

    # MC Step generated by this update

    Δ = Step((x,), (y,))
    apply_step!(c,Δ)

    # get orbitals for reverse update
    oe2 = filter(x -> isunaffected(c.kinks,x), setdiff!(sphere_with_same_spin(y, dk = ex_radius), c.occupations) )


    # quotient of proposal probabilities
    dv = length(oe)/length(oe2)
    @assert dv >= 0
    @assert dw >= 0 "dw=$(dw), x=$(x), y=$(y), delta_di=$(delta_di)"
    @assert !isinf(abs(dv*dw))
    return dv*dw, Δ
end

"""
    equilibrate_diagonal(m::Model, e::Ensemble, c::Configuration)
This function will propose 100 type-a-updates for each particle in c::Configuration.
As move_particle! can only be called on kinkfreee orbitals, calling this before starting the true equilibration
might reduce the error of the kinetic energy in runs with lots of kinks.
"""
function equilibrate_diagonal!(m::Model, e::Ensemble, c::Configuration)
    for _ in 1:e.N*1e2
        dv, Δ = move_particle!(m, e, c)
        if rand() > dv
            #deny step
            apply_back_step!(c, Δ)
        end
    end
end
