include("../src/Util_multi.jl")
using MultipleParticle

function run_smc(N, lbs, ubs, f, ϵ₀, Δϵ_tol, α, Δϵ_window, output;
                 save_interval=0, target=0.234 , parallel=false, N_sample=100, γ=1., nu=1.)
    mkdir(output)
    lap_start = time()
    lap_end = 0.
    stepsize_flag = 0
    particles = generate_init_particles(N, lbs, ubs, f)

    ϵ_old = ϵ₀
    ϵ_new = 0.0
    cnt = 0
    terminate = false
    ESS=N
    Δϵ_mean = zeros(Δϵ_window)
    logfilename = joinpath(output, "log.csv")

    while terminate==false
        lap_end = time()
        lap = lap_end - lap_start
        lap_start = lap_end
        open(logfilename, "a") do fp
            writecsv(fp, Array([cnt, ϵ_old, ESS, lap])')
        end
        cnt += 1

        ϵ_new = calcPA(particles, α)
        Δϵ = ϵ_old - ϵ_new
        pop!(Δϵ_mean);
        Δϵ_mean = vcat(Δϵ, Δϵ_mean)
        if mean(Δϵ_mean) < Δϵ_tol && cnt > (Δϵ_window+100)
            terminate=true
        end

        # resampling
        candidates = deepcopy(particles)
        while ESS < N
            candidates = vcat(candidates, resample(particles, N-ESS))
            ESS = length(candidates)
        end

        F_old = get_obj_func_vals(candidates)
        candidates = reset_weights(candidates)
        accept_flags = zeros(length(candidates))
        if parallel == true
            candidates, accept_flags = MCstep_all_parallel(candidates, ϵ_old)
        else
            candidates, accept_flags = MCstep_all(candidates, ϵ_old)
        end


        particles = update_weights(candidates, F_old, ϵ_old, ϵ_new)
        particles = drop_zero_weights(particles)
        particles = normalize_weights(particles)


        acceptance = length(particles)/length(candidates)
        MH_acceptance = sum(accept_flags)/length(accept_flags)
        open(joinpath(output, "acceptance.csv"), "a") do fp
            writecsv(fp, Array([MH_acceptance acceptance]))
        end

        if MH_acceptance < target
            stepsize_flag = 1
        end
        if stepsize_flag == 0
            particles, γ, nu = update_stepsize_inRKHS_all(particles, γ, nu, target, MH_acceptance,
                                                          parallel=parallel, N=N_sample, adaptive=false)
        else
            particles, γ, nu = update_stepsize_inRKHS_all(particles, γ, nu, target, MH_acceptance,
                                                          parallel=parallel, N=N_sample, adaptive=true)
        end

        if length(particles) == 0
            println("particles were vanished at ϵ=$ϵ_new")
            break
        end

        if save_interval==0
            nothing
        elseif mod(cnt, save_interval) == 0 && terminate == false
            est_params = (lbs.*(ubs./lbs).^(hcat(get_values(particles)...)./10))'
            open(joinpath(output, @sprintf("est_param_%05d.csv", cnt)), "w") do fp
                writecsv(fp, est_params)
            end
        elseif terminate == true
            ESS = length(particles)
            lap_end = time()
            lap = lap_end - lap_start
            open(logfilename, "a") do fp
                writecsv(fp, Array([cnt, ϵ_old, Int(floor(ESS)), lap])')
            end

            est_params = (lbs.*(ubs./lbs).^(hcat(get_values(particles)...)./10))'
            open(joinpath(output, @sprintf("est_param_%05d.csv", cnt)), "w") do fp
                writecsv(fp, est_params)
            end
            break
        end

        ϵ_old = ϵ_new
        ESS = length(particles)
    end
end
