import Pkg
Pkg.activate(".")

using CairoMakie, ColorSchemes, DelimitedFiles

two_column_width = 510
single_column_width = 246
scale=2
height = 150*scale
width = single_column_width*scale
font_size = 10*scale


fig2 = Figure(size = (1.25*width, 2*height))

fontsize_theme = Theme(fontsize=font_size)
set_theme!(fontsize_theme)
Makie.theme(:fonts).:regular = "CMU Serif"

       
ax_sim_D = Axis(fig2[2,1], 
                title=L"$\mathrm{HS\, Simulations}$", 
                ylabel=L"$D$",
                xscale=log10, 
                yscale=log10, 
                xticks = LogTicks(-2:0),
                limits = (10.0^-2, 1.0, nothing, nothing),
                xminorticksvisible=true,
                xminorticks = IntervalsBetween(9),
                yminorticks = IntervalsBetween(9),
                yminorticksvisible=true,
                # xticksmirrored=true,
                xgridvisible=false, 
                ygridvisible=false, 
                xticklabelsvisible=false)


ax_HS_MCT_D = Axis(fig2[2,2], 
                title=L"$\mathrm{HS\, MCT}$", 
                ylabel=L"$D$",
                xscale=log10, 
                yscale=log10, 
                yaxisposition= :right,
                yticks=LogTicks(-6:-1),
                xticks=LogTicks(-2:0),
                limits = (10.0^-3, 10.0^-1, 10.0^-6, 10.0^-1),
                xminorticks = IntervalsBetween(9),
                xminorticksvisible=true,
                yminorticks = IntervalsBetween(9),
                yminorticksvisible=true,
                # xticksmirrored=true,
                xgridvisible=false, 
                ygridvisible=false, 
                xticklabelsvisible=false)
                
ax_sim_NGP_t = Axis(fig2[3,1], 
                xlabel=L"$\varepsilon$", 
                ylabel=L"$t_{\mathrm{peak}}$",
                xscale=log10, 
                yscale=log10, 
                xticks = LogTicks(-2:0),
                limits = (10.0^-2, 1.0, nothing, nothing),
                xminorticksvisible=true,
                xminorticks = IntervalsBetween(9),
                yminorticks = IntervalsBetween(9),
                yminorticksvisible=true,
                xticksmirrored=true,
                xgridvisible=false, 
                ygridvisible=false)

ax_HS_MCT_NGP_t = Axis(fig2[3,2], 
                xlabel=L"$\varepsilon$", 
                ylabel=L"$t_{\mathrm{peak}}$",
                xscale=log10, 
                yscale=log10, 
                yaxisposition= :right,
                limits = (10.0^-3, 10.0^-1, nothing, 10.0^6),
                xticks=LogTicks(-2:0),
                xminorticks = IntervalsBetween(9),
                xminorticksvisible=true,
                yminorticks = IntervalsBetween(9),
                yminorticksvisible=true,
                xticksmirrored=true,
                xgridvisible=false, 
                ygridvisible=false)

# ## extracting from HS MCT 
# MCT_NEP_data = joinpath("DATA", "non_erg_params_PY")
# MCT_sol_data = joinpath("DATA", "MCT_sol_PY")

# MCT_sol_fnames = readdir(MCT_sol_data)

# d_list = list_d(MCT_sol_data)
# sort!(d_list)

colmap1 = cgrad(:matter, 19, categorical = true)
colmap = cgrad([colmap1[i] for i = 2:19], categorical=true)

dir_data = joinpath("PROCESSED_DATA", "ngp_peak_time_diff_const_scaling_data")

fnames = readdir(dir_data)

for fname in fnames 

    if fname == "README.txt"
        continue
    end

    split_fname = split(fname, "_")

    d = parse(Int64, split_fname[end][1:end-4])

    if split_fname[2] == "MCT"

        data = readdlm(joinpath(dir_data, fname), '\t')
        ϵ_arr = data[:,1]
        D_arr = data[:,2]
        t_peak = data[:,3]

        scatter!(ax_HS_MCT_NGP_t, ϵ_arr, t_peak, marker=:utriangle, color=colmap1[d-2])
        scatter!(ax_HS_MCT_D, ϵ_arr, D_arr, marker=:utriangle, color=colmap1[d-2])

    else

        data = readdlm(joinpath(dir_data, fname), '\t')
        ϵ_arr = data[:,1]
        D_arr = data[:,2]
        t_peak = data[:,3]

    scatter!(ax_sim_NGP_t, ϵ_arr, t_peak, color=colmap1[d-2])
    scatter!(ax_sim_D, ϵ_arr, D_arr, color=colmap1[d-2])
    end
end

## Adding decorative straight line for power law divergence of D = ϵ^γ_MF
γ_MF = 2.34
ϵ_arr = [10^(-i) for i in LinRange(1.25,0.1,10)]
lines!(ax_sim_D, ϵ_arr, ϵ_arr.^(γ_MF), color= :black, linestyle= :dash)
text!(ax_sim_D, 0.125, 2*10.0^-3, text=L"$D \propto \varepsilon^{\gamma_{\mathrm{MF}}}$", fontsize=20)

## Adding decorative straight line for power law divergence of t_peak = ϵ^-γ_MF
γ_MF = 2.34
ϵ_arr = [10^(-i) for i in LinRange(1.25,0.1,10)]
lines!(ax_sim_NGP_t, ϵ_arr, 0.11*ϵ_arr.^(-γ_MF), color= :black, linestyle= :dash)
text!(ax_sim_NGP_t, 0.85*ϵ_arr[end-6], 1.1*ϵ_arr[end-4]^(-2.0), text=L"$t_{\mathrm{peak}} \propto \varepsilon^{-\gamma_{\mathrm{MF}}}$", fontsize=font_size)


## Adding Labels
marks = []
labs = []

d_list = [3,4,6,8,10,12]

for d in d_list
    mark = MarkerElement(color = colmap1[d-2], marker = :circle, markersize = 12)
    lab = L"$d=%$(d)$"

    push!(marks, mark)
    push!(labs, lab)
end


Legend(fig2[1, 1:2], marks, labs,
    patchsize = (30, 10), rowgap = 0, framevisible=false, labelsize=font_size, nbanks = 3, tellheight = true)
    
for (ax, label) in zip([ax_sim_D, ax_HS_MCT_D, ax_sim_NGP_t, ax_HS_MCT_NGP_t], ["(a)", "(b)", "(c)", "(d)"])
    text!(
        ax, 1, 1,
        text = label,
        align = (:right, :top),
        offset = (-3, -1),
        space = :relative,
        fontsize=font_size
    )
end

display(fig2)

figname2 = "vanishing_D_growing_t_peak_HS_sim_MCT.pdf"
save(joinpath("PLOTS", figname2), fig2)