text:
using DelimitedFiles
using Plots
using Statistics
using StatsPlots
using Plots.PlotMeasures
using CSV
using Dates
using DocStringExtensions

include("ObjectiveFunctions.jl")

relative_error(future, initial) = (future - initial) ./ initial
# ---------------------------------- ANNUAL MAXIMUM DISCHARGE --------------------------------
function circleShape(h,k,r)
    tau = LinRange(0, 2*pi, 500)
    h .+ r*sin.(tau), k .+ r*cos.(tau)
end

"""
Plots the maximum annual discharge of the past as a polar plot

$(SIGNATURES)

Needs Pyplot (using PyPlot)
"""
function AMF_circular_plot(Timing::Array{Float64,1}, Magnitude::Array{Float64,1}, Timeseries, Catchment_Name, Area_Catchment)
    @assert length(Timing) == length(Magnitude)
    Years = collect(Dates.year(Timeseries[1]): Dates.year(Timeseries[end]))
    Nr_Days_Year = Dates.daysinyear.(Years)
    # in m³/s
    fig = figure(figsize=(10,10)) # Create a new figure
    ax = PyPlot.axes(polar="true") # Create a polar axis
    PyPlot.title("Annual Maximum Discharge [m³/s] "*Catchment_Name*" Past")
    width = 2pi/365
    Farbe = "blue"
    for (i, current_timing) in enumerate(Timing)
        theta = current_timing * 2 * pi / Nr_Days_Year[i]
        b = bar(theta,Magnitude[i],width=width, color=Farbe)
    end # Bar plot

    dtheta = 10
    #Days_in_Month = [31,28,31,30,31,30,31,31,30,31,30,31]

    Days_in_Month = [0,31,59,90,120,151,181,212,243,273,304,334]
    Days_in_Month = Days_in_Month ./ 365 .* 360
    #ax.set_thetagrids(collect(0:dtheta:360-dtheta)) # Show grid lines from 0 to 360 in increments of dtheta
    ax.set_thetagrids(Days_in_Month)

    ax.set_xticklabels(["Jan", "Feb", "Mar", "Apr", "May","Jun","Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
    ax.set_theta_zero_location("N") # Set 0 degrees to the top of the plot
    ax.set_theta_direction(-1) # Switch to clockwise
    fig.canvas.draw() # Update the figure
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/Circular/"*Catchment_Name*"_AMF_Past.png")

    # in mm/d
    Magnitude = convertDischarge(Magnitude, Area_Catchment)

    fig = figure(figsize=(10,10)) # Create a new figure
    ax = PyPlot.axes(polar="true") # Create a polar axis
    PyPlot.title("Annual Maximum Discharge [mm/d] "*Catchment_Name*" Past")
    width = 2pi/365
    Farbe = "blue"
    for (i, current_timing) in enumerate(Timing)
        theta = current_timing * 2 * pi / Nr_Days_Year[i]
        b = bar(theta,Magnitude[i],width=width, color=Farbe)
    end # Bar plot

        dtheta = 10
        #Days_in_Month = [31,28,31,30,31,30,31,31,30,31,30,31]

        Days_in_Month = [0,31,59,90,120,151,181,212,243,273,304,334]
        Days_in_Month = Days_in_Month ./ 365 .* 360
        #ax.set_thetagrids(collect(0:dtheta:360-dtheta)) # Show grid lines from 0 to 360 in increments of dtheta
        ax.set_thetagrids(Days_in_Month)

        ax.set_xticklabels(["Jan", "Feb", "Mar", "Apr", "May","Jun","Jul", "Aug", "Sep", "Oct", "Nov", "Dec"])
        ax.set_theta_zero_location("N") # Set 0 degrees to the top of the plot
        ax.set_theta_direction(-1) # Switch to clockwise
        fig.canvas.draw() # Update the figure
        savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/Circular/"*Catchment_Name*"_AMF_Past_mm.png")
end

"""
Computes the date and magnitude of the yearly maximum discharge using calender year.

$(SIGNATURES)

The function returns the maginutde and Date as day of year of the maximum annual discharge.
    As input a discharge series and the corresponding timeseries is needed.
"""
function max_Annual_Discharge(Discharge, Timeseries)
    Years = collect(Dates.year(Timeseries[1]): Dates.year(Timeseries[end]))
    max_Annual_Discharge = Float64[]
    Date_max_Annual_Discharge = Float64[]
    for (i, Current_Year) in enumerate(Years)
            Dates_Current_Year = filter(Timeseries) do x
                              Dates.Year(x) == Dates.Year(Current_Year)
                          end
            max_Discharge = maximum(Discharge[indexin(Dates_Current_Year, Timeseries)])

            #Date_Max_Discharge = Timeseries[findfirst(x->x == max_Discharge, Discharge)]
            append!(max_Annual_Discharge, max_Discharge)
            append!(Date_max_Annual_Discharge, Dates.dayofyear(Dates_Current_Year[argmax(Discharge[indexin(Dates_Current_Year, Timeseries)])]))
    end
    return max_Annual_Discharge, Date_max_Annual_Discharge
end

"""
Computes the magnitude and timing of the yearly maximum discharge on 7 consecutive days using calender year (uses mean discharge of 7 days)

$(SIGNATURES)

The function returns the magnitude and Date as day of year of the first of 7 days of maximum annual discharge.
    As input a discharge series and the corresponding timeseries is needed.
"""
function max_Annual_Discharge_7days(Discharge, Timeseries)
    days = 7
    Years = collect(Dates.year(Timeseries[1]): Dates.year(Timeseries[end]))
    max_Annual_Discharge = Float64[]
    Date_max_Annual_Discharge = Float64[]
    for (i, Current_Year) in enumerate(Years)
            Dates_Current_Year = filter(Timeseries) do x
                              Dates.Year(x) == Dates.Year(Current_Year)
                          end
            Current_Discharge = Discharge[indexin(Dates_Current_Year, Timeseries)]

            All_Discharges_7days = Float64[]
            for week in 1: length(Current_Discharge) - days
                Current_Discharge_7days = mean(Current_Discharge[week: week+days])
                append!(All_Discharges_7days, Current_Discharge_7days)
            end
            #Date_Max_Discharge = Timeseries[findfirst(x->x == max_Discharge, Discharge)]
            append!(max_Annual_Discharge, maximum(All_Discharges_7days))
            append!(Date_max_Annual_Discharge, Dates.dayofyear(Dates_Current_Year[argmax(All_Discharges_7days)]))
    end
    return max_Annual_Discharge, Date_max_Annual_Discharge
end

"""
Computes the average timing of maximum annual discharges of a multiyear timeseries using circular statistics.

$(SIGNATURES)

The function returns the mean timing of the annual maximum discharge as well as the concentration of the date of occurrence around the average date.
    (Concentration = 1, all events occur on the same day, concentration = 0, events are widely spread)
"""
function average_timing(Dates_Max_Annual_Discharge, Timeseries)
    Years = collect(Dates.year(Timeseries[1]): Dates.year(Timeseries[end]))
    @assert length(Dates_Max_Annual_Discharge) ==  length(Years)
    Days_In_Year = Float64[]
    Formula_x = Float64[]
    Formula_y = Float64[]
    for (i, current_year) in enumerate(Years)
        Current_Circular_Date = Dates_Max_Annual_Discharge[i] * 2 * pi / Dates.daysinyear(current_year)
        x = cos(Current_Circular_Date)
        y = sin(Current_Circular_Date)
        append!(Days_In_Year, Dates.daysinyear(current_year))
        append!(Formula_x, x)
        append!(Formula_y, y)
    end
    mean_x = mean(Formula_x)
    mean_y = mean(Formula_y)
    mean_DaysinYear = mean(Days_In_Year)

    if mean_x > 0 && mean_y >= 0
        Mean_Timing = atan(mean_y / mean_x) * mean_DaysinYear / (2*pi)
    elseif mean_x <= 0
        Mean_Timing = (atan(mean_y / mean_x) + pi) * mean_DaysinYear / (2*pi)
    else
        Mean_Timing = (atan(mean_y / mean_x) + 2*pi) * mean_DaysinYear / (2*pi)
    end
    concentration_floods = sqrt(mean_x^2 + mean_y^2)

    return Mean_Timing, concentration_floods
end

function average_timing_simulations(Dates_Max_Annual_Discharge)
    Days_In_Year = Float64[]
    Formula_x = Float64[]
    Formula_y = Float64[]
    mean_DaysinYear = 365.23
    for i in 1:length(Dates_Max_Annual_Discharge)
        Current_Circular_Date = Dates_Max_Annual_Discharge[i] * 2 * pi / mean_DaysinYear
        x = cos(Current_Circular_Date)
        y = sin(Current_Circular_Date)
        append!(Formula_x, x)
        append!(Formula_y, y)
    end
    mean_x = mean(Formula_x)
    mean_y = mean(Formula_y)
    if mean_x > 0 && mean_y >= 0
        Mean_Timing = atan(mean_y / mean_x) * mean_DaysinYear / (2*pi)
    elseif mean_x <= 0
        Mean_Timing = (atan(mean_y / mean_x) + pi) * mean_DaysinYear / (2*pi)
    else
        Mean_Timing = (atan(mean_y / mean_x) + 2*pi) * mean_DaysinYear / (2*pi)
    end
    concentration_floods = sqrt(mean_x^2 + mean_y^2)
    std = sqrt(-2*log(concentration_floods)) * (mean_DaysinYear/ (2*pi))

    return Mean_Timing, std
end

"""
Converts the timing to circular statistics and calculates the differences in days of occurence between past and future.

$(SIGNATURES)

The function returns and array with the shift in occurences. As input an array with timing of AMF in past and future is needed.
"""
function difference_timing(Timing_Past, Timing_Future)
    Timing_Past = Timing_Past .* (2*pi) ./ 365.23333
    Timing_Future = Timing_Future .* (2*pi) ./ 365.23333
    Difference = Timing_Future - Timing_Past
    for (i,current_Difference) in enumerate(Difference)
        if current_Difference > pi
            Difference[i] = 2*pi - current_Difference
        elseif current_Difference < -pi
            Difference[i] = -(2*pi + current_Difference)
        end
    end
    return Difference .* (365.23333/(2*pi))
end

"""
Computes the average timing and magnitude of AMF for all projections for past and future.

$(SIGNATURES)

The function returns the mean timing of the mean magnitude of past and future timeseries, the mean timing of past and future
    and the concentraton of timing in past and future.

"""
function change_max_Annual_Discharge(path_to_projections, Catchment_Name)
    Name_Projections_45 = readdir(path_to_projections)
    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",',')
    change_all_runs = Float64[]
    average_max_Discharge_past = Float64[]
    average_max_Discharge_future = Float64[]
    Timing_max_Discharge_past = Float64[]
    Timing_max_Discharge_future = Float64[]
    All_Concentration_past = Float64[]
    All_Concentration_future = Float64[]
    if path_to_projections[end-2:end-1] == "45"
        index = 1
        rcp = "45"
        print(rcp, " ", rcp)
    elseif path_to_projections[end-2:end-1] == "85"
        index = 2
        rcp="85"
        print(rcp, " ", rcp)
    end
    for (i, name) in enumerate(Name_Projections_45)
        Timeseries_Future = collect(Date(Timeseries_End[i,index]-29,1,1):Day(1):Date(Timeseries_End[i,index],12,31))
        Past_Discharge = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_discharge_past_2010.csv", ',')
        println(size(Past_Discharge)[1])
        Future_Discharge = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_discharge_future_2100.csv", ',')
        #change_all_runs = Float64[]
        for run in 1:size(Past_Discharge)[1]
            max_Discharge_past, Date_max_Discharge_past = max_Annual_Discharge(Past_Discharge[run,:], Timeseries_Past)
            max_Discharge_future, Date_max_Discharge_future = max_Annual_Discharge(Future_Discharge[run,:], Timeseries_Future)
            append!(average_max_Discharge_past, mean(max_Discharge_past))
            append!(average_max_Discharge_future, mean(max_Discharge_future))
            timing_average_max_Discharge_past, Concentration_past = average_timing(Date_max_Discharge_past, Timeseries_Past)
            timing_average_max_Discharge_future, Concentration_future = average_timing(Date_max_Discharge_future, Timeseries_Future)
            append!(Timing_max_Discharge_past, timing_average_max_Discharge_past)
            append!(Timing_max_Discharge_future, timing_average_max_Discharge_future)
            append!(All_Concentration_past, Concentration_past)
            append!(All_Concentration_future, Concentration_future)
        end
    end
    return average_max_Discharge_past, average_max_Discharge_future, Timing_max_Discharge_past, Timing_max_Discharge_future, All_Concentration_past, All_Concentration_future
end

"""
Computes the timing and magnitude of AMF for all projections for past and future for all years using a probability distribution.

$(SIGNATURES)

"""
function change_max_Annual_Discharge_Prob_Distribution(path_to_projections, Catchment_Name)
    Name_Projections_45 = readdir(path_to_projections)
    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",',')
    change_all_runs = Float64[]
    average_max_Discharge_past = Float64[]
    average_max_Discharge_future = Float64[]
    Exceedance_Probability = Float64[]
    Timing_max_Discharge_past = Float64[]
    Timing_max_Discharge_future = Float64[]
    All_Concentration_past = Float64[]
    All_Concentration_future = Float64[]
    Date_Past = Float64[]
    Date_Future = Float64[]
    if path_to_projections[end-2:end-1] == "45"
        index = 1
        rcp = "45"
        print(rcp, " ", rcp)
    elseif path_to_projections[end-2:end-1] == "85"
        index = 2
        rcp="85"
        print(rcp, " ", rcp)
    end
    for (i, name) in enumerate(Name_Projections_45)
        Timeseries_Future = collect(Date(Timeseries_End[i,index]-29,1,1):Day(1):Date(Timeseries_End[i,index],12,31))
        Past_Discharge = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_discharge_past_2010.csv", ',')
        Future_Discharge = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_discharge_future_2100.csv", ',')
        println(size(Past_Discharge)[1])
        for run in 1:size(Past_Discharge)[1]
            max_Discharge_past, Date_max_Discharge_past = max_Annual_Discharge(Past_Discharge[run,:], Timeseries_Past)
            max_Discharge_future, Date_max_Discharge_future = max_Annual_Discharge(Future_Discharge[run,:], Timeseries_Future)
            # don't take mean of thirty years but probability distirbution
            max_Discharge_past_sorted, Prob_Dis_past = flowdurationcurve(max_Discharge_past)
            max_Discharge_future_sorted, Prob_Dis_future = flowdurationcurve(max_Discharge_future)
            #@assert Prob_Dis_past == Prob_Dis_future
            append!(average_max_Discharge_past, max_Discharge_past_sorted)
            append!(average_max_Discharge_future, max_Discharge_future_sorted)
            append!(Exceedance_Probability, Prob_Dis_past)
            append!(Date_Past, Date_max_Discharge_past)
            append!(Date_Future, Date_max_Discharge_future)
        end
    end
    return average_max_Discharge_past, average_max_Discharge_future, Exceedance_Probability, Date_Past, Date_Future #,Timing_max_Discharge_past, Timing_max_Discharge_future, All_Concentration_past, All_Concentration_future
end

"""
Plots the absolute and relative difference in magnitude of mean annual maximum discharge.
Also plots the abosulte difference in timing of AMF

$(SIGNATURES)

"""
function plot_Max_Flows(Max_Flows_past45, Max_Flows_future45, Max_Flows_past85, Max_Flows_future85, Timing_Max_Flows_past45, Timing_Max_Flows_future45,  Timing_Max_Flows_past85, Timing_Max_Flows_future85, Catchment_Name, Area_Catchment, nr_runs)
    Farben45=palette(:blues)
    Farben85=palette(:reds)
    # plot flows of each projection
    for proj in 1:14
        boxplot(Max_Flows_past45[1+(proj-1)*nr_runs: proj*nr_runs], color=[Farben45[1]])
        boxplot!(Max_Flows_future45[1+(proj-1)*nr_runs: proj*nr_runs],color=[Farben45[2]])
        boxplot!(Max_Flows_past85[1+(proj-1)*nr_runs: proj*nr_runs], color=[Farben85[1]])
        boxplot!(Max_Flows_future85[1+(proj-1)*nr_runs: proj*nr_runs], size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
        xticks!([1:4;], ["Past 4.5", "Future 4.5", "Past 8.5", "Future 8.5"])
        ylabel!("mean annual maximum daily Discharge [m³/s]")
        ylims!((15,35))
        title!("Annual Maximum Discharge")
        savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_max_yearly_discharge_"*string(Name_Projections_45[proj])*".png")
    end
    # plot flows of all projections combined
    boxplot(Max_Flows_past45, notch=true, color=[Farben45[1]])
    boxplot!(Max_Flows_future45,notch=true, color=[Farben45[2]])
    boxplot!(Max_Flows_past85, notch=true, color=[Farben85[1]])
    boxplot!(Max_Flows_future85, notch=true, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:4;], ["Past", "Future 4.5", "Future 8.5"])
    ylabel!("Mean annual maximum yearly Discharge [m³/s]")
    #ylims!((40,100))
    title!("Annual Maximum Discharge")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_max_yearly_discharge_notch.png")

    #absolute and relative decrease
    boxplot(Max_Flows_future45 - Max_Flows_past45,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    boxplot!(Max_Flows_future85 - Max_Flows_past85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("absolute change [m³/s]")
    #title!("Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    absolute_change = boxplot!()
    # relative change
    boxplot(relative_error(Max_Flows_future45, Max_Flows_past45)*100,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    boxplot!(relative_error(Max_Flows_future85, Max_Flows_past85)*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("relative change [%]")
    #title!("Relative Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    relative_change = boxplot!()
    plot(absolute_change, relative_change)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_change_max_yearly_discharge.png")

    #------------ ABSOLUTE AND RELATIVE CHANGE VIOLIN PLOTS ------------------------
    violin(Max_Flows_future45 - Max_Flows_past45,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    violin!(Max_Flows_future85 - Max_Flows_past85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("absolute change [m³/s]")
    #title!("Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    absolute_change = boxplot!()
    # relative change
    violin(relative_error(Max_Flows_future45, Max_Flows_past45)*100,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    violin!(relative_error(Max_Flows_future85, Max_Flows_past85)*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("relative change [%]")
    #title!("Relative Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    relative_change = boxplot!()

    plot(absolute_change, relative_change)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_change_max_yearly_discharge_violin.png")

    #------------ mm/d -------------
    Max_Flows_past45 = convertDischarge(Max_Flows_past45, Area_Catchment)
    Max_Flows_future45 = convertDischarge(Max_Flows_future45, Area_Catchment)
    Max_Flows_past85 = convertDischarge(Max_Flows_past85, Area_Catchment)
    Max_Flows_future85 = convertDischarge(Max_Flows_future85, Area_Catchment)
    boxplot(Max_Flows_past45, notch=true, color=[Farben45[1]])
    boxplot!(Max_Flows_future45, notch=true, color=[Farben45[2]])
    boxplot!(Max_Flows_past85, notch=true, color=[Farben85[1]])
    boxplot!(Max_Flows_future85, notch=true, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:4;], ["Past", "Future 4.5", "Future 8.5"])
    ylabel!("Mean annual maximum yearly Discharge [mm/d]")
    #ylims!((40,100))
    title!("Annual Maximum Discharge")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_max_yearly_discharge_mm_notch.png")

    #absolute and relative decrease
    boxplot(Max_Flows_future45 - Max_Flows_past45,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    boxplot!(Max_Flows_future85 - Max_Flows_past85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("absolute change [mm/d]")
    #title!("Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    absolute_change = boxplot!()
    # relative change
    boxplot(relative_error(Max_Flows_future45, Max_Flows_past45)*100,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    boxplot!(relative_error(Max_Flows_future85, Max_Flows_past85)*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("relative change [%]")
    #title!("Relative Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    relative_change = boxplot!()
    plot(absolute_change, relative_change)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_change_max_yearly_discharge_mm.png")

    #------------ ABSOLUTE AND RELATIVE CHANGE VIOLIN PLOTS ------------------------
    violin(Max_Flows_future45 - Max_Flows_past45,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    violin!(Max_Flows_future85 - Max_Flows_past85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("absolute change [mm/d]")
    #title!("Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    absolute_change = boxplot!()
    # relative change
    violin(relative_error(Max_Flows_future45, Max_Flows_past45)*100,color=[Farben45[2]])
    #boxplot!(, color=[Farben85[1]])
    violin!(relative_error(Max_Flows_future85, Max_Flows_past85)*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("relative change [%]")
    #title!("Relative Change in Summer Low Flows 30 year average of Lowest 7 day runoff from May - Nov (Future - Present)")
    relative_change = boxplot!()

    plot(absolute_change, relative_change)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_change_max_yearly_discharge_violin_mm.png")

    # ----------------- TIMING -----------------
    boxplot([1], Timing_Max_Flows_past45, color=[Farben45[1]])
    boxplot!([2], Timing_Max_Flows_past85, color=[Farben45[1]])
    boxplot!([3], Timing_Max_Flows_future45,color=[Farben45[2]])
    #boxplot!(Timing_Max_Flows_past85, color=[Farben85[1]])
    boxplot!([4], Timing_Max_Flows_future85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    violin!([1], Timing_Max_Flows_past45, color=[Farben45[1]], alpha=0.5)
    violin!([2], Timing_Max_Flows_past85, color=[Farben45[1]], alpha=0.5)
    violin!([3], Timing_Max_Flows_future45,color=[Farben45[2]], alpha=0.5)
    #boxplot!(Timing_Max_Flows_past85, color=[Farben85[1]])
    violin!([4], Timing_Max_Flows_future85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px, alpha=0.5)
    xticks!([1:4;], ["Past45", "Past8.5", "Future 4.5", "Future 8.5"])
    ylabel!("Timing of Maximum Discharge")
    #ylims!((2,10))
    yticks!([1,32,60,91,121,152,182,213,244,274,305,335], ["1.1", "1.2", "1.3", "1.4", "1.5", "1.6", "1.7","1.8", "1.9", "1.10", "1.11", "1.12"])
    title!("Timing of Maximum Discharge")
    title!("past45 "*string(round(median(Timing_Max_Flows_past45)))*"past85 "*string(round(median(Timing_Max_Flows_past85)))*"future45 "*string(round(median(Timing_Max_Flows_future45)))*"future85 "*string(round(median(Timing_Max_Flows_future85))))
    println("percentiles 25 past ", mean([quantile(Timing_Max_Flows_past45, 0.25), quantile(Timing_Max_Flows_past85, 0.25)]))
    println("percentiles 75 past ", mean([quantile(Timing_Max_Flows_past45, 0.75), quantile(Timing_Max_Flows_past85, 0.75)]))
    println("percentiles 25 future ", quantile(Timing_Max_Flows_future45, 0.25), " ", quantile(Timing_Max_Flows_future85, 0.25))
    println("percentiles 75 future ", quantile(Timing_Max_Flows_future45, 0.75), " ", quantile(Timing_Max_Flows_future85, 0.75))
    println("mean past45 "*string(round(mean(Timing_Max_Flows_past45)))*"past85 "*string(round(mean(Timing_Max_Flows_past85)))*"future45 "*string(round(mean(Timing_Max_Flows_future45)))*"future85 "*string(round(mean(Timing_Max_Flows_future85))))
    println("mean past45 "*string(round(std(Timing_Max_Flows_past45)))*"past85 "*string(round(std(Timing_Max_Flows_past85)))*"future45 "*string(round(std(Timing_Max_Flows_future45)))*"future85 "*string(round(std(Timing_Max_Flows_future85))))
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_timing_max_yearly_discharge_new.png")

    println("past 45 ", average_timing_simulations(Timing_Max_Flows_past45))
    println("past 85 ", average_timing_simulations(Timing_Max_Flows_past85))
    println("future 45 ", average_timing_simulations(Timing_Max_Flows_future45))
    println("future 85 ", average_timing_simulations(Timing_Max_Flows_future85))
    #absolute and relative change in timing of low flows
    #   dates have to be transformed to circular coordinates
    Difference_Timing_45 = difference_timing(Timing_Max_Flows_past45, Timing_Max_Flows_future45)
    Difference_Timing_85 = difference_timing(Timing_Max_Flows_past85, Timing_Max_Flows_future85)
    boxplot(Difference_Timing_45,color=[Farben45[2]])
    boxplot!(Difference_Timing_85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("absolute change [days]")
    title!("Change in Occurence of Average Maximum Annual Discharge [days]")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_change_timing_max_yearly_discharge.png")
    #
    violin(Difference_Timing_45,color=[Farben45[2]])
    violin!(Difference_Timing_85, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    xticks!([1:2;], ["RCP 4.5", "RCP 8.5"])
    ylabel!("absolute change [days]")
    title!("Change in Occurence of Average Maximum Annual Discharge [days]")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_change_timing_max_yearly_discharge_violin.png")
end

function plot_Max_Flows_Prob_Distribution(Max_Flows_past45, Max_Flows_future45, Max_Flows_past85, Max_Flows_future85, Exceedance_Probability, Catchment_Name)
    plot()
    Farben45=palette(:blues)
    Farben85=palette(:reds)
    for exceedance in collect(5:5:30)
        index = findall(x->x == Exceedance_Probability[exceedance], Exceedance_Probability)
        boxplot!(relative_error(Max_Flows_future45[index], Max_Flows_past45[index])*100,color=[Farben45[2]])
        #boxplot!(, color=[Farben85[1]])
        boxplot!(relative_error(Max_Flows_future85[index], Max_Flows_past85[index])*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    end
    #xticks!([1.5:2:11.5;], ["5 years", "10 years", "15 years", "20 years", "25 years", "30 years"])
    xticks!([1.5:2:11.5;], ["16 %", "32 %", "48 %", "65 %", "81 %", "97 %"])
    ylabel!("relative change in discharge [%]")
    title!("Relative Change in Maximum Annual Discharge which is exceeded in x % of the years")
    xlabel!("Years in 30 years")
    hline!([0], color=["grey"], linestyle = :dash)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_probability_Distribution_relative_Change.png")
    plot()
    for exceedance in collect(5:5:30)
        index = findall(x->x == Exceedance_Probability[exceedance], Exceedance_Probability)
        #change = relative_error(Max_Flows_future45[exceedance], Max_Flows_past45[exceedance])*100
        #print(mean(change), " ", maximum(change), " ", minimum(change), "\n")
        violin!(relative_error(Max_Flows_future45[index], Max_Flows_past45[index])*100,color=[Farben45[2]])
        #boxplot!(, color=[Farben85[1]])
        violin!(relative_error(Max_Flows_future85[index], Max_Flows_past85[index])*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    end
    #xticks!([1.5:2:11.5;], ["5 years", "10 years", "15 years", "20 years", "25 years", "30 years"])
    xticks!([1.5:2:11.5;], ["16 %", "32 %", "48 %", "65 %", "81 %", "97 %"])
    ylabel!("relative change in discharge [%]")
    xlabel!("Years in 30 years")
    hline!([0], color=["grey"], linestyle = :dash)
    #ylims!(-150,350)
    #title!("Relative Change in Maximum Annual Discharge which is exceeded in x years of the 30 year period")
    title!("Relative Change in Maximum Annual Discharge which is exceeded in x % of the years")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_probability_Distribution_relative_Change_violin.png")

    plot()
    mean_change = Float64[]
    max_change = Float64[]
    min_change = Float64[]
    mean_change_85 = Float64[]
    max_change_85 = Float64[]
    min_change_85 = Float64[]
    std_change_45 = Float64[]
    std_change_85 = Float64[]
    for exceedance in Exceedance_Probability[1:30]
        index = findall(x->x == exceedance, Exceedance_Probability)
        change_45 = relative_error(Max_Flows_future45[index], Max_Flows_past45[index])*100
        change_85 = relative_error(Max_Flows_future85[index], Max_Flows_past85[index])*100
        append!(mean_change, mean(change_45))
        append!(max_change, maximum(change_45))
        append!(min_change, minimum(change_45))
        append!(mean_change_85, mean(change_85))
        append!(max_change_85, maximum(change_85))
        append!(min_change_85, minimum(change_85))
        append!(std_change_45, std(change_45))
        append!(std_change_85, std(change_85))
        #boxplot!(, color=[Farben85[1]])
        #violin!(relative_error(Max_Flows_future85[index], Max_Flows_past85[index])*100, size=(1000,500), leg=false, left_margin = [5mm 0mm], xrotation = 60, color=[Farben85[2]], bottom_margin = 20px)
    end
    plot(collect(1/31:1/31:30/31)*100, mean_change, color=[Farben45[2]], label="RCP 4.5", ribbon = (mean_change - min_change, max_change - mean_change))
    plot!(collect(1/31:1/31:30/31)*100, mean_change_85, color=[Farben85[2]], label="RCP 8.5", ribbon = (mean_change_85 - min_change_85, max_change_85 - mean_change_85), size=(1600,800))
    #xticks!([1.5:2:11.5;], ["5 years", "10 years", "15 years", "20 years", "25 years", "30 years"])
    ylabel!("relative change in discharge [%]")
    #ylims!(-150,400)
    xlabel!("Years in 30 years")
    title!("Relative Change in Maximum Annual Discharge which is exceeded in x % of years in the 30 year period")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_probability_Distribution_relative_Change_magnitude.png")

    plot()
    plot(collect(1/31:1/31:30/31)*100, mean_change, color=[Farben45[2]], label="RCP 4.5", ribbon = std_change_45)
    plot!(collect(1/31:1/31:30/31)*100, mean_change_85, color=[Farben85[2]], label="RCP 8.5", ribbon = std_change_85, size=(1600,800))
    #xticks!([1.5:2:11.5;], ["5 years", "10 years", "15 years", "20 years", "25 years", "30 years"])
    ylabel!("relative change in discharge [%]")
    #ylims!(-150,400)
    xlabel!("Years in 30 years")
    title!("Relative Change in Maximum Annual Discharge which is exceeded in x % of years in the 30 year period, Filled: 1 Std")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_probability_Distribution_relative_Change_magnitude_std.png")

    return mean_change, min_change, max_change, std_change_45

end

"""
For each timeseries it calculates the number of times the maximum annual discharge occurs within a timerange (e.g. 15days)

$(SIGNATURES)

The function returns probability of occurence of AMF within a certain timerange of the year, and an array of this timerange
"""
function get_distributed_dates(Date_Past, Timerange, nr_runs, nr_years)
    nr_yearly_max_period_15_days = Float64[]
    day_range = Float64[]
    for i in 1:14*nr_runs
        Current_Date_Past = Date_Past[1+(i-1)*nr_years:nr_years*i]
        for days in 1:Timerange:366
            current_days = filter(Current_Date_Past) do x
                x >= days && x < days +Timerange
            end
            append!(day_range, days-1)
            if current_days == Float64[]
                append!(nr_yearly_max_period_15_days, 0)
            else
                append!(nr_yearly_max_period_15_days, length(current_days)/nr_years)
            end
        end
    end
    return nr_yearly_max_period_15_days, day_range
end
"""
For each timeseries it calculates the number of times the maximum annual discharge occurs within a timerange (e.g. 15days)

$(SIGNATURES)

The function returns probability of occurence of AMF within a certain timerange of the year, and an array of this timerange
"""
function plot_change_timing_AMF_over_year(Date_Past, Date_Future, Catchment_Name, nr_runs, rcp)
    period_15_days_past, day_range_past = get_distributed_dates(Date_Past, 15, nr_runs, 30)
    period_15_days_future, day_range_future = get_distributed_dates(Date_Future, 15, nr_runs, 30)
    #change = period_15_days_future - period_15_days_past
    plot()
    for i in collect(0:15:366)
        current_past = period_15_days_past[findall(x->x==i, day_range_future)]
        current_future = period_15_days_future[findall(x->x==i, day_range_future)]
        violin!(current_past*100, leg=false, size=(1500,800), color="blue")
        violin!(current_future*100, leg=false, size=(1500,800), color="red", left_margin = [5mm 0mm], bottom_margin = 20px, xrotation = 60)
    end
    ylabel!("Probability of Occurence in Timeseries [%]")
    xlabel!("15 days timesteps in year")
    title!("Timing of Maximum Annual Discharge,Blue=Past Red=Future")
    xticks!([1.5:2:48.5;],["Begin Jan", "End Jan", "Begin Feb", "End Feb", "Begin Mar", "End Mar", "Begin Apr", "End Apr", "Begin May", "End May", "Begin June", "End June","Begin Jul", "End Jul", "Begin Aug", "Eng Aug", "Begin Sep", "End Sep", "Begin Oct", "End Oct", "Begin Nov", "End Nov", "Begin Dec", "End Dec"])
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_probability_Distribution_timing_"*rcp*"_violin.png")
    plot()
    for i in collect(0:15:366)
        current_past = period_15_days_past[findall(x->x==i, day_range_future)]
        current_future = period_15_days_future[findall(x->x==i, day_range_future)]
        boxplot!(current_past*100, leg=false, size=(1500,800), color="blue")
        boxplot!(current_future*100, leg=false, size=(1500,800), color="red", left_margin = [5mm 0mm], bottom_margin = 20px, xrotation = 60)
    end
    ylabel!("Probability of Occurence in Timeseries [%]")
    xlabel!("15 days timesteps in year")
    title!("Timing of Maximum Annual Discharge,Blue=Past Red=Future")
    xticks!([1.5:2:48.5;],["Begin Jan", "End Jan", "Begin Feb", "End Feb", "Begin Mar", "End Mar", "Begin Apr", "End Apr", "Begin May", "End May", "Begin June", "End June","Begin Jul", "End Jul", "Begin Aug", "Eng Aug", "Begin Sep", "End Sep", "Begin Oct", "End Oct", "Begin Nov", "End Nov", "Begin Dec", "End Dec"])
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Annual_Max_Discharge/"*Catchment_Name*"_probability_Distribution_timing_"*rcp*".png")

end

function std_timing(All_Concentration_past_45, All_Concentration_future_45, All_Concentration_past_85, All_Concentration_future_85, Category)
    std_occurence_past_45 = sqrt.(-2*log.(All_Concentration_past_45)).*(365.23333/(2*pi))
    std_occurence_future_45 = sqrt.(-2*log.(All_Concentration_future_45)).*(365.23333/(2*pi))
    std_occurence_past_85 = sqrt.(-2*log.(All_Concentration_past_85)).*(365.23333/(2*pi))
    std_occurence_future_85 = sqrt.(-2*log.(All_Concentration_future_85)).*(365.23333/(2*pi))
    plot()
    scatter(std_occurence_future_45)
    xlabel!("Nr Runs")
    ylabel!("Std in Days")
    title!("Standard Deviation of Timing between years in timeperiod")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/"*Category*"/std_timing_future_45.png")
    scatter(std_occurence_future_85)
    xlabel!("Nr Runs")
    ylabel!("Std in Days")
    title!("Standard Deviation of Timing between years in timeperiod")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/"*Category*"/std_timing_future_85.png")
    plot()
    scatter(std_occurence_past_45)
    xlabel!("Nr Runs")
    ylabel!("Std in Days")
    title!("Standard Deviation of Timing between years in timeperiod")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/"*Category*"/std_timing_past_45.png")
    scatter(std_occurence_past_85)
    xlabel!("Nr Runs")
    ylabel!("Std in Days")
    title!("Standard Deviation of Timing between years in timeperiod")
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/"*Category*"/std_timing_past_85.png")
end

# precipitation statistics monthy

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 Distributions, Random, Test

@testset "truncated exponential" begin
    d = Exponential(1.5)
    l = 1.2
    r = 2.7
    @test mean(d) ≈ mean(truncated(d; lower=-3.0)) # non-binding truncation
    @test mean(d) ≈ mean(truncated(d; lower=-3.0, upper=Inf))
    @test mean(truncated(d; lower=l)) ≈ mean(d) + l
    @test mean(truncated(d; lower=l, upper=Inf)) ≈ mean(d) + l
    # test values below calculated using symbolic integration in Maxima
    @test mean(truncated(d, 0, r)) ≈ 0.9653092084094841
    @test mean(truncated(d, l, r)) ≈ 1.82703493969601

    # all the fun corner cases and numerical quirks
    @test mean(truncated(Exponential(1.0); upper=0)) == 0                   # degenerate
    @test mean(truncated(Exponential(1.0); lower=-Inf, upper=0)) == 0
    @test mean(truncated(Exponential(1.0); upper=0+eps())) ≈ 0 atol = eps() # near-degenerate
    @test mean(truncated(Exponential(1.0); lower=-Inf, upper=0+eps())) ≈ 0 atol=eps()
    @test mean(truncated(Exponential(1.0), 1.0, 1.0+eps())) ≈ 1.0 # near-degenerate
    @test mean(truncated(Exponential(1e308), 1.0, 1.0+eps())) ≈ 1.0 # near-degenerate
end

ext: jl


text:
"Soft Boolean.  Value in [0, 1]"
struct SoftBool{ET <: Real} <: AbstractSoftBool
  logerr::ET
  SoftBool(l::T) where {T} = new{T}(l)
  SoftBool(l::T) where {T <: ForwardDiff.Dual} = new{T}(l) # Resolves type ambiguity (from flux)
end
@invariant 0 <= err(b::SoftBool) <= 1

"Error in [0, 1]"
err(x::SoftBool) = exp(x.logerr)

"Log error in `[-Inf, 0]`"
logerr(x::SoftBool) = x.logerr
Bool(x::SoftBool) = logerr(x) == 0.0
ssofttrue(::Type{T} = Float64) where T = SoftBool(zero(t))
ssoftfalse(::Type{T} = Float64) where T = SoftBool(-inf(T))

# (In)Equalities #

"Kernel return type as function of arguments"
kernelrettype(x::T, y::T) where T = T

softeq(a::Bool, b::Bool) = softeq(float(a), float(b))
softeq(a::Integer, b::Integer) = softeq(float(a), float(b))

"Soft Equality"
function ssofteq(x, y, k = globalkernel())
  r = d(x, y)
  SoftBool(k(r)::typeof(r))
end

"Soft >"
function ssoftgt(x::Real, y::Real, k = globalkernel())
  r = bound_loss(x, y, Inf)
  SoftBool(k(r)::typeof(r))
end

"Soft <"
function ssoftlt(x::Real, y::Real, k = globalkernel())
  r = bound_loss(x, -Inf, y)
  SoftBool(k(r)::typeof(r))
end

# Boolean Operators #
function Base.:&(x::SoftBool, y::SoftBool)
  a = logerr(x)
  b = logerr(y)
  # c = min(a, b)
  c = a + b
  SoftBool(c)
end
Base.:|(x::SoftBool, y::SoftBool) = SoftBool(max(logerr(x), logerr(y)))
# zt: FIXME abstractvector?
Base.all(xs::Vector{<:SoftBool}) = SoftBool(minimum(logerr.(xs)))
Base.all(xs::Vector{<:RandVar}) = RandVar(all, (xs, ))

# Arithmetic #
Base.:*(x::SoftBool{T}, y::T) where T <: Real = SoftBool(x.logerr * y)
Base.:*(x::T, y::SoftBool{T}) where T <: Real = SoftBool(x * y.logerr)

# Show #
Base.show(io::IO, sb::SoftBool) = print(io, "ϵ:$(logerr(sb))")


text:
using AlgebraPDF
using Test

#              _|_|_|
#    _|_|_|  _|          _|_|_|  _|    _|    _|_|_|
#  _|    _|  _|  _|_|  _|    _|  _|    _|  _|_|
#  _|    _|  _|    _|  _|    _|  _|    _|      _|_|
#    _|_|_|    _|_|_|    _|_|_|    _|_|_|  _|_|_|

@testset "Standard functions" begin
    @test AlgebraPDF.standardgauss(0,1) ≈ 1/sqrt(2π)
    @test AlgebraPDF.standardgauss(2,2) ≈ exp(-1/2)/sqrt(2π*4)
    x, m, Γ = 1.0, 1.0, 0.1
    @test AlgebraPDF.amplitudeBW(x, m, Γ) ≈ 1im
    @test AlgebraPDF.amplitudeBWsq(x, m, Γ) ≈ 1.0
    #
    x, σ, r, n = 0, 0.1, 0.3, 5
    @test AlgebraPDF.standarddoublegauss(x,σ,r,n) ≈ 1/sqrt(2π)/σ*(r+(1-r)/n)
    #
end

@testset "Precodded pdfs" begin
    dGauss = Normalized(FGauss((μ1=1.0,σ1=0.3)), (-2,2))
    @test keys(pars(dGauss)) == (:μ1,:σ1)
    exp_gauss_at_zero = 1/sqrt(2π)/pars(dGauss).σ1
    @test abs(dGauss(1.0) - exp_gauss_at_zero) / exp_gauss_at_zero < 1e-3
end

@testset "Precodded pdfs" begin
    d2 = FBreitWigner((mΩb = 6030, Γ=17.0))
    @test freepars(d2) == (mΩb = 6030, Γ=17.0)
    @test real(d2(freepars(d2).mΩb)) ≈ 0
    #
    d3 = FExp((τ = -1.1,))
    @test freepars(d3) === (τ = -1.1,)
    #
    d4 = FPowExp((n=3.3, τ=-4.0))
    @test d4(1.1) != 0.0
    @test freepars(d4) == (n=3.3, τ=-4.0)
    #
    pars3 = (c0=1.1, c1=2.2, c2=3.3, c4=4.4)
    d5 = FPol(pars3)
    @test func(d5, 1.1; p=pars3) == sum(1.1^(i-1)*c for (i,c) in enumerate(pars3))
    @test freepars(d5) == pars3
    #
    d6 = FDoubleGaussFixedRatio((m = 0.77, Γ=0.15, r=0.8, n=3))
    @test d6(0.77) != 0.0
    #
    d7 = FBreitWignerConvGauss((m = 0.77, Γ=0.15, σ=0.03))
    @test d7(0.77) != 0.0
    #
    xv = range(-π, 2π, length=40)
    yv = map(x->x*cos(3x) - 3*sin(x), xv)
    d8 = FTabulated(xv,yv)
    @test length(freepars(d8)) == 0
    @test d8(1.1) != 0.0
    @test d8( 3π) == 0.0
    @test d8( -π) != 0.0
    @test d8( 2π) != 0.0
end

# here is MWE of the normalization problem
# let
#     d = aGauss((a=0.01,b = 1.1), (-3,3))
#     d(1.1, [36.1,0.1])
# end


text:
module MySignalProcessing|

# Write your package code here.

# Dependancies
include("types.jl")
using .types: signal

# imports

# Importiong packages that will then be extended.
import Base.+
import Base.*
import DSP

#exports
export impseq
export signal
export stepseq
export sigshift
export realexp
export sigadd
export sigmult
export sigrand
export sigscale
export sigfold
export energy
export +
export *
export sinseq
export sinseqmn
export conv
export cexp

"""
    δ(start::Int, stop::Int)

Unit impulse signal
This function has the form of :

    δ(n) = 1 if n=0
         = 0 if n!=0

# parameters
    `start:: Int` : start of the impulse sequence n (horizontal axis)
    `stop:: Int` : stop of the impulse sequence n (horizontal axis)

# returns
    s::signal
"""
function impseq(start::Int, stop::Int)
    s = signal(0.0, start,stop); #constructor for initial signal

    for x in 1:length(s.n)

        s.A[x] = Float64(s.n[x] == 0) # elegant way to return 1

    end

    return s

end

"""
    stepseq(start::Int, stop::Int)

Unit Step Sequence
This function has the form of :

    stepseq(n) = 1 if n ≥ 0
                = 0 if n < 0

# parameters
    `start:: Int` : start of the sampple position n (horizontal axis)
    `stop:: Int` : stop of the sample position n (horizontal axis)

# returns
    s::signal

"""
function stepseq(start::Int, stop::Int)
    s = signal(0.0, start,stop); #constructor for initial signal

    for x in 1:length(s.n)

        s.A[x] = Float64(s.n[x] >= 0) # elegant way to return 1

    end

    return s

end

"""
    _sigshift(s::signal, m::Int)

Function to shift signal s by adding m to each element of s.n
It implements y(n) = x(n-n₀)

# parameters
    `s::signal` : signal to be shifted
    `n₀::Int` : shift amount

# returns
    s'::signal shifted
"""
function _sigshift(s::signal, n₀::Int)
    y = s.A;
    n = s.n .+ n₀;
    s = signal(y, n);
    return s

end

"""
    sigshift(n₀::Int)

# parameters
    `n₀::Int` : shift amount

# returns
    s:: anonimous  function that shifts a signal (designed to be composable)

"""
function sigshift(n₀::Int)

    return (sig -> _sigshift(sig, n₀))

end

"""
    realexp(a::Real, n₀::Int , nf::Int)

Function to generate a real exponential signal.

Signal of the form:
```math
r(a, n) = a^n
```

# parameters
    `a::Real` : real-valued parameter of the exponential
    `n₀::Int` : start  the sample positon signal
    `nf::Int` : stop  of the sample positon  signal

# returns
    s::signal

 """
function realexp(a::Real, n₀::Int, nf::Int)
    s = signal(0.0, n₀, nf); #constructor for initial signal

    for x in 1:length(s.n)

        s.A[x] = a^(s.n[x])

    end

    return s

end

"""
    sigadd(s1::signal, s2::signal)

Function to add two signal s1 and s2.
pls note that the signals must have the same length
If sequences are of unequal lenghts or if the sample positions are
    different for equal-length sequences, the we cannot directly add
        We have to first augment x1(n) and x2(n) , so that they have the same
         position vector n (and hence the same length)

# parameters
    `s1::signal` : signal to be added
    `s2::signal` : signal to be added

# returns
    `s::signal` : signal s1 + s2

"""
function sigadd(s₁::signal  , s₂::signal)

    n₀ =min(minimum(s₁.n), minimum(s₂.n)) # find the minimum sample position
    nf = max(maximum(s₁.n), maximum(s₂.n)) # find the maximum sample position
    n=n₀:nf; # create the sample position vector

    y₁ = signal(0.0, n₀, nf); #constructor for initial signal
    y₂ = signal(0.0, n₀, nf); #constructor for initial signal

    filter1 = findall((n .≥ minimum(s₁.n)) .& (n .≤ maximum(s₁.n)) .== 1)
    filter2 = findall((n .≥ minimum(s₂.n)) .& (n .≤ maximum(s₂.n)) .== 1)

    y₁.A[filter1] = s₁.A; # fill the signal with the values of s₁
    y₂.A[filter2] = s₂.A; # fill the signal with the values of s₂

    s = signal(0.0, n₀, nf); #constructor for initial signal
    s.A = y₁.A + y₂.A; # add the two signals

    return s

end

"""
    Base.+(s₁::signal  , s₂::signal)

Extends Base.+ to add two signal objects together.

# Arguments

    `s₁::signal` : signal to be added
    `s₂::signal` : signal to be added

# Returns
    `s::signal` : signal s1 + s2
"""
function +(s₁::signal  , s₂::signal)
    s = sigadd(s₁ , s₂)
    return s
end

"""
    sigmult(s1::signal, s2::signal)

Function to multiply two signal s1 and s2.
pls note that the signals must have the same length
If sequences are of unequal lenghts or if the sample positions are
    different for equal-length sequences, the we cannot directly add
        We have to first augment x1(n) and x2(n) , so that they have the same
         position vector n (and hence the same length)

# parameters
    `s1::signal` : signal to be added
    `s2::signal` : signal to be added

# returns
    `s::signal` : signal s1 * s2

"""
function sigmult(s₁::signal  , s₂::signal)

    n₀ =min(minimum(s₁.n), minimum(s₂.n)) # find the minimum sample position
    nf = max(maximum(s₁.n), maximum(s₂.n)) # find the maximum sample position
    n=n₀:nf; # create the sample position vector

    y₁ = signal(0.0, n₀, nf); #constructor for initial signal
    y₂ = signal(0.0, n₀, nf); #constructor for initial signal

    filter1 = findall((n .≥ minimum(s₁.n)) .& (n .≤ maximum(s₁.n)) .== 1)
    filter2 = findall((n .≥ minimum(s₂.n)) .& (n .≤ maximum(s₂.n)) .== 1)

    y₁.A[filter1] = s₁.A; # fill the signal with the values of s₁
    y₂.A[filter2] = s₂.A; # fill the signal with the values of s₂

    s = signal(0.0, n₀, nf); #constructor for initial signal
    s.A = y₁.A .* y₂.A; # add the two signals

    return s

end

"""
    *(s₁::signal  , s₂::signal)

Multiples two signal objects together.

# Arguments

    `s₁::signal` : signal to be multiplied
    `s₂::signal` : signal to be Multiples

# Returns
    `s::signal` : signal s1 * s2

"""
function *(s₁::signal  , s₂::signal)
    s = sigmult(s₁ , s₂)
    return s
end

"""
    sigrand(start::Int, stop::Int)

Function to generate a random signal starting at sample position start and ending at sample position stop.
The function generates a length N random sequence whose elements are uniformely distributed between 0 and 1.

# parameters
    `start::Int` : start  the sample positon signal
    `stop::Int` : stop  of the sample positon  signal

# returns
    s::signal

"""
function sigrand(start::Int, stop::Int)
    n = start:stop |> collect; # create the sample position vector
    y = rand.(length(n)); # generate a random sequence
    s = signal(y, n); #constructor for initial signal
    return s

end

"""
    _sigscale(s::signal, α::Real)

Each signal sample is multiple by a scalar α

# parameters
    `s::signal` : signal to be scaled
    `α::Real` : scalar

# returns
    `s::signal` : scaled signal
"""
function _sigscale(s::signal, α::Real)

    s1 = signal(s.A .* α, s.n); #constructor for initial signal

    return s1

end

"""
    sigscale(α::Real)

Each signal sample is multiple by a scalar α

# parameters
    `α::Real` : scalar

# returns
    anonimous function that takes a signal as argument and returns
         scaled signal (designed to be composable)

"""
function sigscale(α::Real)

    return (sig -> _sigscale(sig, α))
end

"""
    *(s₁::signal  , α::Real)

    scale signal by scalar.

# Arguments

    `s₁::signal` : signal to be multiplied
    `α::Real` : scalar

# Returns
    `s::signal` : signal s1 * α

"""
function *(s₁::signal  , α::Real)
    s = s₁ |> sigscale(α);
    return s
end

"""
    *(α::Real, s₁::signal)

scale signal by scalar.

# Arguments

    `s₁::signal` : signal to be multiplied
    `α::Real` : scalar

# Returns
    `s::signal` : signal s1 * α

"""
function *(α::Real, s₁::signal)
    s = s₁ |> sigscale(α);
    return s
end

"""
    sigfold(s::signal)

In this operation, each sample of x(n) is flipped around n=0 to obtain a folded sequence  y(n)

y(n) = {x(-n)}

# parameters
    `s::signal` : signal to be folded

# returns
    `s::signal` : folded signal

"""
function  sigfold(s::signal)
    s1 = signal(s.A, s.n); #constructor for initial signal
    s1.A = reverse(s.A);
    s1.n = -reverse(s.n);
    return s1
end

"""
    energy(s::signal)

Calculate the energy of a signal
The energy of a signal is given by the sum of the
    absolute squared values of the signal samples

# parameters
    `s::signal` : input signal

# returns
    `e::Real` : energy of the signal

"""
function energy(s::signal)
    e = sum(abs.(s.A) .^ 2);
    return e
end

"""
    sinSeq(start::Int, stop::Int, ω::Real, φ::Real)

Sinusoidal function

# parameters
    `start::Int` : start  the sample positon signal
    `stop::Int` : stop  of the sample positon  signal
    `ω::Real` : frequency of the sinusoid in radians
    `φ::Real` : phase of the sinusoid, in radians

# returns
    s::sinusoidal signal

"""
function sinseq(start::Int, stop::Int, ω::Real, φ::Real= 0.0)
    n = start:stop |> collect; # create the sample position vector
    y = sin.(n .* ω .+ φ); # generate a sinusoidal sequence
    s = signal(y, n); #constructor for initial signal
    return s
end

"""
    sinSeq(start::Int, stop::Int, M::Int,N::Int, φ::Real)

Sinusoidal function with more specific constructor

# parameters
    `start::Int` : start  the sample positon signal
    `stop::Int` : stop  of the sample positon  signal
    `M::Int` : multiple of 2π
    `N::Int` : number of periods before pattern repeats itself

# returns
    s::sinusoidal signal

"""
function sinseqmn(start::Int, stop::Int, M::Int,N::Int, φ::Real= 0.0)
    n = start:stop |> collect; # create the sample position vector
    y = sin.(n .* ((M / N)*2π .+ φ)); # generate a sinusoidal sequence
    s = signal(y, n); #constructor for initial signal
    return s
end

"""
    conv(s₁::signal, s₂::signal)

Convolution function

# parameters
    `s₁::signal` : input signal
    `s₂::signal` : impulse response signal

# returns
    s::convolution signal output

"""
function conv(s₁::signal, s₂::signal)
    n₁ = s₁.n;
    n₂ = s₂.n;
    A₁ = convert(Array{Float64,1}, s₁.A) #conv requires more specific type
    A₂ = convert(Array{Float64,1}, s₂.A) #conv requires more specific type

    Nyb = n₁[1] + n₂[1];
    Nye = n₁[length(A₁)] + n₂[length(A₂)];

    ny = collect(Nyb: Nye);
    y = DSP.conv(A₁, A₂);

    return signal(y, ny);
  end

  """
  function cexp(A::Vector, k, N)

      params:
          - A: Vector of amplitudes
          - n: element index within each Vector, from 0 .. N-1
          - k: index indicating each vector , from 0 .. N-1
          - N: dimension of signal (number of vectors)

      returns:
          - complex exponential signal

  """
  function cexp(A::Vector, k, N = 64)

      n = [n for n in range(0, N-1)]

      w = [exp(((2pi/64)*im)*n*k) for n in range(0, N-1)]
      s = A .*w

      return signal(s, n)

  end

end

max_stars_repo_name: gjunqueira-sys/MySignalProcessing

max_issues_repo_name: gjunqueira-sys/MySignalProcessing

max_forks_repo_name: gjunqueira-sys/MySignalProcessing


text:
# OPTIMO_FEAS
#
# PART OF OPTIMO

"""
    FeasOptiModel( prob, [x0, with_indicator, xp, dp, wp, name] )
    represents the following structured, proximal, feasibility problem:
    ```
    minimize  Φ(x) + with_indicator ind_g( x )
    ```
    with respect to `x`, starting from `x0`, where
    ```
    Φ(x)  :=  m( c(x) - proj_S( c(x) ) ) + wprox/2 * || dprox .* (x - xprox) ||^2
    ```
    and `m()` is a loss function, e.g., L1-norm, L2-norm, and Huber loss.

    Optional arguments:
    * `x0::Vector` (default: `xprox`) initial guess
    * `with_indicator::Bool` (default: `false`) include proximable term as indicator?
    * `xprox::Vector` (default: `prob.meta.x0`) proximal point
    * `dprox::Vector` (default: `ones(R,nvar)`) vector for nonnegative, diagonal scaling
    * `wprox::Real` (default: `0`) nonnegative scalar value for proximal regularization
    * `name::String` (default: `prob.meta.name * "-feas"`) model name
"""

using LinearAlgebra: norm

export FeasOptiModel
export infeasibility, infeasibilitygrad!
export cviolation, proxdistance, unsproxdistance

mutable struct FeasOptiModel <: AbstractOptiModel
    meta::OptiModelMeta
    base::AbstractOptiModel
    xprox::AbstractVector
    dprox::AbstractVector
    wprox::Real
    with_prox::Bool
    with_indicator::Bool
    loss::Symbol
    huber_rho::Real
    huber_mu::Real
    # additional allocation
    _cx::AbstractVector
    _px::AbstractVector

    function FeasOptiModel( meta::OptiModelMeta,
                            base::TP,
                            xprox::Tx,
                            dprox::Tx,
                            wprox::Real,
                            with_prox::Bool,
                            with_indicator::Bool,
                            loss::Symbol,
                            huber_rho::Real,
                            huber_mu::Real,
    ) where {TP <: AbstractOptiModel, Tx <: AbstractVector}
        R = eltype(base.meta.x0)
        ncon = base.meta.ncon
        new( meta, base, xprox, dprox, wprox, with_prox, with_indicator, loss,
             huber_rho, huber_mu, zeros(R,ncon), zeros(R,ncon) )
    end
end

function FeasOptiModel( prob::TP;
                        xprox::Tx=prob.meta.x0,
                        dprox::Tx=ones(eltype(prob.meta.x0),prob.meta.nvar),
                        wprox::Real=zero(eltype(prob.meta.x0)),
                        x0::Tx=xprox,
                        with_indicator::Bool=false,
                        loss::Symbol=:l2sq,
                        huber_rho::Real=1.0,
                        huber_mu::Real=1.0,
                        name::String=prob.meta.name * "-feas",
        ) where {TP <: AbstractOptiModel, Tx <: AbstractVector}
    @assert length(x0) == prob.meta.nvar
    @assert length(xprox) == prob.meta.nvar
    @assert length(dprox) == prob.meta.nvar
    @assert 0 <= wprox
    @assert all( 0 .<= dprox )
    @assert loss ∈ [:l1, :l2, :l2sq, :huber]
    @assert 0 < huber_rho
    @assert 0 < huber_mu
    meta = OptiModelMeta( prob.meta.nvar, 0, x0=x0, minimize=true, name=name )
    with_prox = wprox > 0 && all( dprox .> 0 )
    fprob = FeasOptiModel( meta, prob, xprox, dprox, wprox, with_prox,
                           with_indicator, loss, huber_rho, huber_mu )
    finalizer( p -> finalize(p.base), fprob )
    return fprob
end

##########################################################
# Methods
##########################################################
# obj, grad!, prox!, objprox!
"""
    obj( fprob, x )
evaluates Φ at x.
"""
function obj(prob::FeasOptiModel, x::AbstractVector)
    Φ = infeasibility( prob, x )
    if prob.with_prox
        Φ += prob.wprox * proxdistance( prob, x )
    end
    return Φ
end

"""
    grad!( fprob, x, dfx )
computes ∇Φ at x, in place.
"""
function grad!(prob::FeasOptiModel, x::AbstractVector, dfx::AbstractVector)
    @lencheck prob.meta.nvar dfx
    infeasibilitygrad!( prob, x, dfx )
    if prob.with_prox
        dfx .+= prob.wprox .* prob.dprox .* (x .- prob.xprox)
    end
    return nothing
end

"""
    prox!( fprob, x, a, z )
computes z := prox_{a g}(x), in place.
"""
function prox!(prob::FeasOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
    @lencheck prob.meta.nvar x z
    if prob.with_indicator
        prox!(prob.base, x, a, z )
    else
        copyto!(z, x)
    end
    return nothing
end

"""
    objprox!( fprob, x, a, z )
computes z := prox_{a g}(x) and g(z), in place, assuming g is an indicator.
"""
function objprox!(prob::FeasOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
    @lencheck prob.meta.nvar x z
    if prob.with_indicator
        prox!(prob.base, x, a, z )
    else
        copyto!(z, x)
    end
    return 0.0
end

##########################################################
"""
    infeasibility( fprob, x )
computes the infeasibility measure Φ at x, without proximal regularization (wprox = 0).
"""
function infeasibility( prob::FeasOptiModel, x::AbstractVector )
    evalinfeasvec( prob, x )
    if prob.loss == :l1
        return norm( prob._px, 1 )
    elseif prob.loss == :l2
        return norm( prob._px, 2 )
    elseif prob.loss == :l2sq
        return 0.5 * norm( prob._px, 2 )^2
    elseif prob.loss == :huber
        return huberloss( prob._px, prob.huber_rho, prob.huber_mu )
    end
end

"""
    infeasibility_grad!( fprob, x, dfx )
computes the gradient ∇Φ of the infeasibility measure Φ at x, without proximal
regularization (wprox = 0).
"""
function infeasibilitygrad!(prob::FeasOptiModel, x::AbstractVector, dfx::AbstractVector)
    evalinfeasvec( prob, x )
    if prob.loss == :l1
        prob._px .= sign.( prob._px )
        jtprod!( prob.base, x, prob._px, dfx )
    elseif prob.loss == :l2
        normx = norm( prob._px, 2 )
        if normx == 0
            dfx .= 0
        else
            prob._px ./= normx
            jtprod!( prob.base, x, prob._px, dfx )
        end
    elseif prob.loss == :l2sq
        jtprod!( prob.base, x, prob._px, dfx )
    elseif prob.loss == :huber
        huberlossgrad!( prob._px, prob.huber_rho, prob.huber_mu, prob._px )
        jtprod!( prob.base, x, prob._px, dfx )
    end
    return nothing
end

"""
    cviolation( fprob, x )
computes the constraint violation.
"""
function cviolation( prob::FeasOptiModel, x::AbstractVector )
    evalinfeasvec( prob, x )
    return norm( prob._px, Inf )
end

"""
    proxdistance( fprob, x )
computes the (scaled) distance to the proximal point.
"""
proxdistance, unsproxdistance
function proxdistance( prob::FeasOptiModel, x::AbstractVector )
    @lencheck prob.meta.nvar x
    return 0.5 * sum( prob.dprox .* (x .- prob.xprox).^2 )
end

"""
    unsproxdistance( fprob, x )
computes the unscaled distance to the proximal point.
"""
function unsproxdistance( prob::FeasOptiModel, x::AbstractVector )
    @lencheck prob.meta.nvar x
    return 0.5 * norm( x .- prob.xprox, 2)^2
end

##########################################################
# INTERNAL FUNCTIONS
##########################################################
function evalinfeasvec( prob::FeasOptiModel, x::AbstractVector )
    @lencheck prob.meta.nvar x
    cons!( prob.base, x, prob._cx )
    proj!( prob.base, prob._cx, prob._px )
    prob._px .= prob._cx .- prob._px
    return nothing
end

"""
    Huber loss, with parameter ρ > 0
        f(x) = { 1/2 ||x||^2           if ||x|| ⩽ ρ
               { ρ ( ||x|| - ρ/2 )     otherwise
"""
function huberloss( x::AbstractVector, rho::Real, mu::Real )
    normx = norm( x, 2 )
    if normx <= rho
        return 0.5 * mu * normx^2
    else
        return rho * mu * ( normx - 0.5 * rho )
    end
end

function huberlossgrad!( x::AbstractVector, rho::Real, mu::Real, dfx::AbstractVector )
    normx = norm( x, 2 )
    if normx <= rho
        dfx .= mu .* x
    else
        dfx .= ( rho * mu / normx ) .* x
    end
    return nothing
end


text:
"""
A `DataFrame` is already materialized, so this is a no-op.
"""
materialize(df::DataFrames.DataFrame) = df

"""
A `GroupedDataFrame` is already materialized, so this is a no-op.
"""
materialize(gdf::DataFrames.GroupedDataFrame) = gdf


text:
using Unreal
using Base.Test: @test

@test begin
    l = unreal(1)
    r = Int[]
    l == r
end

@test begin
    l = unreal(1.0)
    r = Float64[]
    l == r
end


text:
const TIME_ZONE_CACHE = Dict{String,Tuple{TimeZone,Class}}()

"""
    TimeZone(str::AbstractString) -> TimeZone

Constructs a `TimeZone` subtype based upon the string. If the string is a recognized
standard time zone name then data is loaded from the compiled IANA time zone database.
Otherwise the string is parsed as a fixed time zone.

A list of recognized standard and legacy time zones names can is available by running
`timezone_names()`. Supported fixed time zone string formats can be found in docstring for
[`FixedTimeZone(::AbstractString)`](@ref).

## Examples
```jldoctest
julia> TimeZone("Europe/Warsaw")
Europe/Warsaw (UTC+1/UTC+2)

julia> TimeZone("UTC")
UTC
```
"""
TimeZone(::AbstractString)

"""
    TimeZone(str::AbstractString, mask::Class) -> TimeZone

Similar to [`TimeZone(::AbstractString)`](@ref) but allows you to control what time zone
classes are allowed to be constructed with `mask`. Can be used to construct time zones
which are classified as "legacy".

## Examples
```jldoctest
julia> TimeZone("US/Pacific")
ERROR: ArgumentError: The time zone "US/Pacific" is of class `TimeZones.Class(:LEGACY)` which is currently not allowed by the mask: `TimeZones.Class(:FIXED) | TimeZones.Class(:STANDARD)`

julia> TimeZone("US/Pacific", TimeZones.Class(:LEGACY))
US/Pacific (UTC-8/UTC-7)
```
"""
TimeZone(::AbstractString, ::Class)

function TimeZone(str::AbstractString, mask::Class=Class(:DEFAULT))
    # Note: If the class `mask` does not match the time zone we'll still load the
    # information into the cache to ensure the result is consistent.
    tz, class = get!(TIME_ZONE_CACHE, str) do
        tz_path = joinpath(TZData.COMPILED_DIR, split(str, "/")...)

        if isfile(tz_path)
            open(deserialize, tz_path, "r")
        elseif occursin(FIXED_TIME_ZONE_REGEX, str)
            FixedTimeZone(str), Class(:FIXED)
        elseif !isdir(TZData.COMPILED_DIR) || isempty(readdir(TZData.COMPILED_DIR))
            # Note: Julia 1.0 supresses the build logs which can hide issues in time zone
            # compliation which result in no tzdata time zones being available.
            throw(ArgumentError(
                "Unable to find time zone \"$str\". Try running `TimeZones.build()`."
            ))
        else
            throw(ArgumentError("Unknown time zone \"$str\""))
        end
    end

    if mask & class == Class(:NONE)
        throw(ArgumentError(
            "The time zone \"$str\" is of class `$(repr(class))` which is " *
            "currently not allowed by the mask: `$(repr(mask))`"
        ))
    end

    return tz
end

"""
    @tz_str -> TimeZone

Constructs a `TimeZone` subtype based upon the string at parse time. See docstring of
`TimeZone` for more details.

```julia
julia> tz"Africa/Nairobi"
Africa/Nairobi (UTC+3)
```
"""
macro tz_str(str)
    TimeZone(str)
end

"""
    istimezone(str::AbstractString, mask::Class=Class(:DEFAULT)) -> Bool

Check whether a string is a valid for constructing a `TimeZone` with the provided `mask`.
"""
function istimezone(str::AbstractString, mask::Class=Class(:DEFAULT))
    # Start by performing quick FIXED class test
    if mask & Class(:FIXED) != Class(:NONE) && occursin(FIXED_TIME_ZONE_REGEX, str)
        return true
    end

    # Perform more expensive checks against pre-compiled time zones
    tz, class = get(TIME_ZONE_CACHE, str) do
        tz_path = joinpath(TZData.COMPILED_DIR, split(str, "/")...)

        if isfile(tz_path)
            # Cache the data since we're already performing the deserialization
            TIME_ZONE_CACHE[str] = open(deserialize, tz_path, "r")
        else
            nothing, Class(:NONE)
        end
    end

    return tz !== nothing && mask & class != Class(:NONE)
end

# After `broadcastable` was defined but before this was addressed in the `Dates` stdlib
# - Broadcastable introduction: https://github.com/JuliaLang/julia/pull/26601
# - Fixed in Dates: https://github.com/JuliaLang/julia/pull/30159
# Note: The change was backported to 1.1 as well.
if v"0.7.0-DEV.4743" <= VERSION < v"1.1.0-DEV.722" || v"1.2-" <= VERSION < v"1.2.0-DEV.114"
    Base.broadcastable(tz::TimeZone) = Ref(tz)
end


text:
# This file was generated, do not modify it. # hide
using MLJ
X, y = @load_iris
@load DecisionTreeClassifier

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
@testset "calc" begin

# using prev defs from fix.jl

@with_names begin

    sum = Delayed()
    val = E"(" + sum + E")" | PFloat64()

    neg = Delayed()             # allow multiple negations (eg ---3)
    neg.matcher = Nullable{Matcher}(val | (E"-" + neg > Neg))

    mul = E"*" + neg
    div = E"/" + neg > Inv
    prd = neg + (mul | div)[0:end] |> Prd

    add = E"+" + prd
    sub = E"-" + prd > Neg
    sum.matcher = Nullable{Matcher}(prd + (add | sub)[0:end] |> Sum)

    all = sum + Eos()

end

println(all)

@test val.name == :val
@test typeof(val) == Alt
@test length(val.matchers) == 2
@test neg.name == :neg
@test typeof(neg) == Delayed
@test typeof(get(neg.matcher)) == Alt
@test length(get(neg.matcher).matchers) == 3  # flattening
@test prd.name == :prd
@test typeof(prd) == Transform
@test typeof(prd.matcher) == Seq
@test length(prd.matcher.matchers) == 2
@test typeof(prd.matcher.matchers[2]) == Depth
@test sum.name == :sum
@test typeof(sum) == Delayed
@test typeof(get(sum.matcher).matcher) == Seq
@test length(get(sum.matcher).matcher.matchers) == 2
@test typeof(get(sum.matcher).matcher.matchers[2]) == Depth

parse_dbg("1+2*3/4", Trace(all))

for (src, val) in [
                   ("1", 1),
                   ("-1", -1),
                   ("1+1", 2),
                   ("1-1", 0),
                   ("-1-1", -2)
                   ]
#    @test_approx_eq calc(parse_dbg(src, Trace(all))[1]) val
    @test calc(parse_one(src, Trace(all))[1]) ≈ val
    println("$src = $val")
end

for (src, ast, val) in
    [
     ("1.0", Sum([Prd([1.0])]), 1.0)
     ("-1.0", Sum([Prd([-1.0])]), -1.0)
     ("--1.0", Sum([Prd([Neg(-1.0)])]), 1.0)
     ("1+2", Sum([Prd([1.0]),Prd([2.0])]), 3.0)
     ("1+2*3/4", nothing, 2.5)
     ]
    if ast != nothing
        @test parse_one(src, all)[1] == ast
    end
    @test calc(parse_one(src, all)[1]) ≈ val
    println("$src = $val")
end

# some regression tests
@test calc(parse_one("-5.0/7.0+5.0-5.0", all)[1]) ≈ -0.7142857142857144
@test eval(Meta.parse("-5.0/7.0+5.0-5.0")) ≈ -0.7142857142857144
@test calc(parse_one("(0.0-9.0)", all)[1]) ≈ -9.0
@test calc(parse_one("((0.0-9.0))", all)[1]) ≈ -9.0
@test calc(parse_one("-((0.0-9.0))", all)[1]) ≈ 9.0
@test calc(parse_one("(-6.0/5.0)", all)[1]) ≈ -1.2
@test calc(parse_one("3.0*-((0.0-9.0))", all)[1]) ≈ 27
@test calc(parse_one("-9.0*3.0*-((0.0-9.0))*9.0", all)[1]) ≈ -2187.0
@test calc(parse_one("5.0/3.0*(-6.0/5.0)", all)[1]) ≈ -2.0
@test calc(parse_one("3*6-9*1", all)[1]) ≈ 9
@test calc(parse_one("4.0-5.0-0.0/8.0/5.0/3.0*(-6.0/5.0)-9.0*3.0*-((0.0-9.0))*9.0", all)[1]) ≈ -2188.0
@test calc(parse_one("((-6.0/6.0+7.0))*((-1.0-3.0/5.0))+-(9.0)", all)[1]) ≈ -18.6
# this has a large numerical error (~1e-15) and i don't understand why
@test abs(calc(parse_one("7.0/3.0*9.0-5.0-0.0+-(-9.0/7.0)*9.0*-0.0-7.0+-4.0-5.0", all)[1]) - 0.0) < 1e-10

x = Neg(Prd(Any[7.0,
                Inv(0.0),
                Inv(2.0),
                Inv(Neg(0.0))]))
y = calc(x)
z = -7.0/0.0/2.0/-0.0
println("$x $y $z")
@test isequal(y, z)

for x in [Inv(Neg(0.0)),
          Inv(Prd(Any[Neg(Sum(Any[Prd(Any[0.0])]))])),
          Inv(Sum(Any[Prd(Any[Neg(Sum(Any[Prd(Any[0.0])]))])]))
          ]
    y = calc(x)
    println("$x $y")
    @test isequal(y, -Inf)
end

x = Neg(Prd(Any[7.0,
                Inv(0.0),
                Inv(2.0),
                Inv(Sum(Any[Prd(Any[Neg(Sum(Any[Prd(Any[0.0])]))])])),
                3.0]))
y = calc(x)
z = -7.0/0.0/2.0/(-(0.0))*3.0
println("$x $y $z")
@test isequal(y, z)

p = Sum(Any[Prd(Any[-9.0]),
            Neg(Prd(Any[7.0,Inv(0.0),Inv(2.0),Inv(Sum(Any[Prd(Any[Neg(Sum(Any[Prd(Any[0.0])]))])])),3.0])),
            Neg(Prd(Any[7.0,Inv(Neg(Sum(Any[Prd(Any[9.0]),Prd(Any[5.0])])))])),
            Prd(Any[5.0]),
            Neg(Prd(Any[7.0]))])

a = eval(Meta.parse("-9.0-7.0/0.0/2.0/(-(0.0))*3.0-7.0/-(9.0+5.0)+5.0-7.0"))
b = parse_one("-9.0-7.0/0.0/2.0/(-(0.0))*3.0-7.0/-(9.0+5.0)+5.0-7.0", all)[1]
c = calc(b)
println("$a $b $c")
@test isequal(a, c)

# generate random expressions, parse them, and compare the results to
# evaluating in julia

function subexpr()
    expr = string(rand(0:9)) * ".0"
    for k in 1:rand(1:5)
        case = rand(1:6)
        if case < 5
            expr = expr * "+-*/"[case:case] * string(rand(0:9)) * ".0"
        elseif case == 5
            expr = "(" * expr * ")"
        elseif case == 6
            expr = "-" * expr
        end
    end
    expr
end

for i in 1:20
    expr = subexpr()
    for _ in 1:rand(0:4)
        expr = expr * string("+-*/"[rand(1:4)]) * subexpr()
    end
    # julia desn't like repeated -
    expr = replace(expr, r"-+" => "-")
    println("expr $(expr)")
    try
        a = eval(Meta.parse(expr))
        b = calc(parse_one(expr, all)[1])
        println("$a $b")
        if ! isequal(a, b)   # allow for Inf etc
            @test a ≈ b
        end
    catch
        @test_throws Exception parse(expr)
        @test_throws Exception calc(parse_one(expr, all)[1])
    end
end

println("calc ok")

end


text:
function default_sat(concept; encoding = :domain_wall)
    # Invalid binary representation
    check_validity(x, Val(encoding)) && return sqrt(length(X))
    # Return -1 if the configuration is a solution, 1 otherwise
    return concept(integerize(X)) ? -1. : 1.
end

default_opt(x...) = 0.0


text:
module Fox

using Playground
using ArgParse
using URIParser
import Logging
import YAML: load

PLAYGROUND_CONFIG_PATH = joinpath(homedir(), ".playground")
PLAYGROUND_CONFIG_FILE = joinpath(PLAYGROUND_CONFIG_PATH, "config.yml")

function run(cmd_args=ARGS)
    config = load_config(PLAYGROUND_CONFIG_FILE, PLAYGROUND_CONFIG_PATH)
    parse_settings = ArgParseSettings(suppress_warnings=true)

    @add_arg_table parse_settings begin
        "--julia-versions", "-j"
            help = "A list of julia versions to test on. Must available via Playgournd.jl"
            action = :store_arg
            nargs = '*'
            default = []
        "--config", "-c"
            help = "By default the fox.yml or .travis.yml in the local directory will be used."
            action = :store_arg
            default = ""
    end

    args = parse_args(cmd_args, parse_settings)
    config_file = args["config"]
    if config_file == ""
        root_files = readdir(pwd())

        if "fox.yml" in root_files
            config_file = joinpath(pwd(), "fox.yml")
        elseif ".travis.yml" in root_files
            config_file = joinpath(pwd(), ".travis.yml")
        else
            error("$(pwd()) doesn't contain a fox.yml or .travis.yml file. Please provide one.")
        end
    elseif !ispath(config_file)
        error("$config_file is not a valid path")
    end

    settings = load(open(config_file))

    julia_versions = settings["julia"]
    if !isempty(args["julia-versions"])
        julia_versions = intersect(
            julia_versions,
            args["julia-versions"]
        )
    end

    found = AbstractString[]
    for j in julia_versions
        Logging.debug(j)
        found = ""

        for v in readdir(config.dir.bin)
            # Split on "-" cause Playground.jl suggests labelling
            # julia versions as julia-0.4 rather than just 0.4
            Logging.debug(v)
            extracted = split(v, "-")
            Logging.debug(extracted)

            if v == j || (length(extracted) == 2 && extracted[2] == string(j))
                found = v
            end
        end

        if found == ""
            Logging.info("$j not installed. Skipping.")
        else
            # Process tests for julia version
            # Create the playground
            playground_path = joinpath(pwd(), ".fox", found)
            pg_config = Playground.PlaygroundConfig(config, playground_path, "")

            # Make sure the default_shell is set to
            pg_config.default_shell = "/bin/sh"

            # if ispath(playground_path)
            #     Logging.info("Cleaning existing playground $found ...")
            #     rm(config; dir=playground_path)
            # end

            Logging.info("Building test playground $found ...")
            create(
                config;
                dir=playground_path,
                julia=found
            )

            # Start setting up the current package in the playground
            # by 1) making a symlink 2) adding it to the REQUIRE file
            # and 3) calling Pkg.resolve()

            # Get current package name
            pkg_name = replace(basename(pwd()), ".jl", "")

            # Get the playground package version path
            pkg_path = ""
            for vpath in readdir(pg_config.pkg_path)
                if contains(found, strip(vpath, 'v'))
                    pkg_path = joinpath(pg_config.pkg_path, vpath)
                end
            end

            # Update the REQUIRE file
            if pkg_path != ""
                Playground.mklink(pwd(), joinpath(pkg_path, pkg_name))
                open(joinpath(pkg_path, "REQUIRE"), "w+") do fstream
                    existing = readall(fstream)
                    if !contains(existing, pkg_name)
                        write(fstream, "$pkg_name\n")
                    end
                end

                # Call resolve
                execute(config, `julia -e 'Pkg.resolve()'`, dir=playground_path)
            end

            # We'll need to update `execute` in Playground
            # to return the status code.
            for entry in settings["script"]
                cmd_str = entry

                # Handle the gross Pkg.clone conditions cause
                # duplicate runs will result in an error cause of cloning
                # a pkg that already exists.
                # This should probably be done in a more robust way...
                # Once we refactor in the future this should be its own function
                # and tested like crazy.
                if contains(cmd_str, "Pkg.clone(pwd())")
                    Logging.info("Ignoring clone of pwd as it has already been linked.")
                    cmd_str = replace(cmd_str, "Pkg.clone(pwd())", "")
                elseif contains(cmd_str, "Pkg.clone")
                    # Match against all 'Pkg.clone()'s
                    Logging.debug("Has 'Pkg.clone...' in it cmd_str")

                    for sub_str in matchall(r"Pkg.clone\(\"(.*?(?=\"))\"\)", cmd_str)
                        Logging.debug(sub_str)

                        # Extract the url
                        grps = match(r"Pkg.clone\(\"(.*?(?=\"))\"\)", sub_str).captures
                        Logging.debug(grps)

                        # Assuming we've extracted the url, extract the pkg name
                        # and add a check if the Pkg is available before running the line.
                        if grps != nothing && length(grps) == 1
                            remote_pkg_name = replace(basename(URI(grps[1]).path), ".jl", "")
                            new_clone = "try Pkg.available(\"$(remote_pkg_name)\") catch $sub_str end"

                            Logging.debug("Replacing $sub_str with $new_clone")
                            cmd_str = replace(
                                cmd_str,
                                sub_str,
                                new_clone
                            )
                        end
                    end
                end

                if !contains(cmd_str, "git fetch --unshallow")
                    execute(
                        config,
                        `$(Base.shell_split(cmd_str))`;
                        dir=playground_path
                    )
                end
            end
        end
    end
end

end # module


text:
#Import CppSimData or explain how to install.

try
eval(:(import CppSimData))
catch e
msg = "This sample requires data installed in the CppSimData module."
msg *= "\nTo continue demo, install with the following:\n\n"
msg *= "    Pkg.clone(\"git://github.com/ma-laforge/CppSimData.jl\")"
@info(msg)
println();println()
rethrow(e)
end

#Last line


text:
# Factor Graph OS type utilities
#  IIF methods should direclty detect extended types from user import
# of convert in their namespace

function convert{PT <: PackedInferenceType, T <:FunctorInferenceType}(::Type{PT}, ::T)
  getfield(T.name.module, Symbol("Packed$(T.name.name)"))
end
function convert{T <: FunctorInferenceType, PT <: PackedInferenceType}(::Type{T}, ::PT)
  getfield(PT.name.module, Symbol(string(PT.name.name)[7:end]))
end

"""
    convert2packedfunctionnode(fgl::FactorGraph, fsym::Symbol)

Encode complicated function node type with assumed to exist, user supplied convert
function to related 'Packed<type>' format.
"""
function convert2packedfunctionnode(fgl::FactorGraph,
      fsym::Symbol,
      api::DataLayerAPI=localapi  )
  #
  fid = fgl.fIDs[fsym]
  fnc = getfnctype(fgl, fid)
  usrtyp = convert(PackedInferenceType, fnc)
  cfnd = convert(PackedFunctionNodeData{usrtyp},getData(fgl, fid, api=api) )
  return cfnd, usrtyp
end

"""
    encodefg(fgl::FactorGraph)

Make a full memory copy of the graph and encode all complicated function node
types with assumed to exist convert to 'Packed<type>' formats. Same converters
as used for database persistence storage with CloudGraphs.jl.
"""
function encodefg(fgl::FactorGraph;
      api::DataLayerAPI=localapi  )
  #
  fgs = deepcopy(fgl)
  fgs.cg = nothing
  fgs.registeredModuleFunctions = nothing
  fgs.g = Graphs.incdict(Graphs.ExVertex,is_directed=false)
  @showprogress 1 "Encoding variables..." for (vsym,vid) in fgl.IDs
    cpvert = deepcopy(getVert(fgl, vid, api=api))
    api.addvertex!(fgs, cpvert) #, labels=vnlbls)  # currently losing labels
  end

  @showprogress 1 "Encoding factors..." for (fsym,fid) in fgs.fIDs
    data,ftyp = convert2packedfunctionnode(fgl, fsym)
    # data = FunctionNodeData{ftyp}(Int[], false, false, Int[], m, gwpf)
    newvert = ExVertex(fid,string(fsym))
    for (key,val) in getVert(fgl,fid,api=api).attributes
      newvert.attributes[key] = val
    end
    setData!(newvert, data)
    api.addvertex!(fgs, newvert)
  end
  fgs.g.inclist = typeof(fgl.g.inclist)()

  # iterated over all edges
  @showprogress 1 "Encoding edges..." for (eid, edges) in fgl.g.inclist
    fgs.g.inclist[eid] = Vector{typeof(edges[1])}()
    for ed in edges
      newed = Graphs.Edge(ed.index,
          fgs.g.vertices[ed.source.index],
          fgs.g.vertices[ed.target.index]  )
      push!(fgs.g.inclist[eid], newed)
    end
  end

  return fgs
end

"""
    savefgjld(fgl::FactorGraph; file::AbstractString="tempfg.jld")

Save mostly complete Factor Graph type by converting complicated FunctionNodeData
types to 'Packed' types using user supplied converters. Ground truth can also be
saved and recovered by the associated loadjld(file="tempfg.jld") method.
"""
function savejld(fgl::FactorGraph;
      file::AbstractString="tempfg.jld",
      groundtruth=nothing)
  fgs = encodefg(fgl)
  if groundtruth == nothing
    @save file fgs
  else
    @save file fgs groundtruth
  end
  return file
end

"""
    convertfrompackedfunctionnode(fgl, fsym)

If you get unknown type conversion error when loading a .jld, while using your own
FunctorInferenceTypes, you should:
Copy these functions below, and overload in your package with local extented
FunctorInferenceType definitions.
See RoME/src/fgos.jl for example.
"""
function convertfrompackedfunctionnode(fgl::FactorGraph,
      fsym::Symbol,
      api::DataLayerAPI=localapi  )
  #
  fid = fgl.fIDs[fsym]
  fnc = getData(fgl, fid).fnc #getfnctype(fgl, fid)
  usrtyp = convert(FunctorInferenceType, fnc)
  data = getData(fgl, fid, api=api)
  newtype = FunctionNodeData{GenericWrapParam{usrtyp}}
  cfnd = convert(newtype, data)
  return cfnd, usrtyp
end

"""
    decodefg(fgs::FactorGraph)

Unpack PackedFunctionNodeData formats back to regular FunctonNodeData.
"""
function decodefg(fgs::FactorGraph; api::DataLayerAPI=localapi)
  fgu = deepcopy(fgs)
  fgu.cg = nothing
  fgu.registeredModuleFunctions = nothing
  fgu.g = Graphs.incdict(Graphs.ExVertex,is_directed=false)
  @showprogress 1 "Decoding variables..." for (vsym,vid) in fgs.IDs
    cpvert = deepcopy(getVert(fgs, vid, api=api))
    api.addvertex!(fgu, cpvert) #, labels=vnlbls)  # currently losing labels
  end

  @showprogress 1 "Decoding factors..." for (fsym,fid) in fgu.fIDs
    data,ftyp = convertfrompackedfunctionnode(fgs, fsym)
    # data = FunctionNodeData{ftyp}(Int[], false, false, Int[], m, gwpf)
    newvert = ExVertex(fid,string(fsym))
    for (key,val) in getVert(fgs,fid,api=api).attributes
      newvert.attributes[key] = val
    end
    setData!(newvert, data)
    api.addvertex!(fgu, newvert)
  end
  fgu.g.inclist = typeof(fgs.g.inclist)()

  # iterated over all edges
  @showprogress 1 "Decoding edges..." for (eid, edges) in fgs.g.inclist
    fgu.g.inclist[eid] = Vector{typeof(edges[1])}()
    for ed in edges
      newed = Graphs.Edge(ed.index,
          fgu.g.vertices[ed.source.index],
          fgu.g.vertices[ed.target.index]  )
      push!(fgu.g.inclist[eid], newed)
    end
  end

  return fgu
end

"""
    loadjld(file="tempfg.jld")

Opposite of savejld(fg, gt=gt, file="tempfg.jl") to load data from file. This function
uses the unpacking converters for converting all PackedInferenceType to FunctorInferenceType.
"""
function loadjld(;file::AbstractString="tempfg.jld")
  dd = JLD.load(file)
  # fgs = jldopen(file,"r") do file
  #   read(file, "fgs")
  # end
  fgs, gt = nothing, nothing
  if haskey(dd, "fgs")
    fgs = dd["fgs"]
  else
    error("No factor graph (fgs) data found in this file, only found $(keys(dd))")
  end
  if haskey(dd, "groundtruth")
    gt = dd["groundtruth"]
    println("Also found ground truth data")
  end
  fgd = decodefg(fgs)
  return fgd, gt
end

function ls(fgl::FactorGraph, lbl::Symbol; api::DataLayerAPI=dlapi, ring::Int=1)
  # TODO ring functionality must still be implemented
  lsa = Symbol[]
  # v = nothing
  if haskey(fgl.IDs, lbl)
    id = fgl.IDs[lbl]
  else
    return lsa
  end
  # this is unnecessary
  v = getVert(fgl,id, api=api)
  for outn in api.outneighbors(fgl, v)
    # if outn.attributes["ready"] = 1 && outn.attributes["backendset"]=1
      push!(lsa, Symbol(outn.label))
    # end
  end
  return lsa
end
ls(fgl::FactorGraph, lbl::T) where {T <: AbstractString} = ls(fgl, Symbol(lbl))

function ls(fgl::FactorGraph)
  k = collect(keys(fgl.IDs))
  l = Int[]
  x = Int[]
  for kk in k
    kstr = string(kk)
    val = parse(Int,kstr[2:end]) # kk
    if kstr[1] == 'l'
      push!(l,val)
    elseif kstr[1] == 'x'
      push!(x,val)
    end
  end
  l = sort(l)
  x = sort(x)
  ll = Array{Symbol,1}(length(l))
  xx = Array{Symbol,1}(length(x))
  for i in 1:length(l)
    ll[i] = Symbol(string("l",l[i]))
  end
  for i in 1:length(x)
    xx[i] = Symbol(string("x",x[i]))
  end
  return xx,ll
end

function lsf(fgl::FactorGraph, lbl::Symbol; api::DataLayerAPI=dlapi)
  lsa = Symbol[]
  # v = Union{}
  if haskey(fgl.fIDs, lbl)
    id = fgl.fIDs[lbl]
  else
    return lsa
  end
  v = getVert(fgl, id, api=api) # fgl.g.vertices[id] #fgl.f[id]
  for outn in api.outneighbors(fgl, v) # out_neighbors(v, fgl.g)
    push!(lsa, Symbol(outn.label))
  end
  return lsa
end
lsf{T <: AbstractString}(fgl::FactorGraph, lbl::T) = lsf(fgl,Symbol(lbl))

function lsf{T <: FunctorInferenceType}(fgl::FactorGraph,
      mt::Type{T};
      api::DataLayerAPI=dlapi  )
  #
  syms = Symbol[]
  for (fsym,fid) in fgl.fIDs
    if typeof(getfnctype(fgl, fid, api=api))==T
      push!(syms, fsym)
    end
  end
  return syms
end

function ls2(fgl::FactorGraph, vsym::Symbol)
  xxf = ls(fgl, vsym)
  xlxl = Symbol[]
  for xf in xxf
    xx = lsf(fgl,xf)
    xlxl = union(xlxl, xx)
  end
  xlxl = setdiff(xlxl, [vsym])
  return xlxl
end

"""
    landmarks(fgl::FactorGraph, vsym::Symbol)

Return Vector{Symbol} of landmarks attached to vertex vsym in fgl.
"""
function landmarks(fgl::FactorGraph, vsym::Symbol)
  fsyms = ls(fgl, vsym)
  lms = Symbol[]
  for fs in fsyms
    for varv = lsf(fgl, fs)
      if string(varv)[1] == 'l'
        push!(lms, varv)
      end
    end
  end
  lms
end

#


text:
using Test
using LinearAlgebra

using CholeskyDecomposition

@testset "Cholesky decomposition" begin
    for i in 1:5
        n = rand(1:25)
        # Generates a random hermitian positive definite matrix
        A = begin
            L = rand(Complex{Float64}, n, n)

            # Makes it diagonally dominant
            (L + L') / sqrt(8.0) + I(n) * n
        end

        @assert ishermitian(A) && isposdef(A)

        @test begin
            cholesky_outer_product_decompose(A) ≈
            cholesky_banachiewicz_decompose(A) ≈
            cholesky_crout_decompose(A) ≈
            cholesky(A).L
        end
    end
end

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using LightGraphs.SimpleGraphs

function isvalid(g::SimpleGraph, colors)
    degrees = degree(g)
    usedc_nodes = [Vector{Int}(undef, d) for d in degrees]
    empty!.(usedc_nodes)

    for (c, e) in zip(colors, edges(g))
        usedc_src = usedc_nodes[e.src]
        usedc_dst = usedc_nodes[e.dst]
        if c ∈ usedc_src || c ∈ usedc_dst
            return false
        else
            push!(usedc_src, c)
            push!(usedc_dst, c)
        end
    end
    return true
end

function color_edges_greedy(g::SimpleGraph)
    degrees = degree(g)
    colors = zeros(Int, ne(g))
    usedc_nodes = [Vector{Int}(undef, d) for d in degrees]
    empty!.(usedc_nodes)

    maxd = maximum(degrees)
    cstat = zeros(Int, maxd)
    corder = collect(1:maxd+1)

    for (i, e) in enumerate(edges(g))
        usedc_src = usedc_nodes[e.src]
        usedc_dst = usedc_nodes[e.dst]

        # somewhat arbitrary: for big exclude lists it is beneficial to sort them first
        sortfirst = length(usedc_src) + length(usedc_dst) > 1000
        c = pickfrom_without(corder, usedc_src, usedc_dst; sortfirst)

        colors[i] = c
        push!(usedc_src, c)
        push!(usedc_dst, c)

        if c > length(cstat)
            @assert c == length(cstat) + 1
            resize!(cstat, c)
            cstat[c] = 0
            resize!(corder, c+1)
            corder[c+1] = c+1
        end
        cstat[c] += 1
        sortperm!(view(corder, 1:length(cstat)), cstat)
    end
    @assert isvalid(g, colors) "Could not find valid edge coloring 😱"
    return colors
end

@inline function pickfrom_without(X, exA, exB; sortfirst=false)
    if sortfirst
        ex = sort!(vcat(exA, exB))
        for x in X
            insorted(x, ex) || return x
        end
    else
        for x in X
            if x ∉ exA && x ∉ exB
                return x
            end
        end
    end
    error("did not find any element")
end


text:
using Switching
using Test

@testset "Switching.jl" begin
    # Write your own tests here.
end


text:
# Linear Regression

using CSV
using Statistics
using Missings
using Plots
using StatsBase
using HypothesisTests

# Read data
dataset = CSV.read("C:/Users/Andrew.Bannerman/Desktop/Machine Learning/Machine-Learning/Machine-Learning/linear_regression/Salary_Data.csv", header=true, types=[Float64, Float64])
dataset[:index] = collect(1:1:size(dataset,1)) # index the data
dataset = dataset[sample(axes(dataset, 1), size(dataset,1); replace = false, ordered = false), :] # Random shuffle the dataframe

# Feature Scaling
dataset[:Salary] = standardization(Float64.(dataset[:Salary]))
dataset[:YearsExperience] = standardization(Float64.(dataset[:YearsExperience]))

# Split data to train and test sets
split_ratio = 2/3
len = size(dataset,1)
split_no = Int64.(round(split_ratio * len))
train_set = dataset[1:split_no,:]
test_set = dataset[split_no+1:len,:]

# Visualize the data
gr() # We will continue onward using the GR backend
plot(train_set[:YearsExperience],train_set[:Salary],seriestype=:scatter,title="My Scatter Plot (Training Data)")

# Find the regression formula
x = Float64.(train_set[:YearsExperience])
y = Float64.(train_set[:Salary])

# Regression coefficients (y = b0 .+ (b1 .* x))
b1 = cov(x,y) / var(x) # b (slope)
b0 = mean(y) - b1 * mean(x) # a (intercept)

# Apply the regression formula to the training data
y_hat_train = fill(0.0,size(train_set,1)) # initialize array
y_hat_train = b0 .+ (b1 .* x)
plot!(train_set[:YearsExperience],y_hat_train, smooth=true)
plot_annotation = "y = ",round(b0,digits=2)," + (",round(b1,digits=2)," * x)"
annotate!(-.5, 1.0, text(plot_annotation, :red, :left, 10))
training_set_residuals = y .- y_hat_train
r_2 = cor(x,y) ^2
n = size(train_set,1)
standard_error = sqrt((sum(y_hat_train .- y)^2) / (n-2))

# Using the same regression coefficients found on training data - plot the regression line on unseen test data
x = Float64.(test_set[:YearsExperience])
y = Float64.(test_set[:Salary])
y_hat_test = fill(0.0,size(test_set,1))
y_hat_test = b0 .+ (b1 .* x)
plot(test_set[:YearsExperience],test_set[:Salary],seriestype=:scatter,title="My Scatter Plot (Test Data)")
plot!(test_set[:YearsExperience],y_hat_test, smooth=true)
plot_annotation = "y = ",round(b0,digits=2)," + (",round(b1,digits=2)," * x)"
annotate!(-.5, 1.0, text(plot_annotation, :red, :left, 10))
test_set_residuals = y .- y_hat_test
r_2 = cor(x,y) ^2
n = size(test_set,1)
standard_error = sqrt((sum(y_hat_test .- y)^2) / (n-2))

"""
```
lin_reg_ols(x_train::Array{T},y_train::Array{T},x_test::Array{T},y_test::Array{T})
```
Train the regression model on training data - derive coefficients b0 and b1
Apply training coefficients on test set data:: y_hat_test = b0 .+ (b1 .* x_test)
"""
function lin_reg_ols(x_train::Array{T},y_train::Array{T},x_test::Array{T},y_test::Array{T})::Array{T} where {T<:Real}
    let y_hat_train = size(x_train,1), y_hat_test = size(x_test,1)
    b1 = cov(x_train,y_train) / var(x_train) # b (slope)
    b0 = mean(y_train) - b1 * mean(x_train) # a (intercept)
    coefficients = b1,b0
    y_hat_train = b0 .+ (b1 .* x_train)
    # Test set
    y_hat_test = b0 .+ (b1 .* x_test)
    return y_hat_test
end
end

# Test function
x_train = Float64.(train_set[:YearsExperience])
y_train= Float64.(train_set[:Salary])
x_test = Float64.(test_set[:YearsExperience])
y_test= Float64.(test_set[:Salary])

test_predictions = lin_reg_ols(x_train,y_train,x_test,y_test)
test_set_residuals = y_test .- test_predictions
r_2 = cor(x_test,y_test) ^2

# How to calculate the p-value
# https://www.wikihow.com/Calculate-P-Value

max_stars_repo_name: JuliaTagBot/Machine-Learning

max_issues_repo_name: JuliaTagBot/Machine-Learning

max_forks_repo_name: JuliaTagBot/Machine-Learning


text:
module OpenModels
import Base: ==
using SemanticModels.ModelTools
import SemanticModels.ModelTools: model
using Catlab.Graphics.Graphviz
import Catlab.Graphics.Graphviz: Graph, Edge

export OpenModel

struct OpenModel{V,M} <: AbstractModel
    dom::V
    model::M
    codom::V
end

"""    Graph(f::OpenModel)

convert an OpenModel into a GraphViz Graph. This calls Graph(::Model) and then adds vertices (and edges) for the domain and codomain of the open model.
"""
function Graph(f::OpenModel)
    g = Graph(f.model)
    A, M, B = f.dom, f.model, f.codom
    stmts_dom = map(enumerate(A)) do (i,a)
        m = M.S[a]
        Edge(["I$i", "$m"], Attributes(:style=>"dashed"))
    end
    stmts_codom = map(enumerate(B)) do (i,a)
        m = M.S[a]
        Edge(["$m", "O$i"], Attributes(:style=>"dashed"))
    end
    append!(g.stmts, append!(stmts_dom, stmts_codom))
    return g
end

function ==(f::OpenModel,g::OpenModel)
    all(isequal.(f.dom, g.dom)) && all(isequal.(f.codom, g.codom)) && f.model == g.model
end

# otimes and compose general that use disjoin union functions

end


text:
using StaticArrays, Test
using OrdinaryDiffEq
using RecursiveArrayTools

u0 = [fill(2, MVector{2,Float64}), ones(MVector{2,Float64})]
f = (du,u,p,t) -> du .= u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = [fill(2, SVector{2,Float64}), ones(SVector{2,Float64})]
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, SSPRK22(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = ones(MVector{2,Float64})
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = ones(SVector{2,Float64})
f = (u,p,t) -> u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
integrator = init(ode, ImplicitEuler())
@test OrdinaryDiffEq.get_W(integrator.cache.nlsolver) isa StaticArrays.LU
sol = solve(ode, ImplicitEuler())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, ImplicitEuler(nlsolve=NLAnderson()))
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

#https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/373
function lorenz_static(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 @SVector [dx,dy,dz]
end

u0 = @SVector [1.0,0.0,0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz_static,u0,tspan)
solve(prob,dt=0.1,Rosenbrock23(autodiff=false))

# Check that ArrayPartitions of static vectors work
#https://github.com/SciML/OrdinaryDiffEq.jl/issues/1308
function lorenz_static(u::ArrayPartition,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 du1 = @SVector [dx,dy]
 du2 = @SVector [dz]
 ArrayPartition(du1,du2)
end

u01 = @SVector [1.0,0.0]
u02 = @SVector [0.0]
u0ap = ArrayPartition(u01,u02)
probap = ODEProblem(lorenz_static,u0ap,tspan)

sol = solve(prob,dt=1e-2,Heun())
solap = solve(probap,dt=1e-2,Heun())
@test sol(30) == solap(30)

sol = solve(prob,dt=1e-2,Tsit5())
solap = solve(probap,dt=1e-2,Tsit5())
@test sol(30) == solap(30)