text:
using Plots
using StatsPlots
using DelimitedFiles
using Plots.PlotMeasures
using DocStringExtensions
using StatsBase
relative_error(future, initial) = (future - initial) ./ initial
include("compare_Present_Future_low_flows.jl")
include("loadfunctions.jl")

Area_Catchment_Gailtal = sum([98227533.0, 184294158.0, 83478138.0, 220613195.0])
Area_Catchment_Palten = sum([198175943.0, 56544073.0, 115284451.3])
Area_Catchment_Pitten = 115496400.
Area_Catchment_Silbertal = 100139168.
Area_Catchment_Defreggental = sum([235811198.0, 31497403.0])
Area_Catchment_Pitztal = sum([20651736.0, 145191864.0])

Catchment_Names = ["Pitten", "Palten", "Gailtal", "IllSugadin", "Defreggental", "Pitztal"]
Catchment_Height = [917, 1315, 1476, 1776, 2233, 2558]
Area_Catchments = [Area_Catchment_Pitten, Area_Catchment_Palten, Area_Catchment_Gailtal, Area_Catchment_Silbertal, Area_Catchment_Defreggental, Area_Catchment_Pitztal]
nr_runs = [300,300,298,300, 300, 300]

function plot_changes_monthly_discharge_all_catchments_absolute_1proj(All_Catchment_Names, Elevation, Area_Catchments, nr_runs, proj_nr)
    xaxis_45 = collect(1:2:23)
    xaxis_85 = collect(2:2:24)
    # ----------------- Plot Absolute Change ------------------
    plot()
    Farben_85 = palette(:reds)
    Farben_45 = palette(:blues)
    all_boxplots = []
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        first_index = nr_runs[i] * 12 * (proj_nr-1) + 1
        last_index = nr_runs[i] *12* proj_nr
        monthly_changes_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_8.5.txt", ',')
        months_85 = monthly_changes_85[first_index:last_index,1]
        Monthly_Discharge_past_85 = monthly_changes_85[first_index:last_index,2]
        Monthly_Discharge_future_85  = monthly_changes_85[first_index:last_index,3]
        monthly_changes_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_4.5.txt", ',')
        months_45 = monthly_changes_45[first_index:last_index,1]
        Monthly_Discharge_past_45 = monthly_changes_45[first_index:last_index,2]
        Monthly_Discharge_future_45  = monthly_changes_45[first_index:last_index,3]
        Monthly_Discharge_Change_45  = monthly_changes_45[first_index:last_index,4]
        Monthly_Discharge_past_45 = convertDischarge(Monthly_Discharge_past_45, Area_Catchments[i])
        Monthly_Discharge_past_85 = convertDischarge(Monthly_Discharge_past_85, Area_Catchments[i])
        Monthly_Discharge_future_45 = convertDischarge(Monthly_Discharge_future_45, Area_Catchments[i])
        Monthly_Discharge_future_85 = convertDischarge(Monthly_Discharge_future_85, Area_Catchments[i])
        plot()
        box = []
        for month in 1:12
            boxplot!([xaxis_45[month]],Monthly_Discharge_future_45[findall(x-> x == month, months_45)] - Monthly_Discharge_past_45[findall(x-> x == month, months_45)] , size=(2000,800), leg=false, color=["blue"], alpha=0.8, label="RCP 4.5", outlier=false)
            boxplot!([xaxis_85[month]],Monthly_Discharge_future_85[findall(x-> x == month, months_45)] - Monthly_Discharge_past_85[findall(x-> x == month, months_45)], size=(2000,800), leg=false, color=["red"], left_margin = [5mm 0mm], minorticks = true, gridlinewidth=2, framestyle = :box, label="RCP 8.5", outlier=false)
        end
        ylabel!("[mm/d]", yguidefontsize=20)
        #title!("Relative Change in Discharge RCP 4.5 =blue, RCP 4.5 = red")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        boxplot!(left_margin = [5mm 0mm], bottom_margin = 20px, xtickfont = font(20), ytickfont = font(20), outlier=false)

        ylims!((-3.5,3.5))
        yticks!([-3:1:3;])

        hline!([0], color=["grey"], linestyle = :dash)
        #hline!([100], color=["grey"], linestyle = :dash)
        #hline!([50], color=["grey"], linestyle = :dash)
        #hline!([-25], color=["grey"], linestyle = :dash)
        xticks!([1.5:2:23.5;], ["Jan", "Feb", "Mar", "Apr", "May","Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
        box = boxplot!()
        push!(all_boxplots, box)
    end
    plot(all_boxplots[1], all_boxplots[2], all_boxplots[3], all_boxplots[4], all_boxplots[5], all_boxplots[6], layout= (3,2), legend = false, size=(2000,1500), left_margin = [5mm 0mm], bottom_margin = 20px, yguidefontsize=20, xtickfont = font(20), ytickfont = font(20), dpi=300, minorgrid=true, gridlinewidth=4, minorgridlinewidth=2)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/monthly_discharges_all_catchments_absolute_change_proj_"*string(proj_nr)*".png")
end

function plot_changes_monthly_discharge_all_catchments_1proj(All_Catchment_Names, Elevation, nr_runs, proj_nr)
    xaxis_45 = collect(1:2:23)
    xaxis_85 = collect(2:2:24)
    # ----------------- Plot Absolute Change ------------------
    plot()
    Farben_85 = palette(:reds)
    Farben_45 = palette(:blues)
    all_boxplots = []
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        first_index = nr_runs[i] * 12 * (proj_nr-1) + 1
        last_index = nr_runs[i] *12* proj_nr
        monthly_changes_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_8.5.txt", ',')
        months_85 = monthly_changes_85[first_index:last_index,1]
        Monthly_Discharge_past_85 = monthly_changes_85[first_index:last_index,2]
        Monthly_Discharge_future_85  = monthly_changes_85[first_index:last_index,3]
        monthly_changes_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_4.5.txt", ',')
        months_45 = monthly_changes_45[first_index:last_index,1]
        Monthly_Discharge_past_45 = monthly_changes_45[first_index:last_index,2]
        Monthly_Discharge_future_45  = monthly_changes_45[first_index:last_index,3]
        Monthly_Discharge_Change_45  = monthly_changes_45[first_index:last_index,4]
        Monthly_Discharge_past_45 = convertDischarge(Monthly_Discharge_past_45, Area_Catchments[i])
        Monthly_Discharge_past_85 = convertDischarge(Monthly_Discharge_past_85, Area_Catchments[i])
        Monthly_Discharge_future_45 = convertDischarge(Monthly_Discharge_future_45, Area_Catchments[i])
        Monthly_Discharge_future_85 = convertDischarge(Monthly_Discharge_future_85, Area_Catchments[i])
        plot()
        box = []
        for month in 1:12
            boxplot!([xaxis_45[month]],relative_error(Monthly_Discharge_future_45[findall(x-> x == month, months_45)], Monthly_Discharge_past_45[findall(x-> x == month, months_45)])*100, size=(2000,800), leg=false, color=["blue"], alpha=0.8)
            boxplot!([xaxis_85[month]],relative_error(Monthly_Discharge_future_85[findall(x-> x == month, months_45)], Monthly_Discharge_past_85[findall(x-> x == month, months_45)])*100, size=(2000,800), leg=false, color=["red"], left_margin = [5mm 0mm], minorticks = true, gridlinewidth=2, framestyle = :box)
        end
        ylabel!("[%]", yguidefontsize=20)
        #title!("Relative Change in Discharge RCP 4.5 =blue, RCP 4.5 = red")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        boxplot!(left_margin = [5mm 0mm], bottom_margin = 20px, xtickfont = font(20), ytickfont = font(20))
        if Catchment_Name == Catchment_Name == "Pitztal"# || "Defreggental"
            ylims!((-100,850))
            yticks!([-100:100:750;])
        elseif Catchment_Name == "Pitten" || Catchment_Name == "Palten"
            ylims!((-100,100))
            yticks!([-100:50:100;])
        else
            ylims!((-100,275))
            yticks!([-100:50:275;])
        end
        hline!([0], color=["grey"], linestyle = :dash)
        #hline!([100], color=["grey"], linestyle = :dash)
        #hline!([50], color=["grey"], linestyle = :dash)
        #hline!([-25], color=["grey"], linestyle = :dash)
        xticks!([1.5:2:23.5;], ["Jan", "Feb", "Mar", "Apr", "May","Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
        box = boxplot!()
        push!(all_boxplots, box)
    end
    plot(all_boxplots[1], all_boxplots[2], all_boxplots[3], all_boxplots[4], all_boxplots[5], all_boxplots[6], layout= (3,2), legend = false, size=(2000,1500), left_margin = [5mm 0mm], bottom_margin = 20px, yguidefontsize=20, xtickfont = font(20), ytickfont = font(20), dpi=300, minorgrid=true, gridlinewidth=4, minorgridlinewidth=2)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/monthly_discharges_all_catchments_proj"*string(proj_nr)*".png")
end

# plot_changes_monthly_discharge_all_catchments_absolute_1proj(Catchment_Names, Catchment_Height, Area_Catchments, nr_runs, 12)
# plot_changes_monthly_discharge_all_catchments_1proj(Catchment_Names, Catchment_Height, nr_runs, 1)

function plot_ecdf_timing_AMF_proj(All_Catchment_Names, Elevation, time)
    Farben_proj = palette(:tab20)
    all_plots = []
    highest_temp_proj_85 = [12,12,12,12,12,12]
    lowest_temp_proj_85 = [1,1,1,6,1,1]
    highest_prec_proj_85 = [8,8,8,8,12,2]
    lowest_prec_proj_85 = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        annual_max_flow_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_4.5.txt",',')
        annual_max_flow_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_8.5.txt",',')
        average_max_Discharge_past_45 = annual_max_flow_45[:,1]
        average_max_Discharge_future_45 = annual_max_flow_45[:,2]
        Timing_max_Discharge_past_45 = annual_max_flow_45[:,3]
        Timing_max_Discharge_future_45 = annual_max_flow_45[:,4]
        average_max_Discharge_past_85 = annual_max_flow_85[:,1]
        average_max_Discharge_future_85 = annual_max_flow_85[:,2]
        Timing_max_Discharge_past_85 = annual_max_flow_85[:,3]
        Timing_max_Discharge_future_85 = annual_max_flow_85[:,4]
        # now for each projection
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection
            first_index = nr_runs[i] * (proj-1) + 1
            last_index = nr_runs[i] * proj
            if proj == highest_temp_proj_85[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj_85[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj_85[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj_85[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(Timing_max_Discharge_past_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Timing_max_Discharge_future_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "change"
                plot!(StatsBase.ecdf(Timing_max_Discharge_future_85[first_index:last_index] - Timing_max_Discharge_past_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Day of Year")
        if time == "past" || time =="future"
            xlims!((1,365))
        elseif time == "change" && Catchment_Name == "Gailtal"
            xlims!((-300,150))
        elseif time == "change" && Catchment_Name != "Gailtal"
            xlims!((-100,150))
        end
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15), minorticks=true, minorgrid=true, minorgridlinewidth=2)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_timing_AMF_"*time*"_85.png")
end

function plot_ecdf_magnitude_AMF_proj(All_Catchment_Names, Area_Catchments, Elevation, time)
    Farben_proj = palette(:tab20)
    all_plots = []
    highest_temp_proj_85 = [12,12,12,12,12,12]
    lowest_temp_proj_85 = [1,1,1,6,1,1]
    highest_prec_proj_85 = [8,8,8,8,12,2]
    lowest_prec_proj_85 = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        annual_max_flow_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_4.5.txt",',')
        annual_max_flow_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_8.5.txt",',')
        average_max_Discharge_past_45 = convertDischarge(annual_max_flow_45[:,1], Area_Catchments[i])
        average_max_Discharge_future_45 = convertDischarge(annual_max_flow_45[:,2], Area_Catchments[i])
        average_max_Discharge_past_85 = convertDischarge(annual_max_flow_85[:,1], Area_Catchments[i])
        average_max_Discharge_future_85 = convertDischarge(annual_max_flow_85[:,2], Area_Catchments[i])
        # now for each projection
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection
            first_index = nr_runs[i] * (proj-1) + 1
            last_index = nr_runs[i] * proj
            if proj == highest_temp_proj[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(average_max_Discharge_past_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(average_max_Discharge_future_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "change"
                plot!(StatsBase.ecdf(relative_error(average_max_Discharge_future_45[first_index:last_index], average_max_Discharge_past_45[first_index:last_index])*100), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("[%]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15), minorticks=true, minorgrid=true, minorgridlinewidth=2)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/plot_ecdf_magnitude_AMF_rel_"*time*"_45.png")
end

function plot_ecdf_timing_AMF_proj_all_years(All_Catchment_Names, Elevation, time)
    Farben_proj = palette(:tab20)
    all_plots = []
    all_plots_85 = []
    highest_temp_proj_85 = [12,12,12,12,12,12]
    lowest_temp_proj_85 = [1,1,1,6,1,1]
    highest_prec_proj_85 = [8,8,8,8,12,2]
    lowest_prec_proj_85 = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        max_discharge_prob_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_4.5.txt", ',')
        max_Discharge_Past_45 = max_discharge_prob_45[:,1]
        Max_Discharge_Future_45 = max_discharge_prob_45[:,2]
        Exceedance_Probability_45 = max_discharge_prob_45[:,3]
        Date_Past_45 = max_discharge_prob_45[:,4]
        Date_Future_45 = max_discharge_prob_45[:,5]

        max_discharge_prob_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_8.5.txt", ',')
        max_Discharge_Past_85 = max_discharge_prob_85[:,1]
        Max_Discharge_Future_85 = max_discharge_prob_85[:,2]
        Exceedance_Probability_85 = max_discharge_prob_85[:,3]
        Date_Past_85 = max_discharge_prob_85[:,4]
        Date_Future_85 = max_discharge_prob_85[:,5]
        # now for each projection
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *30* (proj-1) + 1
            last_index = nr_runs[i] *30* proj
            if proj == highest_temp_proj[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(Date_Past_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Date_Future_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Day of Year")
        if time == "past" || time =="future"
            xlims!((1,365))
        elseif time == "change" && Catchment_Name == "Gailtal"
            xlims!((-300,150))
        elseif time == "change" && Catchment_Name != "Gailtal"
            xlims!((-100,150))
        end
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
        # for RCP 8.5
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *30* (proj-1) + 1
            last_index = nr_runs[i] *30* proj
            if proj == highest_temp_proj_85[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj_85[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj_85[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj_85[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(Date_Past_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Date_Future_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Day of Year")
        if time == "past" || time =="future"
            xlims!((1,365))
        elseif time == "change" && Catchment_Name == "Gailtal"
            xlims!((-300,150))
        elseif time == "change" && Catchment_Name != "Gailtal"
            xlims!((-100,150))
        end
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots_85, plot!(framestyle = :box))

    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_timing_AMF_all_years"*time*"_45.png")
    plot(all_plots_85[1], all_plots_85[2], all_plots_85[3], all_plots_85[4], all_plots_85[5], all_plots_85[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_timing_AMF_all_years"*time*"_85.png")
end

function plot_ecdf_magnitude_AMF_proj_all_years(All_Catchment_Names, Area_Catchments, Elevation, time)
    Farben_proj = palette(:tab20)
    all_plots = []
    all_plots_85 = []
    highest_temp_proj_85 = [12,12,12,12,12,12]
    lowest_temp_proj_85 = [1,1,1,6,1,1]
    highest_prec_proj_85 = [8,8,8,8,12,2]
    lowest_prec_proj_85 = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        max_discharge_prob_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_4.5.txt", ',')
        max_Discharge_Past_45 = convertDischarge(max_discharge_prob_45[:,1], Area_Catchments[i])
        Max_Discharge_Future_45 = convertDischarge(max_discharge_prob_45[:,2], Area_Catchments[i])

        max_discharge_prob_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/change_max_Annual_Discharge_prob_distr_8.5.txt", ',')
        max_Discharge_Past_85 = convertDischarge(max_discharge_prob_85[:,1], Area_Catchments[i])
        Max_Discharge_Future_85 = convertDischarge(max_discharge_prob_85[:,2], Area_Catchments[i])
        # now for each projection
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *30* (proj-1) + 1
            last_index = nr_runs[i] *30* proj
            if proj == highest_temp_proj[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(max_Discharge_Past_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Max_Discharge_Future_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Magnitude [mm/d]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
        # for RCP 8.5
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *30* (proj-1) + 1
            last_index = nr_runs[i] *30* proj
            if proj == highest_temp_proj_85[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj_85[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj_85[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj_85[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(max_Discharge_Past_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Max_Discharge_Future_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Magnitude [mm/d]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots_85, plot!(framestyle = :box))

    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_magnitude_AMF_all_years_"*time*"_45.png")
    plot(all_plots_85[1], all_plots_85[2], all_plots_85[3], all_plots_85[4], all_plots_85[5], all_plots_85[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_magnitude_AMF_all_years_"*time*"_85.png")
end

#plot_ecdf_timing_AMF_proj(Catchment_Names, Catchment_Height, "change")
#plot_ecdf_magnitude_AMF_proj(Catchment_Names, Area_Catchments, Catchment_Height, "change")
# plot_ecdf_timing_AMF_proj_all_years(Catchment_Names, Catchment_Height, "future")
# plot_ecdf_magnitude_AMF_proj_all_years(Catchment_Names, Area_Catchments, Catchment_Height, "future")

function plot_ecdf_timing_min_flow_proj_all_years(All_Catchment_Names, Elevation, time)
    Farben_proj = palette(:tab20)
    all_plots = []
    all_plots_85 = []
    #for RCP 8.5
    highest_temp_proj_85 = [12,12,12,12,12,12]
    lowest_temp_proj_85 = [1,1,1,6,1,1]
    highest_prec_proj_85 = [8,8,8,8,12,2]
    lowest_prec_proj_85 = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        max_discharge_prob_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_prob_distr_4.5.txt", ',')
        Farben45 = palette(:blues)
        max_discharge_prob_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_prob_distr_8.5.txt", ',')
        Farben85 = palette(:reds)
        # shift dates so that 1. June (day 152) is first day, therefore all days below 151 have to be shifted by 216 and all above have to be shifted by -152
        Date_Past_45 = max_discharge_prob_45[:,4]
        Date_Future_45 = max_discharge_prob_45[:,5]
        Date_Past_85 = max_discharge_prob_85[:,4]
        Date_Future_85 = max_discharge_prob_85[:,5]
        Date_Past_45_shifted = zeros(nr_runs[i] *29* 14)
        Date_Past_45_shifted[findall(x->x >= 152, Date_Past_45)] = Date_Past_45[findall(x->x >= 152, Date_Past_45)] .- 151
        Date_Past_45_shifted[findall(x->x < 152, Date_Past_45)] = Date_Past_45[findall(x->x < 152, Date_Past_45)] .+ 215
        Date_Future_45_shifted = zeros(nr_runs[i] *29* 14)
        Date_Future_45_shifted[findall(x->x >= 152, Date_Future_45)] = Date_Future_45[findall(x->x >= 152, Date_Future_45)] .- 151
        Date_Future_45_shifted[findall(x->x < 152, Date_Future_45)] = Date_Future_45[findall(x->x < 152, Date_Future_45)] .+ 215
        Date_Past_85_shifted = zeros(nr_runs[i] *29* 14)
        Date_Past_85_shifted[findall(x->x >= 152, Date_Past_85)] = Date_Past_85[findall(x->x >= 152, Date_Past_85)] .- 151
        Date_Past_85_shifted[findall(x->x < 152, Date_Past_85)] = Date_Past_85[findall(x->x < 152, Date_Past_85)] .+ 215
        Date_Future_85_shifted = zeros(nr_runs[i] *29* 14)
        Date_Future_85_shifted[findall(x->x >= 152, Date_Future_85)] = Date_Future_85[findall(x->x >= 152, Date_Future_85)] .- 151
        Date_Future_85_shifted[findall(x->x < 152, Date_Future_85)] = Date_Future_85[findall(x->x < 152, Date_Future_85)] .+ 215
        # now for each projection
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *29* (proj-1) + 1
            last_index = nr_runs[i] *29* proj
            if proj == highest_temp_proj[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(Date_Past_45_shifted[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Date_Future_45_shifted[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Day of Year")
        xlims!((1,365))
        xticks!([49,99,149,199,265,315, 355], ["200", "250", "300", "350", "50", "100", "150"])
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
        # for RCP 8.5
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *29* (proj-1) + 1
            last_index = nr_runs[i] *29* proj
            if proj == highest_temp_proj_85[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj_85[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj_85[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj_85[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(Date_Past_85_shifted[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Date_Future_85_shifted[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Day of Year")
        xlims!((1,365))
        xticks!([49,99,149,199,265,315, 355], ["200", "250", "300", "350", "50", "100", "150"])
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots_85, plot!(framestyle = :box))

    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/Low_Flows/ECDF_timing_min_flow_all_years"*time*"_45_new.png")
    plot(all_plots_85[1], all_plots_85[2], all_plots_85[3], all_plots_85[4], all_plots_85[5], all_plots_85[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/Low_Flows/ECDF_timing_min_flow_all_years"*time*"_85_new.png")
end

function plot_ecdf_magnitude_min_flow_proj_all_years(All_Catchment_Names, Area_Catchments, Elevation, time)
    Farben_proj = palette(:tab20)
    all_plots = []
    all_plots_85 = []
    highest_temp_proj_85 = [12,12,12,12,12,12]
    lowest_temp_proj_85 = [1,1,1,6,1,1]
    highest_prec_proj_85 = [8,8,8,8,12,2]
    lowest_prec_proj_85 = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        max_discharge_prob_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_prob_distr_4.5.txt", ',')
        Farben45 = palette(:blues)
        max_discharge_prob_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/LowFlows/low_flows_prob_distr_8.5.txt", ',')
        Farben85 = palette(:reds)
        max_Discharge_Past_45 = max_discharge_prob_45[:,1]
        Max_Discharge_Future_45 = max_discharge_prob_45[:,2]
        Exceedance_Probability_45 = max_discharge_prob_45[:,3]
        max_Discharge_Past_85 = max_discharge_prob_85[:,1]
        Max_Discharge_Future_85 = max_discharge_prob_85[:,2]
        Exceedance_Probability_85 = max_discharge_prob_85[:,3]
        # now for each projection
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *29* (proj-1) + 1
            last_index = nr_runs[i] *29* proj
            if proj == highest_temp_proj[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(max_Discharge_Past_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Max_Discharge_Future_45[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Magnitude [mm/d]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
        # for RCP 8.5
        plot()
        for proj in 1:14
            # get the indices corresponding to this projection, values for 30 years
            first_index = nr_runs[i] *29* (proj-1) + 1
            last_index = nr_runs[i] *29* proj
            if proj == highest_temp_proj_85[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj_85[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj_85[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj_85[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            if time == "past"
                plot!(StatsBase.ecdf(max_Discharge_Past_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Max_Discharge_Future_85[first_index:last_index]), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Magnitude [mm/d]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots_85, plot!(framestyle = :box))

    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/Low_Flows/ECDF_magnitude_min_flow_all_years_"*time*"_45.png")
    plot(all_plots_85[1], all_plots_85[2], all_plots_85[3], all_plots_85[4], all_plots_85[5], all_plots_85[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/Low_Flows/ECDF_magnitude_min_flow_all_years_"*time*"_85.png")
end

#plot_ecdf_timing_min_flow_proj_all_years(Catchment_Names, Catchment_Height, "future")
#plot_ecdf_magnitude_min_flow_proj_all_years(Catchment_Names, Area_Catchments, Catchment_Height, "future")

function plot_ecdf_prec_temp_all_catchments(All_Catchment_Names, Elevation, nr_runs, time)
    plot()
    Farben_85 = palette(:reds)
    Farben_45 = palette(:blues)
    Farben_proj = palette(:tab20)
    all_plots = []
    all_plots_85 = []
    for (h,Catchment_Name) in enumerate(All_Catchment_Names)
        Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
        Timeseries_End = readdlm("/home/sarah/Master/Thesis/Data/Projektionen/End_Timeseries_45_85.txt",',')
        path_45 = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/"
        path_85 = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp85/rcp85/"
        Name_Projections_45 = readdir(path_45)
        Name_Projections_85 = readdir(path_85)
        #all_months_all_runs = Float64[]
        average_monthly_Precipitation_past45 = Float64[]
        average_monthly_Precipitation_future45 = Float64[]
        average_monthly_Precipitation_past85 = Float64[]
        average_monthly_Precipitation_future85 = Float64[]
        average_monthly_Temperature_past45 = Float64[]
        average_monthly_Temperature_past85 = Float64[]
        average_monthly_Temperature_future45 = Float64[]
        average_monthly_Temperature_future85 = Float64[]

        if Catchment_Name == "Gailtal"
            ID_Prec_Zones = [113589, 113597, 113670, 114538]
            # size of the area of precipitation zones
            Area_Zones = [98227533.0, 184294158.0, 83478138.0, 220613195.0]
            Temp_Elevation = 1140.0
            Mean_Elevation_Catchment = 1500
            ID_temp = 113597
            Elevations_Catchment = Elevations(200.0, 400.0, 2800.0, Temp_Elevation, Temp_Elevation)
        elseif Catchment_Name == "Palten"
            ID_Prec_Zones = [106120, 111815, 9900]
            Area_Zones = [198175943.0, 56544073.0, 115284451.3]
            ID_temp = 106120
            Temp_Elevation = 1265.0
            Mean_Elevation_Catchment = 1300 # in reality 1314
            Elevations_Catchment = Elevations(200.0, 600.0, 2600.0, Temp_Elevation, Temp_Elevation)
        elseif Catchment_Name == "Pitten"
            ID_Prec_Zones = [109967]
            Area_Zones = [115496400.]
            ID_temp = 10510
            Mean_Elevation_Catchment = 900 # in reality 917
            Temp_Elevation = 488.0
            Elevations_Catchment = Elevations(200.0, 400.0, 1600.0, Temp_Elevation, Temp_Elevation)
        elseif Catchment_Name == "Defreggental"
            ID_Prec_Zones = [17700, 114926]
            Area_Zones = [235811198.0, 31497403.0]
            ID_temp = 17700
            Mean_Elevation_Catchment =  2300 # in reality 2233.399986
            Temp_Elevation = 1385.
            Elevations_Catchment = Elevations(200.0, 1000.0, 3600.0, Temp_Elevation, Temp_Elevation)
        elseif Catchment_Name == "IllSugadin"
            ID_Prec_Zones = [100206]
            Area_Zones = [100139168.]
            ID_temp = 14200
            Mean_Elevation_Catchment = 1700
            Temp_Elevation = 670.
            Elevations_Catchment = Elevations(200.0, 600.0, 2800.0, Temp_Elevation, Temp_Elevation)
        elseif Catchment_Name == "Pitztal"
            ID_Prec_Zones = [102061, 102046]
            Area_Zones = [20651736.0, 145191864.0]
            ID_temp = 14620
            Mean_Elevation_Catchment =  2500 # in reality 2233.399986
            Temp_Elevation = 1410.
            Elevations_Catchment = Elevations(200.0, 1200.0, 3800.0, Temp_Elevation, Temp_Elevation)
        end
        Area_Catchment = sum(Area_Zones)
        Area_Zones_Percent = Area_Zones / Area_Catchment
        plot()
        for (i, name) in enumerate(Name_Projections_45)
            Timeseries_Future = collect(Date(Timeseries_End[i,1]-29,1,1):Day(1):Date(Timeseries_End[i,1],12,31))
            #print(size(Timeseries_Past), size(Timeseries_Future))
            Timeseries_Proj = readdlm(path_45*name*"/"*Catchment_Name*"/pr_model_timeseries.txt")
            Timeseries_Proj = Date.(Timeseries_Proj, Dates.DateFormat("y,m,d"))
            Temperature = readdlm(path_45*name*"/"*Catchment_Name*"/tas_"*string(ID_temp)*"_sim1.txt", ',')[:,1]
            Elevation_Zone_Catchment, Temperature_Elevation_Catchment, Total_Elevationbands_Catchment = gettemperatureatelevation(Elevations_Catchment, Temperature)
            # get the temperature data at the mean elevation to calculate the mean potential evaporation
            Temperature = Temperature_Elevation_Catchment[:,findfirst(x-> x==Mean_Elevation_Catchment, Elevation_Zone_Catchment)]

            indexstart_past = findfirst(x-> x == Dates.year(Timeseries_Past[1]), Dates.year.(Timeseries_Proj))[1]
            indexend_past = findlast(x-> x == Dates.year(Timeseries_Past[end]), Dates.year.(Timeseries_Proj))[1]
            Temperature_Past = Temperature[indexstart_past:indexend_past] ./ 10
            #print(Dates.year(Timeseries_Future[end]), Dates.year.(Timeseries_Proj[end]))
            indexstart_future = findfirst(x-> x == Dates.year(Timeseries_Future[1]), Dates.year.(Timeseries_Proj))[1]
            indexend_future = findlast(x-> x == Dates.year(Timeseries_Future[end]), Dates.year.(Timeseries_Proj))[1]
            Temperature_Future = Temperature[indexstart_future:indexend_future] ./ 10
            # calculate monthly mean temperature
            Monthly_Temperature_Past, Month = monthly_discharge(Temperature_Past, Timeseries_Past)
            Monthly_Temperature_Future, Month_future = monthly_discharge(Temperature_Future, Timeseries_Future)

            #-------- PRECIPITATION ------------------
            Precipitation_All_Zones = Array{Float64, 1}[]
            Total_Precipitation_Proj = zeros(length(Timeseries_Proj))
            for j in 1: length(ID_Prec_Zones)
                    # get precipitation projections for the precipitation measurement
                    Precipitation_Zone = readdlm(path_45*name*"/"*Catchment_Name*"/pr_"*string(ID_Prec_Zones[j])*"_sim1.txt", ',')[:,1]
                    #print(size(Precipitation_Zone), typeof(Precipitation_Zone))
                    push!(Precipitation_All_Zones, Precipitation_Zone ./10)
                    Total_Precipitation_Proj += Precipitation_All_Zones[j].*Area_Zones_Percent[j]
            end
            #Total_Precipitation_Proj = Precipitation_All_Zones[1].*Area_Zones_Percent[1] + Precipitation_All_Zones[2].*Area_Zones_Percent[2] + Precipitation_All_Zones[3].*Area_Zones_Percent[3] + Precipitation_All_Zones[4].*Area_Zones_Percent[4]
            Precipitation_Past = Total_Precipitation_Proj[indexstart_past:indexend_past]
            Precipitation_Future = Total_Precipitation_Proj[indexstart_future:indexend_future]

            # Monthly_Precipitation_Past, Month = monthly_precipitation(Precipitation_Past, Timeseries_Past)
            # Monthly_Precipitation_Future, Month_future = monthly_precipitation(Precipitation_Future, Timeseries_Future)
            # append!(average_monthly_Temperature_past45, mean(Monthly_Temperature_Past))
            # append!(average_monthly_Temperature_future45, mean(Monthly_Temperature_Future))
            # append!(average_monthly_Precipitation_past45, mean(Monthly_Precipitation_Past)*12)
            # append!(average_monthly_Precipitation_future45, mean(Monthly_Precipitation_Future)*12)

            # make ecdf plots
            if time == "past"
                plot!(StatsBase.ecdf(Precipitation_Past), label="Proj "*string(i), color=[Farben_proj[i]], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Precipitation_Future), label="Proj "*string(i), color=[Farben_proj[i]], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Precipitation [mm/d]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[h])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[h])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[h])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[h])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
        plot()

        for (i, name) in enumerate(Name_Projections_85)
            Timeseries_Future = collect(Date(Timeseries_End[i,2]-29,1,1):Day(1):Date(Timeseries_End[i,2],12,31))
            #print(size(Timeseries_Past), size(Timeseries_Future))
            Timeseries_Proj = readdlm(path_85*name*"/"*Catchment_Name*"/pr_model_timeseries.txt")
            Timeseries_Proj = Date.(Timeseries_Proj, Dates.DateFormat("y,m,d"))
            Temperature = readdlm(path_85*name*"/"*Catchment_Name*"/tas_"*string(ID_temp)*"_sim1.txt", ',')[:,1]
            Elevation_Zone_Catchment, Temperature_Elevation_Catchment, Total_Elevationbands_Catchment = gettemperatureatelevation(Elevations_Catchment, Temperature)
            # get the temperature data at the mean elevation to calculate the mean potential evaporation
            Temperature = Temperature_Elevation_Catchment[:,findfirst(x-> x==Mean_Elevation_Catchment, Elevation_Zone_Catchment)]

            indexstart_past = findfirst(x-> x == Dates.year(Timeseries_Past[1]), Dates.year.(Timeseries_Proj))[1]
            indexend_past = findlast(x-> x == Dates.year(Timeseries_Past[end]), Dates.year.(Timeseries_Proj))[1]
            Temperature_Past = Temperature[indexstart_past:indexend_past] ./ 10
            #print(Dates.year(Timeseries_Future[end]), Dates.year.(Timeseries_Proj[end]))
            indexstart_future = findfirst(x-> x == Dates.year(Timeseries_Future[1]), Dates.year.(Timeseries_Proj))[1]
            indexend_future = findlast(x-> x == Dates.year(Timeseries_Future[end]), Dates.year.(Timeseries_Proj))[1]
            Temperature_Future = Temperature[indexstart_future:indexend_future] ./ 10
            # calculate monthly mean temperature
            Monthly_Temperature_Past, Month = monthly_discharge(Temperature_Past, Timeseries_Past)
            Monthly_Temperature_Future, Month_future = monthly_discharge(Temperature_Future, Timeseries_Future)

            #-------- PRECIPITATION ------------------
            Precipitation_All_Zones = Array{Float64, 1}[]
            Total_Precipitation_Proj = zeros(length(Timeseries_Proj))
            for j in 1: length(ID_Prec_Zones)
                    # get precipitation projections for the precipitation measurement
                    Precipitation_Zone = readdlm(path_85*name*"/"*Catchment_Name*"/pr_"*string(ID_Prec_Zones[j])*"_sim1.txt", ',')[:,1]
                    #print(size(Precipitation_Zone), typeof(Precipitation_Zone))
                    push!(Precipitation_All_Zones, Precipitation_Zone ./10)
                    Total_Precipitation_Proj += Precipitation_All_Zones[j].*Area_Zones_Percent[j]
            end
            #Total_Precipitation_Proj = Precipitation_All_Zones[1].*Area_Zones_Percent[1] + Precipitation_All_Zones[2].*Area_Zones_Percent[2] + Precipitation_All_Zones[3].*Area_Zones_Percent[3] + Precipitation_All_Zones[4].*Area_Zones_Percent[4]
            Precipitation_Past = Total_Precipitation_Proj[indexstart_past:indexend_past]
            Precipitation_Future = Total_Precipitation_Proj[indexstart_future:indexend_future]

            # Monthly_Precipitation_Past, Month = monthly_precipitation(Precipitation_Past, Timeseries_Past)
            # Monthly_Precipitation_Future, Month_future = monthly_precipitation(Precipitation_Future, Timeseries_Future)
            #
            # Monthly_Precipitation_Past, Month = monthly_precipitation(Precipitation_Past, Timeseries_Past)
            # Monthly_Precipitation_Future, Month_future = monthly_precipitation(Precipitation_Future, Timeseries_Future)

            # make ecdf plots
            if time == "past"
                plot!(StatsBase.ecdf(Precipitation_Past), label="Proj "*string(i), color=[Farben_proj[i]], size=(800,800), dpi=300, linewidth=2)
            elseif time == "future"
                plot!(StatsBase.ecdf(Precipitation_Future), label="Proj "*string(i), color=[Farben_proj[i]], size=(800,800), dpi=300, linewidth=2)
            end
        end
        xlabel!("Precipitation [mm/d]")
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[h])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[h])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[h])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[h])*"m)", titlefont = font(20))
        end
        push!(all_plots_85, plot!(framestyle = :box))
    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15), minorticks=true, minorgrid=true, minorgridlinewidth=2)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_precipitation_flow_all_years_"*time*"_45.png")
    plot(all_plots_85[1], all_plots_85[2], all_plots_85[3], all_plots_85[4], all_plots_85[5], all_plots_85[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15), minorticks=true, minorgrid=true, minorgridlinewidth=2)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/ECDF_precipitation_flow_all_years_"*time*"_85.png")
end

function plot_ecdf_monthly_discharge_all_catchments(All_Catchment_Names, Elevation, nr_runs,current_month)
    Farben_proj = palette(:tab20)
    all_plots = []
    #for RCP 8.5
    # highest_temp_proj = [12,12,12,12,12,12]
    # lowest_temp_proj = [1,1,1,6,1,1]
    # highest_prec_proj = [8,8,8,8,12,2]
    # lowest_prec_proj = [10,10,10,10,10,10]
    # for RCP 4.5
    highest_temp_proj = [12,6,6,12,4,6]
    lowest_temp_proj = [13,13,13,13,13,13]
    highest_prec_proj = [7,2,12,8,12,2]
    lowest_prec_proj = [10,10,10,10,10,10]
    for (i,Catchment_Name) in enumerate(All_Catchment_Names)
        monthly_changes_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_8.5.txt", ',')
        monthly_changes_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Monthly_Discharge/discharge_months_4.5.txt", ',')
        plot()
        for proj in 1:14
            first_index = nr_runs[i] * 12 * (proj-1) + 1
            last_index = nr_runs[i] *12* proj
            months_85 = monthly_changes_85[first_index:last_index,1]
            Monthly_Discharge_past_85 = monthly_changes_85[first_index:last_index,2]
            Monthly_Discharge_future_85  = monthly_changes_85[first_index:last_index,3]
            months_45 = monthly_changes_45[first_index:last_index,1]
            Monthly_Discharge_past_45 = monthly_changes_45[first_index:last_index,2]
            Monthly_Discharge_future_45  = monthly_changes_45[first_index:last_index,3]
            Monthly_Discharge_Change_45  = monthly_changes_45[first_index:last_index,4]
            Monthly_Discharge_past_45 = convertDischarge(Monthly_Discharge_past_45, Area_Catchments[i])
            Monthly_Discharge_past_85 = convertDischarge(Monthly_Discharge_past_85, Area_Catchments[i])
            Monthly_Discharge_future_45 = convertDischarge(Monthly_Discharge_future_45, Area_Catchments[i])
            Monthly_Discharge_future_85 = convertDischarge(Monthly_Discharge_future_85, Area_Catchments[i])
            if proj == highest_temp_proj[i]
                farbe = Farben_proj[7]
            elseif proj == lowest_temp_proj[i]
                farbe = Farben_proj[8]
            elseif proj == highest_prec_proj[i]
                farbe = Farben_proj[1]
            elseif proj == lowest_prec_proj[i]
                farbe = Farben_proj[2]
            else
                farbe = Farben_proj[16]
            end
            # plot ecdf of specific month
            plot!(StatsBase.ecdf(relative_error(Monthly_Discharge_future_45[findall(x-> x == current_month, months_45)], Monthly_Discharge_past_45[findall(x-> x == current_month, months_45)])*100), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
            #plot!(StatsBase.ecdf(relative_error(Monthly_Discharge_future_85[findall(x-> x == current_month, months_85)], Monthly_Discharge_past_85[findall(x-> x == current_month, months_85)])*100), label="Proj "*string(proj), color=[farbe], size=(800,800), dpi=300, linewidth=2)
        end
        xlabel!("[%]", yguidefontsize=20)
        if Catchment_Name == "Pitten"
            title!("Feistritztal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "IllSugadin"
            title!("Silbertal ("*string(Elevation[i])*"m)", titlefont = font(20))
        elseif Catchment_Name == "Palten"
            title!("Paltental ("*string(Elevation[i])*"m)", titlefont = font(20))
        else
            title!(Catchment_Name*" ("*string(Elevation[i])*"m)", titlefont = font(20))
        end
        push!(all_plots, plot!(framestyle = :box))
    end
    plot(all_plots[1], all_plots[2], all_plots[3], all_plots[4], all_plots[5], all_plots[6], layout= (3,2), legend = false, size=(1000,1400), left_margin = [10mm 0mm],right_margin = 20px, bottom_margin = 20px, yguidefontsize=15, xtickfont = font(15), ytickfont = font(15), xguidefontsize = 15, dpi=300, titlefont = font(15))#, minorticks=true, minorgrid=true)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/Comparison_Proj/monthly_discharge/ECDF_discharge_month"*string(current_month)*"_45_new.png")
end

#plot_ecdf_prec_temp_all_catchments(Catchment_Names, Catchment_Height, nr_runs, "future")
# for month in 1:12
#     plot_ecdf_monthly_discharge_all_catchments(Catchment_Names, Catchment_Height, nr_runs,month)
# end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
using TrussMorph: SectionProperties, MaterialProperties
using LinearAlgebra: norm, cholesky, I
using SparseArrays

function local_coord_rotation(u::Vector{Float64}, v::Vector{Float64})
    local L = norm(u - v)
    c = (v - u) ./ L
    R = zeros(Float64, 3, 3)
    R[1,1] = c[1]
    R[1,2] = c[2]
    R[2,1] = -c[2]
    R[2,2] = c[1]
    R[3,3] = 1
    return R
end

"""
local_stiffness_matrix: compute elemental FRAME stiffness matrix
in the local coordinates.

Return: 6 x 6 matrix K_le (local elemental)

Note: only support 2d right now.
"""
function local_stiffness_matrix(L::Float64, mp::MaterialProperties, sp::SectionProperties)
    E = mp.E
    G = mp.G
    A = sp.A
    I = sp.Iz
    fs = 1.0 # shear deformation const
    beta = (12*E*I*fs)/(G*A*L^2)

    K_beam = zeros(Float64, 4, 4)
    K_beam[1,2] = 6*L;
    K_beam[1,3] = -12;
    K_beam[1,4] = 6*L;
    K_beam[2,3] = -6*L;
    K_beam[2,4] = L^2*(2-beta);
    K_beam[3,4] = -6*L;
    K_beam = K_beam' + K_beam;

    K_beam += Diagonal([12, L^2*(4+beta), 12, L^2*(4+beta)])
    K_beam *= (E*I)/(L^3*(1+beta))

    K_le = zeros(Float64, 6, 6)
    K_le[1,1] = E * A / L
    K_le[1,4] = -E * A / L
    K_le[4,1] = -E * A / L
    K_le[4,4] = E * A / L
    K_le[2:3,2:3] = K_beam[1:2,1:2];
    K_le[2:3,5:6] = K_beam[1:2,3:4];
    K_le[5:6,2:3] = K_beam[3:4,1:2];
    K_le[5:6,5:6] = K_beam[3:4,3:4];
    return K_le
end

# TODO: can be optimized using symbolic substitution, only X is changing here
# NOTE: r might be set to be uniform, if we are using ∑ f_e * l_e formula?
function assemble_global_stiffness_matrix(X::Matrix{Float64}, T::Matrix{Int64},
                                          r::Vector{Float64}, mp::MaterialProperties,
                                          node_dof::Int, full_node_dof::Int)
    n_elements = size(T,1)
    n_nodes = size(X,1)

    # element -> dof id map
    sys_dof::Int = node_dof * n_nodes
    id_map = zeros(Int, n_elements, 2*node_dof)
    back_dof_lin = collect(node_dof-1:-1:0)
    for i=1:n_elements
        uid = T[i,1]
        vid = T[i,2]
        id_map[i, 1:node_dof] = uid * node_dof * ones(node_dof) - back_dof_lin
        id_map[i, node_dof+1:2*node_dof] = vid * node_dof * ones(node_dof) - back_dof_lin
    end

    # cross sec properties
    sp_vec = Array{SectionProperties}(undef, n_elements)
    for e=1:n_elements
        sp_vec[e] = compute_round_section_properties(r[e])
    end

    # for elemental force calc
    KR_es = Array{Matrix{Float64},1}(undef, n_elements)

    if 2 == node_dof
        ex_id = [1, 4]
        xy_id = [1, 2, 4, 5]
    else
        ex_id = collect(1:2*node_dof)
        xy_id = collect(1:2*node_dof)
    end

    # TODO: array size node_dof^2 * n_elements
    I = Int[]
    J = Int[]
    V = Float64[]
    # K ∈ (node_dof*n_node)^2
    for e=1:n_elements
        u = X[T[e,1],:]
        v = X[T[e,2],:]
        L = norm(u - v)

        # element local stiffness matrix
        K_le = local_stiffness_matrix(L, mp, sp_vec[e])

        # local -> global rotation matrix
        R_b = local_coord_rotation(u, v)
        R = zeros(Float64, 2*full_node_dof, 2*full_node_dof)
        for k=1:(Int).(full_node_dof/3) * 2
            R[k*3-3+1:k*3, k*3-3+1:k*3] = R_b
        end

        K_Ge = R[ex_id, xy_id]' * K_le[ex_id, ex_id] * R[ex_id, xy_id]
        KR_es[e] = K_le[ex_id, ex_id] * R[ex_id, xy_id]

        for i=1:2*node_dof
            for j=1:2*node_dof
                push!(I, id_map[e,i])
                push!(J, id_map[e,j])
                push!(V, K_Ge[i,j])
            end
        end
    end
    return sparse(I, J, V, sys_dof, sys_dof), KR_es, id_map
end

function assemble_load_vector(Load::Matrix{Float64}, n_nodes::Int, node_dof::Int)
    @assert(size(Load,2) == node_dof + 1)
    sys_dof::Int = node_dof * n_nodes
    n_load_nodes = size(Load, 1)
    I = Int[]
    V = Float64[]
    for i=1:n_load_nodes
        v_id = Int.(Load[i,1])
        append!(I, node_dof*(v_id-1)+1 : node_dof*v_id)
        append!(V, Load[i, 2:1:node_dof+1])
    end
    return sparsevec(I, V, sys_dof)
end

function dof_permutation(S::Matrix{Int}, n_nodes::Int, node_dof::Int)
    sys_dof::Int = node_dof * n_nodes
    n_fix = size(S,1)
    n_dof_fix = sum(sum(S[:,2:1+node_dof] .== 1))

    full_dof_fix = zeros(Bool, sys_dof)
    for i=1:n_fix
        v_id = Int.(S[i,1])
        full_dof_fix[node_dof*(v_id-1)+1 : node_dof*v_id] = S[i, 2:1+node_dof]
        #TODO: 3D node_dof here might be incompatible, 2D is okay
    end

    perm = zeros(Int, sys_dof)
    n_dof_free = sys_dof - n_dof_fix
    free_tail = 1
    fix_tail = n_dof_free + 1
    for i=1:sys_dof
        if !full_dof_fix[i]
            perm[free_tail] = i
            free_tail += 1
        else
            perm[fix_tail] = i
            fix_tail += 1
        end
    end
    I = collect(1:sys_dof)
    perm_spm = sparse(I, perm, ones(sys_dof), sys_dof, sys_dof)
    return perm, perm_spm, n_dof_free
end

function get_weight_calculation_fn(X_full::Matrix{Float64}, r::Vector{Float64}, T::Matrix{Int64},
    F_perm_m::Vector{Float64}, perm::SparseMatrixCSC{Float64, Int}, n_dof_free::Int, node_dof::Int, full_node_dof::Int, mp::MaterialProperties, design_var_id::Vector{Int}; opt_compliance::Bool=false, kwargs...)::Function

    # assemble stiffness matrix
    n_v = size(X_full, 1)
    n_e = size(T, 1)
    sys_dof = n_v * node_dof

    # split design variable + the rest
    X_template = copy(X_full)
    design_var_id_map = zeros(Int, length(design_var_id), 2)

    function id_convert(id::Int)::Vector{Int}
        return [Int.((id + id%2) / 2), 2 - id%2]
    end

    for i=1:length(design_var_id)
        design_var_id_map[i,:] = id_convert(design_var_id[i])
    end

    # sym_var_ids = undef
    # if :sym_var_ids in keys(kwargs) && kwargs[:sym_var_ids] != undef
    #     sym_var_ids = kwargs[:sym_var_ids]
    #     @assert(size(sym_var_ids, 2) == 2)
    # end

    function calc_weight(X_var::Vector{Float64})
        for i=1:length(design_var_id)
            X_template[id_convert(design_var_id[i])[1], id_convert(design_var_id[i])[2]] = X_var[i]
        end
        # if sym_var_ids != undef
        #     for i=1:size(sym_var_ids,1)
        #         X_template[id_convert(sym_var_ids[i,2])[1], id_convert(sym_var_ids[i,2])[2]] = X_template[id_convert(sym_var_ids[i,1])[1], id_convert(sym_var_ids[i,1])[2]]
        #     end
        # end

        K, KR_es, id_map = assemble_global_stiffness_matrix(X_template, T, r, mp, node_dof, full_node_dof)
        K_perm = perm * K * perm'
        K_mm = K_perm[1:n_dof_free, 1:n_dof_free]

        # solve linear system
        # U_m = K_mm\F_perm_m
        K_mm = (K_mm + K_mm')/2
        Kpf = cholesky(K_mm, check=false)
        if !LinearAlgebra.issuccess(Kpf)
            return 1e10
        end
        U_m = Kpf\F_perm_m

        # U_I, U_V = findnz(U_m)
        U_perm = vcat(U_m, zeros(sys_dof - n_dof_free))

        # compliance
        if opt_compliance
            return F_perm_m' * U_m
        end

        U = perm' * U_perm

        # TODO: buckling sizing is ignored for now
        # calculate ∑ F_e * l_e
        weight = 0.0
        if node_dof == 2
            e_react_dof = 1 # truss
        else
            e_react_dof = 3 # frame
        end

        eF = zeros(n_e, e_react_dof*2)
        for e=1:n_e
            eF[e,:] = KR_es[e] * U[id_map[e,:]]
            eL = norm(X_template[T[e,1], :] - X_template[T[e,2], :])
            if 2 == node_dof
                weight += abs(eF[e,:][1]) * eL
            else
                # axial force only
                # weight += abs(eF[e,:][1]) * eL

                # moment only
                # weight += 4*max(abs(eF[e,:][3]), abs(eF[e,:][6]))/r[e]) * eL

                # 4M + F
                weight += (abs(eF[e,:][1]) + 4*max(abs(eF[e,:][3]), abs(eF[e,:][6]))) * eL
                # weight += (abs(eF[e,:][1])/(pi * r[e]^2) + abs(eF[e,:][3])/(pi * r[e]^3 / 4)) * eL
            end
        end
        # @show eF
        # @show U
        return weight
    end

    return calc_weight
end

ext: jl


text:
module StringBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools
using Compat

import Compat: UTF8String, view

const SUITE = BenchmarkGroup()

####################
# replace (#12224) #
####################

str = join(samerand('a':'d', 10^4))

SUITE["replace"] = @benchmarkable replace($str, "a", "b")
SUITE["join"] = @benchmarkable join($str, $str) time_tolerance=0.40

str = "Gf6FPPWevqer3di13haDSzrRrSiThqmV3k02dALLu7OHdYRR0dfrKf4iCMcDvgZBawx"

g = addgroup!(SUITE, "search")
g["Char"] = @benchmarkable search($str,  $('x'))
g["String"] = @benchmarkable search($str, $("x"))

g = addgroup!(SUITE, "searchindex")
g["Char"] = @benchmarkable searchindex($str,  $('x'))
g["String"] = @benchmarkable searchindex($str, $("x"))

end # module


text:
module BioFetch

using FASTX
using GenomicAnnotations
using BioServices.EUtils
using BioServices.EBIProteins

export fetchseq, SeqFormat, fasta, gb

@enum SeqFormat fasta gb

"""
    fetchseq(ids::AbstractString::AbstractVector{<:AbstractString}, range::Union{Nothing, AbstractRange} = nothing, revstrand = false; format::SeqFormat = fasta)

Fetches sequence data from a database by accession number in either FASTA format or GenBank flatfile format.
Nucleotide and protein records may be mixed. Results will be returned in the order provided.
Supports NCBI, Ensembl, and UniProt accession numbers.

GenBank format may not work right now.

```julia
fetchseq("AH002844")                # retrive one FASTA NCBI nucleotide record
fetchseq(["CAA41295.1", "NP_000176"]) # retrieve two FASTA NCBI protein records
fetchseq("NC_036893.1", 81_775_230 .+ (1:1_000_000))         # retrive a 1 Mb segment of a FASTA NCBI genomic record
fetchseq("NC_036893.1", 81000000:81999999, true) # retrive a 1 Mb segment of a FASTA NCBI genomic record on the reverse strand
```

`range` and `revstrand` are only valid for nucleotide queries
"""
function fetchseq(ids::AbstractVector{<:AbstractString}, range::Union{Nothing, AbstractRange} = nothing, revstrand = false; format::SeqFormat = fasta)
    ncbinucleotides = String[]
    ncbiproteins = String[]
    ebiuniprots = String[]
    ebiensembls = String[]
    order = IdDict(ncbinucleotides => [], ncbiproteins => [], ebiuniprots => [], ebiensembls => [])
    for (i, id) ∈ enumerate(ids)
        ebiensembl = startswith(id, r"ENS[A-Z][0-9]{11}")
        ncbiprotein = startswith(id, r"[NX]P_|[A-Z]{3}[0-9]")
        ncbinucleotide = startswith(id, r"[NX][CGRTMW]_|[A-Z]{2}[0-9]|[A-Z]{4,6}[0-9]") && !ebiensembl
        ebiuniprot = startswith(id, r"[A-Z][0-9][A-Z0-9]{4}")
        ncbiprotein || ncbinucleotide || ebiuniprot || ebiensembl || throw("could not infer database for $id")

        array = ncbinucleotide ? ncbinucleotides :
                ncbiprotein ? ncbiproteins :
                ebiuniprot ? ebiuniprots :
                ebiensembl ? ebiensembls :
                error()
        push!(array, id)
        push!(order[array], i)
    end

    results = [fetchseq_ncbi(ncbinucleotides, "nuccore"; format, range, revstrand);
               fetchseq_ncbi(ncbiproteins, "protein"; format);
               fetchseq_uniprot(ebiuniprots; format);
               fetchseq_ensembl(ebiensembls; format)]
    order = [order[ncbinucleotides]; order[ncbiproteins]; order[ebiuniprots]; order[ebiensembls]]

    return results[order]
end

function fetchseq(id::AbstractString, range::Union{Nothing, AbstractRange} = nothing, revstrand = false; format::SeqFormat = fasta)
    ebiensembl = startswith(id, r"ENS[A-Z][0-9]{11}")
    ncbiprotein = startswith(id, r"[NX]P_|[A-Z]{3}[0-9]")
    ncbinucleotide = startswith(id, r"[NX][CGRTMW]_|[A-Z]{2}[0-9]|[A-Z]{4,6}[0-9]") && !ebiensembl
    ebiuniprot = startswith(id, r"[A-Z][0-9][A-Z0-9]{4}")

    result = ncbinucleotide ? fetchseq_ncbi(id, "nuccore"; format, range, revstrand) :
           ncbiprotein ? fetchseq_ncbi(id, "protein"; format) :
           ebiuniprot ? fetchseq_uniprot([id]; format) :
           ebiensembl ? fetchseq_ensembl([id]; format) :
           error("could not infer database for $id")

    return result
end

function fetchseq_ncbi(ids, db::AbstractString; format::SeqFormat = fasta, range::Union{Nothing, AbstractRange} = nothing, revstrand = false)
    isempty(ids) && return []
    if isnothing(range)
        response = efetch(; db, id = ids, rettype = String(Symbol(format)), retmode="text", strand = revstrand ? 2 : 1)
    else
        response = efetch(; db, id = ids, rettype = String(Symbol(format)), retmode="text", seq_start = first(range), seq_stop = last(range), strand = revstrand ? 2 : 1)
    end
    body = IOBuffer(response.body)
    reader = format == fasta ? FASTA.Reader(body) :
             format == gb ? GenBank.Reader(body) :
             nothing
    records = [record for record ∈ reader]
    return records
end

function fetchseq_uniprot(ids::AbstractVector; format::SeqFormat = fasta)
    isempty(ids) && return []
    records = []
    contenttype = format == fasta ? "text/x-fasta" : format == gb ? "text/flatfile" : error("unknown format")
    for id ∈ ids
        response = ebiproteins(; accession = id, contenttype)
        body = IOBuffer(response.body)
        reader = format == fasta ? FASTA.Reader(body) :
                 format == gb ? GenBank.Reader(body) :
                 nothing
        records = [record for record ∈ reader]
        append!(records, records)
    end
    return records
end

function fetchseq_ensembl(ids::AbstractVector; format::SeqFormat = fasta)
    isempty(ids) && return []
    records = []
    contenttype = format == fasta ? "text/x-fasta" : format == gb ? "text/flatfile" : error("unknown format")
    for id ∈ ids
        response = ebiproteins(; dbtype = "Ensembl", dbid = id, contenttype)
        body = IOBuffer(response.body)
        reader = format == fasta ? FASTA.Reader(body) :
                 format == gb ? GenBank.Reader(body) :
                 nothing
        records = [record for record ∈ reader]
        append!(records, records)
    end
    return records
end

end


text:
module CYCLOPS_v6_2a_Seed
export clean_data!
export getseed
export dispersion
export dispersion!
export getseed_mca
#export setdiff
export getseed_homologuesymbol
export getseed_homologueprobe
export getseed_homologuesymbol_brain

function clean_data!(data::Array{Float64,2},bluntpercent=0.99)
	ngenes,nsamples=size(data)
	nfloor=Int(1+floor((1-bluntpercent)*nsamples))
	nceiling=Int(ceil(bluntpercent*nsamples))
	for row=1:ngenes
		sorted=sort(vec(data[row,:]))
		vfloor=sorted[nfloor]
		vceil=sorted[nceiling]
		for sample =1:nsamples
			data[row,sample]=max(vfloor,data[row,sample])
			data[row,sample]=min(vceil ,data[row,sample])
		end
	end
	data
end

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

function getseed(data::Array{Any,2},symbol_list,maxcv=.75,mincv=.07,minmean=500,blunt=.99)
	data_symbols=data[2:end ,2]
	data_data=data[2:end , 4:end];
	data_data=Array{Float64}(data_data)
	data_data=clean_data!(data_data,blunt)
	ngenes,namples=size(data_data)

	gene_means=mean(data_data,2)
	gene_sds=std(data_data,2)
	gene_cvs= gene_sds ./ gene_means

	criteria1=findin(data_symbols,symbol_list)
	criteria2=findin((gene_means .> minmean),true)
	criteria3=findin((gene_cvs .> mincv),true)
	criteria4=findin((gene_cvs .< maxcv),true)

	allcriteria=intersect(criteria1,criteria2,criteria3,criteria4)
	seed_data=data_data[allcriteria,:]
	seed_symbols=data_symbols[allcriteria,:]
	seed_symbols, seed_data
end

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

function getseed_mca(data::Array{Any,2},probe_list,maxcv=.75,mincv=.07,minmean=500,blunt=.99)

	data_symbols=data[2:end ,2];
	data_probes=data[2:end ,1];

	data_data=data[2:end , 4:end];
	data_data=Array{Float64}(data_data)
	data_data=clean_data!(data_data,blunt)
	ngenes,namples=size(data_data)

	gene_means=mean(data_data,2)
	gene_sds=std(data_data,2)
	gene_cvs= gene_sds ./ gene_means

	criteria1=findin(data_probes,probe_list)
	criteria2=findin((gene_means .> minmean),true)
	criteria3=findin((gene_cvs .> mincv),true)
	criteria4=findin((gene_cvs .< maxcv),true)

	allcriteria=intersect(criteria1,criteria2,criteria3,criteria4)
	seed_data=data_data[allcriteria,:]
	seed_symbols=data_symbols[allcriteria,:]
	seed_symbols, seed_data
end
###############

function getseed_homologuesymbol(data::Array{Any,2},symbol_list,maxcv=.75,mincv=.07,minmean=500,blunt=.99)

	data_symbols=data[2:end ,2];
	data_probes=data[2:end ,1];

	data_data=data[2:end , 3:end];
	data_data=Array{Float64}(data_data)
	data_data=clean_data!(data_data,blunt)
	ngenes,namples=size(data_data)

	gene_means=mean(data_data,2)
	gene_sds=std(data_data,2)
	gene_cvs= gene_sds ./ gene_means

	criteria1=findin(data_symbols,symbol_list)
	criteria2=findin((gene_means .> minmean),true)
	criteria3=findin((gene_cvs .> mincv),true)
	criteria4=findin((gene_cvs .< maxcv),true)

	allcriteria=intersect(criteria1,criteria2,criteria3,criteria4)
	seed_data=data_data[allcriteria,:]
	seed_symbols=data_symbols[allcriteria,:]
	seed_symbols, seed_data
end

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

function getseed_homologueprobe(data::Array{Any,2},probe_list,maxcv=.75,mincv=.07,minmean=500,blunt=.99)

	data_symbols=data[2:end ,2];
	data_probes=data[2:end ,1];

	data_data=data[2:end , 3:end];
	data_data=Array{Float64}(data_data)
	data_data=clean_data!(data_data,blunt)
	ngenes,namples=size(data_data)

	gene_means=mean(data_data,2)
	gene_sds=std(data_data,2)
	gene_cvs= gene_sds ./ gene_means

	criteria1=findin(data_probes,probe_list)
	criteria2=findin((gene_means .> minmean),true)
	criteria3=findin((gene_cvs .> mincv),true)
	criteria4=findin((gene_cvs .< maxcv),true)

	allcriteria=intersect(criteria1,criteria2,criteria3,criteria4)
	seed_data=data_data[allcriteria,:]
	seed_symbols=data_symbols[allcriteria,:]
	seed_symbols, seed_data
end

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

function dispersion!(data::Array{Float64,2})
	ngenes,nsamples=size(data)
	for gene=1:ngenes
		genemean=mean(data[gene,:])
		for sample=1:nsamples
			data[gene,sample]=(data[gene,sample]-genemean)/genemean
		end
	end
	data
end

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

function dispersion(data::Array{Float64,2})
	ngenes,nsamples=size(data)
	ndata=zeros(ngenes,nsamples)
	for gene=1:ngenes
		genemean=mean(data[gene,:])
		for sample=1:nsamples
			ndata[gene,sample]=(data[gene,sample]-genemean)/genemean
		end
	end
	ndata
end

##############################################
#function setdiff(x,y)
#       	z= [x[1]]
#       	for element in x
#       		if !in(element,y)
#       			push!(z,element)
#       		end
#       	end
#    z=z[2:length(z)]
#    z
#end

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

function getseed_homologuesymbol_brain(data::Array{Any,2},symbol_list,maxcv=.75,mincv=.07,minmean=500,blunt=.99)

	data_symbols=data[4:end ,2];
	data_probes=data[4:end ,1];

	data_data=data[4:end , 4:end];
	data_data=Array{Float64}(data_data)
	data_data=clean_data!(data_data,blunt)
	ngenes,namples=size(data_data)

	gene_means=mean(data_data,2)
	gene_sds=std(data_data,2)
	gene_cvs= gene_sds ./ gene_means

	criteria1=findin(data_symbols,symbol_list)
	criteria2=findin((gene_means .> minmean),true)
	criteria3=findin((gene_cvs .> mincv),true)
	criteria4=findin((gene_cvs .< maxcv),true)

	allcriteria=intersect(criteria1,criteria2,criteria3,criteria4)
	seed_data=data_data[allcriteria,:]
	seed_symbols=data_symbols[allcriteria,:]
	seed_symbols, seed_data
end

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

end

max_stars_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_issues_repo_name: AdithyaVenkateshMohan/CYCLOPS

max_forks_repo_name: AdithyaVenkateshMohan/CYCLOPS


text:
include("current_search.jl")
using SpikingNeuralNetworks
SNN = SpikingNeuralNetworks
using Evolutionary, Test, Random
using Distributed
using SharedArrays
using SignalAnalysis

function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::GA, options)
    idx = sortperm(state.fitpop)
    record["fitpop"] = state.fitpop[:]#idx[1:last(idx)]]
    record["pop"] = population[:]
end

###
# Code appropriated from:
# https://github.com/JuliaML/MLPlots.jl/blob/master/src/optional/onlineai.jl
###

function SpikeTrains(n::Integer; kw...)
    plt = plot(n; leg=false, yticks=nothing, kw...)
    SpikeTrains(n, plt)
end

const _halfheight = 0.6
function Base.push!(spiketrains::SpikeTrains, idx::Integer, t::Real)
    append!(spiketrains.plt, idx, Float64[NaN, t, t],
            Float64[NaN, idx + _halfheight, idx - _halfheight])
    spiketrains
end
function Base.push!{T<:Real}(spiketrains::SpikeTrains, ts::AbstractVector{T})
    for i in 1:length(ts)
        push!(spiketrains, i, ts[i])
    end
    spiketrains
end

function get_ranges(ranges)

    lower = Float32[]
    upper = Float32[]
    for (k,v) in ranges
        append!(lower,v[1])
        append!(upper,v[2])
    end
    lower,upper
end

function init_b(lower,upper)
    gene = []#Float32
    for (i,(l,u)) in enumerate(zip(lower,upper))
        p1 = rand(l:u, 1)
        append!(gene,p1)
    end
    gene
end

function initf(n)
    genesb = []#Float32[]
    for i in 1:n
        genes = init_b(lower,upper)
        append!(genesb,[genes])
    end
    genesb
end

function custom_raster(P1::Array,P2::Array)
    y0 = Int32[0]
    X = Float32[];
	Y = Float32[]
	N=1
	y=0
	for p in [P1,P2]
		y+=1
		append!(X, P1)
		append!(Y, 1 .+ sum(y0))
    	push!(y0, 1)
		append!(X, P2)
		append!(Y, 1 .+ sum(y0))
	    push!(y0, 1)
		N+=1
    end
    plt = scatter(X, Y,w=50, m = (1, :black), leg = :none, marker=:vline,
                  xaxis=("t", (0, Inf)), yaxis = ("neuron",))
    y0 = y0[2:end-1]
    !isempty(y0) && hline!(plt, cumsum(y0), linecolor = :red)
    return plt
end

function custom_raster2(P1::Array,P2::Array)
	#neuron1=1*zeros(length(P1))
	#neuron2 = [2*i for i in neuron1]
	#@show(P1)
	#@show(P2)
    plt3=plot!(P1, seriestype="vline", label="")
	plt4=plot!(P2, seriestype="vline", label="")
	display(plt3)
	display(plt4)

	#plt3 = vline!(P1)#,P1)#,w=50)
	#display(plt3)
	#plt4 = vline!(P2)
	#display(plt4)
	#plot!(plt3,plt3)#,P2)#,w=50)
    #return plt1
end

function get_vm(p,simulation_duration)
    vm = p.records[:v]
	vm = [i[1] for i in vm]
	vm = signal(vm, length(vm)/simulation_duration)
	vm
end

function checkmodel(param,cell_type,ngt_spikes)
    if cell_type=="IZHI"
        pp = SNN.IZParameter(;a = param[1], b = param[2], c = param[3], d = param[4])
        E = SNN.IZ(;N = 1, param = pp)
    end
    if cell_type=="ADEXP"
		adparam = SNN.ADEXParameter(;a = param[1],
            b = param[2],
            cm = param[3],
            v0 = param[4],
            τ_m = param[5],
            τ_w = param[6],
            θ = param[7],
            delta_T = param[8],
            v_reset = param[9],
            spike_delta = param[10])
			E = SNN.AD(;N = 1, param=adparam)
    end

    ALLEN_DURATION = 2000 * ms
    ALLEN_DELAY = 1000 * ms
	current = current_search(cell_type,param,ngt_spikes)
    E.I = [current*nA]#_search(param,ngt_spikes)*nA]

    SNN.monitor(E, [:v])
	SNN.monitor(E, [:fire])

	###
	# TODO buid t into neuron models
	# SNN.monitor(E, [:t])
    ###
	simulation_duration= ALLEN_DURATION+ALLEN_DELAY+343ms
    SNN.sim!([E]; dt =1*ms, delay=ALLEN_DELAY,stimulus_duration=ALLEN_DURATION,simulation_duration = simulation_duration)
	spikes = get_spikes(E)
    spikes = [s/1000.0 for s in spikes]
	vm = get_vm(E,simulation_duration)
	#vm = E.records[:v]
	#vm = [i[1] for i in vm]

	return (vm,spikes)
end
function Evolutionary.value!(::Val{:multiproc}, fitness, objfun, population::AbstractVector{IT}) where {IT}
	fitness = SharedArrays.SharedArray{Float32}(fitness)
	@time @sync @distributed for i in 1:length(population)
		fitness[i] = value(objfun, population[i])
	end
	#@show(fitness)
	fitness
end
function Evolutionary.value!(::Val{:serial}, fitness, objfun, population::AbstractVector{IT}) where {IT}
	Threads.@threads for i in 1:length(population)
		fitness[i] = value(objfun, population[i])
	end
end
function get_data()
	#print(pwd())
	file="../JLD/ground_truth.jld"
    if false#isfile(file)
        vmgtv = load(file,"vmgtv")
        ngt_spikes = load(file,"ngt_spikes")
        gt_spikes = load(file,"gt_spikes")

        ground_spikes = gt_spikes
        ngt_spikes = size(gt_spikes)[1]
        vmgtt = load(file,"vmgtt")

        plot(plot(vmgtv,vmgtt,w=1))

    else
        py"""
        from neo import AnalogSignal
        from neuronunit.allenapi import make_allen_tests_from_id

        """

        py"""
        specimen_id = (
            325479788,
            324257146,
            476053392,
            623893177,
            623960880,
            482493761,
            471819401
        )
        specimen_id = specimen_id[1]
        target_num_spikes=7
        sweep_numbers, data_set, sweeps = make_allen_tests_from_id.allen_id_to_sweeps(specimen_id)
        (vmm,stimulus,sn,spike_times) = make_allen_tests_from_id.get_model_parts_sweep_from_spk_cnt(
            target_num_spikes, data_set, sweep_numbers, specimen_id
        )
        """
        gt_spikes = py"spike_times"
        ground_spikes = gt_spikes

        ngt_spikes = size(gt_spikes)[1]
        vmgtv = py"vmm.magnitude"
        vmgtt = py"vmm.times"

		#s_a = signal(opt_vec, length(opt_vec)/3.5)
		julia_version_vm = signal(vmgtv, length(vmgtt)/last(vmgtt))

		plot(plot(vmgtv,vmgtt,w=1))

        save("ground_truth.jld", "vmgtv", vmgtv,"vmgtt",vmgtt, "ngt_spikes", ngt_spikes,"gt_spikes",gt_spikes,"julia_version_vm",julia_version_vm)
        #filename = string("ground_truth: ", py"target_num_spikes")#,py"specimen_id)
        #filename = string(filename,py"specimen_id")
        #filename = string(filename,".jld")
        #save(filename, "vmgtv", vmgtv,"vmgtt",vmgtt, "ngt_spikes", ngt_spikes,"gt_spikes",gt_spikes)

    end

        vmgtv = load("ground_truth.jld","vmgtv")
        ngt_spikes = load("ground_truth.jld","ngt_spikes")
        ngt_spikes = size(gt_spikes)[1]
        ground_spikes = load("ground_truth.jld","gt_spikes")

        vmgtt = load("ground_truth.jld","vmgtt")
		julia_version_vm = load("ground_truth.jld","julia_version_vm")
    return (vmgtv,vmgtt,ngt_spikes,ground_spikes,julia_version_vm)
end
function get_izhi_ranges()
    ranges_izhi = DataStructures.OrderedDict{Char,Float32}()
    ranges_izhi = ("a"=>(0.002,0.3),"b"=>(0.02,0.36),"c"=>(-75,-35),"d"=>(0.005,16))#,"I"=>[100,9000])
    lower,upper = get_ranges(ranges_izhi)
    return lower,upper
end

function initd()
    population = initf(50)
    garray = zeros((length(population)[1], length(population[1])))
    for (i,p) in enumerate(population)
        garray[i,:] = p
    end
    garray[1,:]
end
function get_adexp_ranges()
    ranges_adexp = DataStructures.OrderedDict{String,Tuple{Float32,Float32}}()
    ranges_adexp[:"a"] = (2.0, 10)
    ranges_adexp[:"b"] = (5.0, 10)
    ranges_adexp[:"cm"] = (700.0, 983.5)
    ranges_adexp[:"v0"] = (-70, -55)
    ranges_adexp[:"τ_m"] = (10.0, 42.78345)
    ranges_adexp[:"τ_w"] = (300.0, 454.0)  # Tau_w 0, means very low adaption
    ranges_adexp[:"θ"] = (-45.0,-10)
    ranges_adexp[:"delta_T"] = (1.0, 5.0)
    ranges_adexp[:"v_reset"] = (-70.0, -15.0)
    ranges_adexp[:"spike_delta"] = (1.25, 20.0)
    lower,upper = get_ranges(ranges_adexp)
    return lower,upper
end

function vecplot(p, sym)
    v = SNN.getrecord(p, sym)
    y = hcat(v...)'
    x = 1:length(v)
    plot(x, y, leg = :none,
    xaxis=("t", extrema(x)),
    yaxis=(string(sym), extrema(y)))
end

function vecplot(P::Array, sym)
    plts = [vecplot(p, sym) for p in P]
    N = length(plts)
    plot(plts..., size = (600, 400N), layout = (N, 1))
end