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


fig = Figure(size = (1.25*width, height))

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

ax_HS_MCT_NGP_peak = Axis(fig[1,2], 
            xlabel=L"$\varepsilon$", 
            ylabel=L"$\alpha_2^*(\varphi ; d)$",
            xscale=log10, 
            yscale=log10,
            limits=(10.0^-3, 10.0^-1, nothing, nothing),
            xticks=LogTicks(-2:0),
            title=L"\mathrm{HS\, \, MCT}", 
            xminorticks = IntervalsBetween(9),
            yminorticks = IntervalsBetween(9),
            yticks = LogTicks(-2:0),
            yaxisposition= :right,
            xminorticksvisible=true,
            yminorticksvisible=true,
            xgridvisible=false, 
            ygridvisible=false)


ax_sim_NGP_peak = Axis(fig[1,1], 
                xlabel=L"$\varepsilon$", 
                ylabel=L"$\alpha_2^*(\varphi ; d)$",
                xscale=log10, 
                yscale=log10,
                title=L"\mathrm{HS\, \,Simulations}", 
                limits = (10.0^-2, 1.0, nothing, nothing),
                xticks = LogTicks(-2:0),
                yticks = LogTicks(-2:0),
                xminorticks = IntervalsBetween(9),
                yminorticks = IntervalsBetween(9),
                xminorticksvisible=true,
                yminorticksvisible=true,
                xgridvisible=false, 
                ygridvisible=false)

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

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

fnames = readdir(dir_data)

for fname in fnames 
    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]
        α2_peak = data[:,2]

        scatter!(ax_HS_MCT_NGP_peak, ϵ_arr, α2_peak, marker=:utriangle, color=colmap1[d-2])
    else

        data = readdlm(joinpath(dir_data, fname), '\t')
        ϵ_arr = data[:,1]
        α2_peak = data[:,2]

        scatter!(ax_sim_NGP_peak, ϵ_arr, α2_peak, color=colmap1[d-2])
    end
end

## Adding decorative straight line for power law divergence of α2_sat = ϵ^-1
ϵ_arr = [10^(-i) for i in LinRange(0.9,0.33,10)]
lines!(ax_sim_NGP_peak, ϵ_arr, 0.03*ϵ_arr.^(-1), color= :black, linestyle= :dash)
text!(ax_sim_NGP_peak, 0.1, 3.5*10.0^-1, 
      text=L"$\alpha_2^{*}(\varphi ; d) \propto \varepsilon^{-1}$", 
      fontsize=font_size, rotation = -0.8)

## 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

for (ax, label) in zip([ax_sim_NGP_peak,ax_HS_MCT_NGP_peak], ["(a)", "(b)"])
    text!(
        ax, 1, 1,
        text = label,
        align = (:right, :top),
        offset = (-3, -1),
        space = :relative,
        fontsize=font_size
    )
end
    
display(fig)

figname = "HS_NGP_scaling.pdf"
save(joinpath("PLOTS", figname), fig)

########################################
########################################

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

fig = Figure(size = (width, height))

ax_HS_NGP_peak = Axis(fig[1,1], 
                        xlabel=L"$\varepsilon$", 
                        # xticklabelsvisible=false,
                        ylabel=L"$\hat{\alpha}_2^*(\varphi ; d)$",
                        # title=L"$\mathrm{HS\, Simulations}$",
                        xscale=log10, 
                        yscale=log10,
                        yticks=LogTicks(-1:1:2),
                        xticks=LogTicks(-2:1:2),
                        xminorticks = IntervalsBetween(9),
                        yminorticks = IntervalsBetween(9),
                        xminorticksvisible=true,
                        yminorticksvisible=true,
                        xgridvisible=false, 
                        ygridvisible=false, 
                        limits=(nothing, nothing, 0.1, 20.0))

dir = joinpath("PROCESSED_DATA", "dimscaled_ngp_peak_scaling_data")

fnames = readdir(dir)   

for fname in fnames
    split_fname = split(fname, "_")

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

    if split_fname[1] == "HS"
        
        data = readdlm(joinpath(dir, fname), '\t')

        ϵ_arr_data = data[:,1]
        dimscaled_α2_peak = data[:,2]

        scatter!(ax_HS_NGP_peak, ϵ_arr_data, dimscaled_α2_peak, marker=:circle, color=colmap1[d-2])

    elseif split_fname[1] == "MK"

        data = readdlm(joinpath(dir, fname), '\t')

        ϵ_arr_data = data[:,1]
        dimscaled_α2_peak = data[:,2]

        scatter!(ax_HS_NGP_peak, ϵ_arr_data, dimscaled_α2_peak, marker=:circle, color=(colmap1[d-2],0.0), strokewidth=2, strokecolor=colmap1[d])
    end

end

d_list = [3,4,6,8,10,12]
marks = []
labs = []
for d in d_list
    mark = MarkerElement(color = colmap1[d-2], marker = :circle, markersize = 9)
    lab = L"$d=%$(d)$"

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


Legend(fig[1, 2], marks, labs, rowgap = 0, framevisible=false, labelsize=font_size, nbanks = 1)

marks = [MarkerElement(marker=:circle, color=(:black, 0), strokewidth=2, strokecolor=:black), MarkerElement(marker=:circle, color=:black)]
labs = [L"$\mathrm{MK}$", L"$\mathrm{HS}$"]

axislegend(ax_HS_NGP_peak, marks, labs, position=:lb, framevisible=false, labelsize=20)

display(fig)

figname = "MK_NGP_scaling.pdf"
save(joinpath("PLOTS", figname), fig)