module Demo
    using Dates
    using Plots
    #using LaTeXStrings

    import ..System_parameters:Params_set,parameterloading
    import ..LTK_universe:Universe
    import  ..MD:construct_MD_parameters
    import ..Measurements:Measurement_set,measurements
    import ..Verbose_print:println_verbose1,println_verbose2,Verbose_1
    import ..Heatbath:heatbath!


    system = Dict()
    actions = Dict()
    md = Dict()
    cg = Dict()
    wilson = Dict()
    staggered = Dict()
    measurement = Dict()
    # - - parameters - - - - - - - - - - - 
# - - parameters - - - - - - - - - - - 
    system["saveU_dir"] = ""
    system["verboselevel"] = 1
    #system["L"] = (6, 6, 6, 6)
    system["L"] = (4, 4, 4, 4)
    system["Nwing"] = 1
    system["Nsteps"] = 10000
    system["quench"] = true
    system["logfile"] = "Heatbath_L06060606_beta8.0_quenched.txt"
    system["initial"] = "hot"
    system["Dirac_operator"] = nothing
    system["log_dir"] = "./logs"
    system["Nthermalization"] = 10
    system["update_method"] = "Heatbath"
    system["randomseed"] = 111
    system["NC"] = 3
    system["saveU_every"] = 1
    system["BoundaryCondition"] = [1, 1, 1, -1]
    system["saveU_format"] = nothing
    #system["β"] = 6.93015
    system["β"] = 6.75850661032668353
    actions["use_autogeneratedstaples"] = false
    actions["couplingcoeff"] = Any[]
    actions["couplinglist"] = Any[]
    md["Δτ"] = 0.05
    md["SextonWeingargten"] = false
    md["MDsteps"] = 20
    cg["eps"] = 1.0e-19
    cg["MaxCGstep"] = 3000
    wilson["Clover_coefficient"] = 0
    wilson["r"] = 1
    wilson["hop"] = 0
    staggered["Nf"] = 0
    staggered["mass"] = 0
    measurement["measurement_methods"] = Dict[Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Polyakov_loop"), Dict{Any,Any}("fermiontype" => nothing,"measure_every" => 1,"methodname" => "Plaquette")]
    measurement["measurement_dir"] = "Heatbath_L06060606_beta8.0_quenched"
    measurement["measurement_basedir"] = "./measurements"
    # - - - - - - - - - - - - - - - - - - -

    function demo()

        params_set = Params_set(system,actions,md,cg,wilson,staggered,measurement)

        parameters = parameterloading(params_set)
        univ = Universe(parameters)

        mdparams = construct_MD_parameters(parameters)

        meas = Measurement_set(univ,parameters.measuredir,measurement_methods=parameters.measurement_methods)
        run_demo!(parameters,univ,meas)

    end

    function run_demo!(parameters,univ,meas)
        plt1 = histogram([0],label=nothing) #plot1
        plt2 = plot([], [],label=nothing) #plot2
        ylabel!("Plaquette")
        xlabel!("MC time")
        plt3 = histogram([0],label=nothing) #plot3
        plt4 = plot([], [],label=nothing) #plot4
        ylabel!("|Polyakov loop|")
        xlabel!("MC time")
        plt5 = scatter([], [],label=nothing,title="Polyakov loop") #plot5 scatter
        ylabel!("Im")
        xlabel!("Re")
        plt6 = plot([], [],label=nothing) #plot6
        ylabel!("Arg(Polyakov loop)")
        xlabel!("MC time")
        plot(plt1, plt2, plt3, plt4, plt5, plt6, layout = 6)

        hist_plaq=[]
        hist_poly=[]
        hist_poly_θ=[]

        function plot_refresh!(plt1,plt2,plt3,plt4,plt5,plt6,hist_plaq,hist_poly,hist_poly_θ, plaq, poly, itrj)
            bins=round(Int,log(itrj)*4+1)
            if itrj < 500
                bins=round(Int,log(itrj)*1.5+1)
            elseif itrj < 1000
                bins=round(Int,log(itrj)*2.5+1)
            elseif itrj < 2000
                bins=round(Int,log(itrj)*3.5+1)
            end
            append!(hist_plaq, plaq)
            append!(hist_poly, abs(poly))
            # omit un-thermalized part
            if (10 < itrj < 1000) & (itrj%5==0) 
                print("remove unthermalize part $(length(hist_plaq)) -> ")
                popfirst!(hist_plaq)
                popfirst!(hist_poly)
                println(" $(length(hist_plaq))")
            end
            #
            plt1 = histogram(hist_plaq,bins=bins,label=nothing) #plot1
            #plot!(plt1,title="SU(3), Quenched, L=6^4,")
            plot!(plt1,title="SU(3), Quenched, L=4^4,")
            xlabel!("Plaquette")
            plot!(plt2,title="Heatbath")
            plt3 = histogram(hist_poly,bins=bins,label=nothing) #plot3
            xlabel!("|Polyakov loop|")
            #
            if itrj<500
                push!(plt2, 1, itrj, plaq)
            elseif 500<itrj<1000
                if itrj%2==0
                push!(plt2, 1, itrj, plaq)
                end
            elseif 1000<itrj
                if itrj%5==0
                push!(plt2, 1, itrj, plaq)
                end
            end
            push!(plt4, 1, itrj, abs(poly))
            if itrj<500
                push!(plt5, 1, real(poly), imag(poly))
            end
            if 500<itrj<1000
                if itrj%10==0
                    push!(plt5, 1, real(poly), imag(poly))
                end
            end
            if 1000<itrj
                if itrj%50==0
                    push!(plt5, 1, real(poly), imag(poly))
                end
            end
            #
            push!(plt6, 1, itrj, angle(poly))
            #
            plot(plt1, plt2, plt3, plt4, plt5, plt6, layout = 6)
            gui()
        end


        @assert parameters.update_method == "Heatbath"
        verbose = Verbose_1()
        Nsteps = parameters.Nsteps
        numaccepts = 0
        plaq,poly = measurements(0,univ.U,univ,meas;verbose = verbose) # check consistency of preparation.

        for itrj=1:Nsteps
            @time heatbath!(univ)
            plaq,poly = measurements(itrj,univ.U,univ,meas;verbose = verbose)
            plot_refresh!(plt1,plt2,plt3,plt4,plt5,plt6,hist_plaq,hist_poly,hist_poly_θ, plaq, poly, itrj)

            println_verbose1(verbose,"-------------------------------------")
            #println("-------------------------------------")
            flush(stdout)
            flush(verbose)
        end
    end


end