text:
using DocStringExtensions
"""
Runs the model for several precipitation zones of a catchment

$(SIGNATURES)

The function returns the sum of the discharges of the precipitation zones as well as the mean difference in snow coverage compared to observed data.
As input potential evaporation, precipitation and temperature data is necessary.
- Inputs, Storages and Paramaters for each elevation zone
- Area_Zones : Area of all precipitation zones in m²
- Elevation_Percentage: areal extend of each elevation zone in each precipitation zone
- Elevation_Zone_Catchment
- ID_Prec_Zones: IDs of rain gauges of each precipitation zone
- Nr_Elevationbands_All_Zones
- observed_snow_cover: the observed snow cover of each elevation zone in each precipitation zone
- year_snow_observations: the day of the timeseries when the first observed snow cover is assessed
"""
function runmodelprecipitationzones(Potential_Evaporation::Array{Float64,1}, Precipitation_All_Zones::Array{Array{Float64,2},1}, Temperature_Elevation_Catchment::Array{Float64,2}, Inputs_All_Zones::Array{Array{HRU_Input,1},1}, Storages_All_Zones::Array{Array{Storages,1},1}, SlowStorage::Float64, parameters::Array{Parameters,1}, slow_parameters::Slow_Paramters, Area_Zones::Array{Float64,1}, Area_Zones_Percent::Array{Float64,1}, Elevation_Percentage::Array{Array{Float64,1},1}, Elevation_Zone_Catchment::Array{Float64,1}, ID_Prec_Zones::Array{Int64,1}, Nr_Elevationbands_All_Zones::Array{Int64,1}, observed_snow_cover::Array{Array{Float64,2},1}, year_snow_observations::Int64)
        Total_Discharge = zeros(length(Precipitation_All_Zones[1][:,1]))
        count = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        Snow_Overall_Objective_Function = 0
        for i in 1: length(ID_Prec_Zones)
                # take the storages and input of the specific precipitation zone
                Inputs_HRUs = Inputs_All_Zones[i]
                Storages_HRUs = Storages_All_Zones[i]
                # run the model for the specific precipitation zone
                Discharge, Snow_Extend, Waterbalance = run_model(Area_Zones[i], Potential_Evaporation, Precipitation_All_Zones[i], Temperature_Elevation_Catchment,
                        Inputs_HRUs[1], Inputs_HRUs[2], Inputs_HRUs[3], Inputs_HRUs[4],
                        Storages_HRUs[1], Storages_HRUs[2], Storages_HRUs[3], Storages_HRUs[4], SlowStorage,
                        parameters[1], parameters[2], parameters[3], parameters[4], slow_parameters, Nr_Elevationbands_All_Zones[i], Elevation_Percentage[i])
                # sum up the discharge of all precipitation zones
                @assert -10^(-10) <= Waterbalance <= 10^(-10)
                Total_Discharge += Discharge
                #snow extend is given as 0 or 1 for each elevation zone at each timestep)
                elevations = size(observed_snow_cover[i],2)
                # only use the modeled snow cover data that is in line with the observed snow cover data
                snow_cover_modelled = Snow_Extend[year_snow_observations: year_snow_observations + length(observed_snow_cover[i][:,1]) - 1, :]
                Mean_difference = 0
                #calculate the mean difference for all elevation zones
                for h in 1: elevations
                        Difference = snowcover(snow_cover_modelled[:,h], observed_snow_cover[i][:,h])
                        # take the area weighted average mean difference in snow cover
                        Mean_difference += Difference * Elevation_Percentage[i][h]
                end
                #take the area weighted mean difference in snow cover
                Snow_Overall_Objective_Function += Mean_difference * Area_Zones_Percent[i]
        end
        # calculate the mean difference over all precipitation zones
        #Snow_Overall_Objective_Function = Snow_Overall_Objective_Function / length(ID_Prec_Zones)
        return Total_Discharge::Array{Float64,1}, Snow_Overall_Objective_Function::Float64
end

"""
Runs the model for several precipitation zones of a catchment

$(SIGNATURES)

The function returns the sum of the discharges of the precipitation zones as well as the mean difference in snow coverage compared to observed data.
As input potential evaporation, precipitation and temperature data is necessary.
- Inputs, Storages and Paramaters for each elevation zone
- Area_Zones : Area of all precipitation zones in m²
- Elevation_Percentage: areal extend of each elevation zone in each precipitation zone
- Elevation_Zone_Catchment
- ID_Prec_Zones: IDs of rain gauges of each precipitation zone
- Nr_Elevationbands_All_Zones
- observed_snow_cover: the observed snow cover of each elevation zone in each precipitation zone
- year_snow_observations: the day of the timeseries when the first observed snow cover is assessed
"""
function runmodelprecipitationzones_glacier(Potential_Evaporation::Array{Float64,1}, Glacier::Array{Array{Float64,2},1}, Precipitation_All_Zones::Array{Array{Float64,2},1}, Temperature_Elevation_Catchment::Array{Float64,2}, Inputs_All_Zones::Array{Array{HRU_Input,1},1}, Storages_All_Zones::Array{Array{Storages,1},1}, SlowStorage::Float64, parameters::Array{Parameters,1}, slow_parameters::Slow_Paramters, Area_Zones::Array{Float64,1}, Area_Zones_Percent::Array{Float64,1}, Elevation_Percentage::Array{Array{Float64,1},1}, Elevation_Zone_Catchment::Array{Float64,1}, ID_Prec_Zones::Array{Int64,1}, Nr_Elevationbands_All_Zones::Array{Int64,1}, observed_snow_cover::Array{Array{Float64,2},1}, year_snow_observations::Int64)
        Total_Discharge = zeros(length(Precipitation_All_Zones[1][:,1]))
        count = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        Snow_Overall_Objective_Function = 0
        for i in 1: length(ID_Prec_Zones)
                # take the storages and input of the specific precipitation zone
                Inputs_HRUs = Inputs_All_Zones[i]
                Storages_HRUs = Storages_All_Zones[i]
                # run the model for the specific precipitation zone
                Discharge, Snow_Extend, Waterbalance = run_model_glacier(Area_Zones[i], Potential_Evaporation, Glacier[i], Precipitation_All_Zones[i], Temperature_Elevation_Catchment,
                        Inputs_HRUs[1], Inputs_HRUs[2], Inputs_HRUs[3], Inputs_HRUs[4],
                        Storages_HRUs[1], Storages_HRUs[2], Storages_HRUs[3], Storages_HRUs[4], SlowStorage,
                        parameters[1], parameters[2], parameters[3], parameters[4], slow_parameters, Nr_Elevationbands_All_Zones[i], Elevation_Percentage[i])
                # sum up the discharge of all precipitation zones
                #@assert -10^(-10) <= Waterbalance <= 10^(-10)
                Total_Discharge += Discharge
                #snow extend is given as 0 or 1 for each elevation zone at each timestep)
                elevations = size(observed_snow_cover[i],2)
                # only use the modeled snow cover data that is in line with the observed snow cover data
                snow_cover_modelled = Snow_Extend[year_snow_observations: year_snow_observations + length(observed_snow_cover[i][:,1]) - 1, :]
                Mean_difference = 0
                #calculate the mean difference for all elevation zones
                for h in 1: elevations
                        Difference = snowcover(snow_cover_modelled[:,h], observed_snow_cover[i][:,h])
                        # take the area weighted average mean difference in snow cover
                        Mean_difference += Difference * Elevation_Percentage[i][h]
                end
                #take the area weighted mean difference in snow cover
                Snow_Overall_Objective_Function += Mean_difference * Area_Zones_Percent[i]
        end
        # calculate the mean difference over all precipitation zones
        #Snow_Overall_Objective_Function = Snow_Overall_Objective_Function / length(ID_Prec_Zones)
        return Total_Discharge::Array{Float64,1}, Snow_Overall_Objective_Function::Float64
end

function runmodelprecipitationzones_future(Potential_Evaporation::Array{Float64,1}, Precipitation_All_Zones::Array{Array{Float64,2},1}, Temperature_Elevation_Catchment::Array{Float64,2}, Inputs_All_Zones::Array{Array{HRU_Input,1},1}, Storages_All_Zones::Array{Array{Storages,1},1}, SlowStorage::Float64, parameters::Array{Parameters,1}, slow_parameters::Slow_Paramters, Area_Zones::Array{Float64,1}, Area_Zones_Percent::Array{Float64,1}, Elevation_Percentage::Array{Array{Float64,1},1}, Elevation_Zone_Catchment::Array{Float64,1}, ID_Prec_Zones::Array{Int64,1}, Nr_Elevationbands_All_Zones::Array{Int64,1}, Elevations_Each_Precipitation_Zone::Array{Array{Float64,1},1})
        Total_Discharge = zeros(length(Precipitation_All_Zones[1][:,1]))
        Total_Snow_Melt = zeros(length(Precipitation_All_Zones[1][:,1]))
        count = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        Snow_Overall_Objective_Function = 0
        snow_cover_modelled_allprec =  zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        Denominator_All_Prec = zeros(length(Elevation_Zone_Catchment))
        for i in 1: length(ID_Prec_Zones)
                # take the storages and input of the specific precipitation zone
                Inputs_HRUs = Inputs_All_Zones[i]
                Storages_HRUs = Storages_All_Zones[i]
                # run the model for the specific precipitation zone
                Discharge, Snow_Extend, Snow_Storage, Waterbalance, Snow_Melt = run_model_future(Area_Zones[i], Potential_Evaporation, Precipitation_All_Zones[i], Temperature_Elevation_Catchment,
                        Inputs_HRUs[1], Inputs_HRUs[2], Inputs_HRUs[3], Inputs_HRUs[4],
                        Storages_HRUs[1], Storages_HRUs[2], Storages_HRUs[3], Storages_HRUs[4], SlowStorage,
                        parameters[1], parameters[2], parameters[3], parameters[4], slow_parameters, Nr_Elevationbands_All_Zones[i], Elevation_Percentage[i])
                # sum up the discharge of all precipitation zones
                Total_Discharge += Discharge
                Total_Snow_Melt += Snow_Melt * Area_Zones_Percent[i]
                snow_cover_modelled = Snow_Extend
                Denominator_All = Float64[]
                for h in 1: length(snow_cover_modelled[1,:])
                        snow_cover_modelled[:,h] = snow_cover_modelled[:,h] * Elevation_Percentage[i][h] * Area_Zones_Percent[i]
                        Denominator = Elevation_Percentage[i][h] * Area_Zones_Percent[i]
                        #push!(Nominator_All, Nominator)
                        append!(Denominator_All, Denominator)
                end
                # make snow cover data same length for all prec zones, therefore min and max elevation of prec zone needs to be known
                snow_cover_modelled_matrix = zeros((length(Discharge), length(Elevation_Zone_Catchment)))
                Denominator_All_matrix = zeros(length(Elevation_Zone_Catchment))
                #print(Elevation_Zone_Catchment)
                #print(Elevations_Each_Precipitation_Zone[i], "\n")
                for (j,current_elevation) in enumerate(Elevations_Each_Precipitation_Zone[i])
                        #print(((current_elevation - Elevation_Zone_Catchment[1]) / 200)+1, "\n")
                        snow_cover_modelled_matrix[:,Int((current_elevation - Elevation_Zone_Catchment[1]) / 200)+1] = snow_cover_modelled[:,j]
                        Denominator_All_matrix[Int((current_elevation - Elevation_Zone_Catchment[1]) / 200)+1] = Denominator_All[j]

                end
                snow_cover_modelled_allprec += snow_cover_modelled_matrix
                Denominator_All_Prec += Denominator_All_matrix
                # #take the area weighted mean difference in snow cover
                # Snow_Overall_Objective_Function += Mean_difference * Area_Zones_Percent[i]

        end
        #print(size(snow_cover_modelled_allprec), size(Denominator_All_Prec))
        Snow_Cover = transpose(snow_cover_modelled_allprec) ./ Denominator_All_Prec
        # calculate the mean difference over all precipitation zones
        #Snow_Overall_Objective_Function = Snow_Overall_Objective_Function / length(ID_Prec_Zones)
        return Total_Discharge::Array{Float64,1}, Snow_Cover::Array{Float64,2}, Total_Snow_Melt::Array{Float64,1}
end

function runmodelprecipitationzones_glacier_future(Potential_Evaporation::Array{Float64,1}, Glacier::Array{Array{Float64,2},1}, Precipitation_All_Zones::Array{Array{Float64,2},1}, Temperature_Elevation_Catchment::Array{Float64,2}, Inputs_All_Zones::Array{Array{HRU_Input,1},1}, Storages_All_Zones::Array{Array{Storages,1},1}, SlowStorage::Float64, parameters::Array{Parameters,1}, slow_parameters::Slow_Paramters, Area_Zones::Array{Float64,1}, Area_Zones_Percent::Array{Float64,1}, Elevation_Percentage::Array{Array{Float64,1},1}, Elevation_Zone_Catchment::Array{Float64,1}, ID_Prec_Zones::Array{Int64,1}, Nr_Elevationbands_All_Zones::Array{Int64,1}, Elevations_Each_Precipitation_Zone::Array{Array{Float64,1},1})
        Total_Discharge = zeros(length(Precipitation_All_Zones[1][:,1]))
        Total_Snow_Melt = zeros(length(Precipitation_All_Zones[1][:,1]))
        count = zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        Snow_Overall_Objective_Function = 0
        snow_cover_modelled_allprec =  zeros(length(Precipitation_All_Zones[1][:,1]), length(Elevation_Zone_Catchment))
        Denominator_All_Prec = zeros(length(Elevation_Zone_Catchment))
        for i in 1: length(ID_Prec_Zones)
                # take the storages and input of the specific precipitation zone
                Inputs_HRUs = Inputs_All_Zones[i]
                Storages_HRUs = Storages_All_Zones[i]
                # run the model for the specific precipitation zone
                Discharge, Snow_Extend, Snow_Storage, Waterbalance, Snow_Melt = run_model_glacier_future(Area_Zones[i], Potential_Evaporation, Glacier[i], Precipitation_All_Zones[i], Temperature_Elevation_Catchment,
                        Inputs_HRUs[1], Inputs_HRUs[2], Inputs_HRUs[3], Inputs_HRUs[4],
                        Storages_HRUs[1], Storages_HRUs[2], Storages_HRUs[3], Storages_HRUs[4], SlowStorage,
                        parameters[1], parameters[2], parameters[3], parameters[4], slow_parameters, Nr_Elevationbands_All_Zones[i], Elevation_Percentage[i])
                # sum up the discharge of all precipitation zones
                Total_Discharge += Discharge
                Total_Snow_Melt += Snow_Melt * Area_Zones_Percent[i]
                snow_cover_modelled = Snow_Extend
                Denominator_All = Float64[]
                for h in 1: length(snow_cover_modelled[1,:])
                        snow_cover_modelled[:,h] = snow_cover_modelled[:,h] * Elevation_Percentage[i][h] * Area_Zones_Percent[i]
                        Denominator = Elevation_Percentage[i][h] * Area_Zones_Percent[i]
                        #push!(Nominator_All, Nominator)
                        append!(Denominator_All, Denominator)
                end
                # make snow cover data same length for all prec zones, therefore min and max elevation of prec zone needs to be known
                snow_cover_modelled_matrix = zeros((length(Discharge), length(Elevation_Zone_Catchment)))
                Denominator_All_matrix = zeros(length(Elevation_Zone_Catchment))
                #print(Elevation_Zone_Catchment)
                #print(Elevations_Each_Precipitation_Zone[i], "\n")
                for (j,current_elevation) in enumerate(Elevations_Each_Precipitation_Zone[i])
                        #print(((current_elevation - Elevation_Zone_Catchment[1]) / 200)+1, "\n")
                        snow_cover_modelled_matrix[:,Int((current_elevation - Elevation_Zone_Catchment[1]) / 200)+1] = snow_cover_modelled[:,j]
                        Denominator_All_matrix[Int((current_elevation - Elevation_Zone_Catchment[1]) / 200)+1] = Denominator_All[j]

                end
                snow_cover_modelled_allprec += snow_cover_modelled_matrix
                Denominator_All_Prec += Denominator_All_matrix
                # #take the area weighted mean difference in snow cover
                # Snow_Overall_Objective_Function += Mean_difference * Area_Zones_Percent[i]

        end
        #print(size(snow_cover_modelled_allprec), size(Denominator_All_Prec))
        Snow_Cover = transpose(snow_cover_modelled_allprec) ./ Denominator_All_Prec
        # calculate the mean difference over all precipitation zones
        #Snow_Overall_Objective_Function = Snow_Overall_Objective_Function / length(ID_Prec_Zones)
        return Total_Discharge::Array{Float64,1}, Snow_Cover::Array{Float64,2}, Total_Snow_Melt::Array{Float64,1}
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
module TestReorderByModule
  using Base.Test
  using Brim

  # Test 1 - simple case

  A = [1 0 0; 0 0 1; 0 1 0]
  S = [0 1 0; 1 0 0; 0 0 1; 0 1 0; 0 0 1; 0 0 1]
  C = [1 0 0; 0 1 0; 0 0 1; 0 1 0; 0 0 1; 0 0 1]

  M = no_empty_modules!(Modular(A, S))

  reorder_by_module!(M)
  @test sum(abs(M.S .- C)) == 0

  # Test 2 - more complex case

  A = [0 0 0 1 1 1; 1 1 1 0 0 0; 0 0 0 1 1 1; 1 1 1 0 0 0; 0 0 0 1 1 1; 1 1 1 0 0 0]
  M = recursive_brim!(partition_random(A))
  reorder_by_module!(M)

  target = [1 0; 1 0; 1 0; 0 1; 0 1; 0 1; 1 0; 1 0; 1 0; 0 1; 0 1; 0 1]

  @test hash(M.S) == hash(target)

end

ext: jl


text:
parameters = [
    ("Temperatur","temperature"),
    ("Luftfeuchtigkeit","humidity"),
    ("Windgeschwindigkeit","wind_speed"),
    ("Windrichtung","wind_direction"),
    ("Regenzähler","rain")]

function server(;  basedir = get(ENV,"WEATHERSTATION_DIR","/var/lib/WeatherStation.jl/"),
                port = parse(Int,get(ENV,"WEATHERSTATION_PORT","8081")))

    app = dash("Wetterstation",
               external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"],
#               meta_tags=[
#                   Dict("name" => "viewport",
#                        "content" => "width=device-width, initial-scale=1")
#               ],
               )

    function make_graph_figure(parameter,taverage,timerange = (typemin(DateTime),typemax(DateTime)))
        df = loadavg(basedir,parameter,taverage; timerange = timerange)
        label = Dict([p[2] => p[1] for p in parameters]...)[parameter]

        (
            data = [(
                x = df[:,:time],
                y = df[:,Symbol(parameter)]
            )],
            layout = (title = "$(label) $(taverage/(1000*60)) min",)
        )
    end

    app.layout = html_div() do

        html_h2("Wetterstation"),
        dcc_dropdown(
            id="parameter",
            options=[Dict("label" => p[1], "value" => p[2]) for p in parameters],
            value="temperature"
        ),
        dcc_dropdown(
            id="taverage",
            options=[
                Dict("label" => "1 min", "value" => 1000 * 60 * 1),
                Dict("label" => "10 min", "value" => 1000 * 60 * 10),
                Dict("label" => "1 hour", "value" => 1000 * 60 * 60),
                Dict("label" => "24 hour", "value" => 1000 * 60 * 60 * 24),
            ],
            value=1000 * 60 * 10
        ),
        dcc_datepickerrange(
            id = "date-picker-range",
            min_date_allowed = DateTime(2020, 4, 1),
            #max_date_allowed = Dates.now() + Dates.Day(1),
            start_date = Dates.now() - Dates.Day(7),
            end_date = Dates.now() + Dates.Day(1),
        ),
        dcc_graph(
            id = "graph",
            #figure = make_graph_figure("temperature",1000 * 60 * 10),
            style = Dict("padding-top" => "20px")
        )
    end

    callback!(app, callid"parameter.value,taverage.value,date-picker-range.start_date,date-picker-range.end_date => graph.figure") do value, taverage, t0, t1
        timerange = (DateTime(t0),DateTime(t1))
        @show timerange
        return make_graph_figure(value,taverage,timerange)
    end

    @info("open http://localhost:$port/")
    run_server(app, "0.0.0.0", port)
end


text:
doc"""
    Pareto(α,θ)

The *Pareto distribution* with shape `α` and scale `θ` has probability density function

$f(x; \alpha, \theta) = \frac{\alpha \theta^\alpha}{x^{\alpha + 1}}, \quad x \ge \theta$

```julia
Pareto()            # Pareto distribution with unit shape and unit scale, i.e. Pareto(1.0, 1.0)
Pareto(a)           # Pareto distribution with shape a and unit scale, i.e. Pareto(a, 1.0)
Pareto(a, b)        # Pareto distribution with shape a and scale b

params(d)        # Get the parameters, i.e. (a, b)
shape(d)         # Get the shape parameter, i.e. a
scale(d)         # Get the scale parameter, i.e. b
```

External links
 * [Pareto distribution on Wikipedia](http://en.wikipedia.org/wiki/Pareto_distribution)

"""
immutable Pareto <: ContinuousUnivariateDistribution
    α::Float64
    θ::Float64

    function Pareto(α::Real, θ::Real)
        @check_args(Pareto, α > zero(α) && θ > zero(θ))
        new(α, θ)
    end
    Pareto(α::Real) = Pareto(α, 1.0)
    Pareto() = new(1.0, 1.0)
end

@distr_support Pareto d.θ Inf

#### Parameters

shape(d::Pareto) = d.α
scale(d::Pareto) = d.θ

params(d::Pareto) = (d.α, d.θ)

#### Statistics

mean(d::Pareto) = ((α, θ) = params(d); α > 1.0 ? α * θ / (α - 1.0) : Inf)
median(d::Pareto) = ((α, θ) = params(d); θ * 2.0 ^ (1.0 / α))
mode(d::Pareto) = d.θ

function var(d::Pareto)
    (α, θ) = params(d)
    α > 2.0 ? (θ^2 * α) / ((α - 1.0)^2 * (α - 2.0)) : Inf
end

function skewness(d::Pareto)
    α = shape(d)
    α > 3.0 ? ((2.0 * (1.0 + α)) / (α - 3.0)) * sqrt((α - 2.0) / α) : NaN
end

function kurtosis(d::Pareto)
    α = shape(d)
    α > 4.0 ? (6.0 * (α^3 + α^2 - 6.0 * α - 2.0)) / (α * (α - 3.0) * (α - 4.0)) : NaN
end

entropy(d::Pareto) = ((α, θ) = params(d); log(θ / α) + 1.0 / α + 1.0)

#### Evaluation

function pdf(d::Pareto, x::Float64)
    (α, θ) = params(d)
    x >= θ ? α * (θ / x)^α * (1.0 / x) : 0.0
end

function logpdf(d::Pareto, x::Float64)
    (α, θ) = params(d)
    x >= θ ? log(α) + α * log(θ) - (α + 1.0) * log(x) : -Inf
end

function ccdf(d::Pareto, x::Float64)
    (α, θ) = params(d)
    x >= θ ? (θ / x)^α : 1.0
end

cdf(d::Pareto, x::Float64) = 1.0 - ccdf(d, x)

function logccdf(d::Pareto, x::Float64)
    (α, θ) = params(d)
    x >= θ ? α * log(θ / x) : 0.0
end

logcdf(d::Pareto, x::Float64) = log1p(-ccdf(d, x))

cquantile(d::Pareto, p::Float64) = d.θ / p^(1.0 / d.α)
quantile(d::Pareto, p::Float64) = cquantile(d, 1.0 - p)

#### Sampling

rand(d::Pareto) = d.θ * exp(randexp() / d.α)

## Fitting

function fit_mle{T <: Real}(::Type{Pareto}, x::AbstractArray{T})
    # Based on
    # https://en.wikipedia.org/wiki/Pareto_distribution#Parameter_estimation

    θ = minimum(x)

    n = length(x)
    lθ = log(θ)
    temp1 = zero(T)
    for i=1:n
        temp1 += log(x[i]) - lθ
    end
    α = n/temp1

    return Pareto(α, θ)
end


text:
prob(p::LogProb) = p
prob(s::Product) = first(s)

#####################################################
### Classes of Traversal, Completion, and ItemKey ###
#####################################################

struct Traversal{S}
    edgeid  :: Int
    consid  :: Int
    score   :: S
    inloop  :: Bool
end
make_accesses(Traversal)

Traversal(edgeid, consid, score) = Traversal(edgeid, consid, score, false)
prob(trav::Traversal) = prob(score(trav))

abstract type Completion{S} end
prob(comp::Completion) = prob(score(comp))

struct EdgeCompletion{CR,S} <: Completion{S}
    edgeid :: Int
    rule   :: CR
    score  :: S
    inloop :: Bool
end
make_accesses(EdgeCompletion)

EdgeCompletion(edgeid, rule, score) = EdgeCompletion(edgeid, rule, score, false)
EdgeCompletion{CR,S}(edgeid, rule, score) where {CR,S} = EdgeCompletion{CR,S}(edgeid, rule, score, false)

struct TerminalCompletion{T,TR,S} <: Completion{S}
    terminal :: T
    rule     :: TR
    score    :: S
end
make_accesses(TerminalCompletion)

# is not a bitstype
@auto_hash_equals struct EdgeKey{St}
    start :: Int
    tail  :: Int
    state :: St
end
make_accesses(EdgeKey)

# is a bitstype
struct ConsKey{C}
    start :: Int
    tail  :: Int
    cat   :: C
end
make_accesses(ConsKey)

function ==(k1::ConsKey{T}, k2::ConsKey{T}) where {T <: AbstractString}
    k1.start == k2.start && k1.tail == k2.tail && k1.cat == k2.cat
end

function Base.hash(k::ConsKey{T}) where {T <: AbstractString}
    hash(k.start, hash(k.tail, hash(k.cat)))
end

##################
### Item class ###
##################

abstract type Item end
Item(key::EdgeKey, trav, id, grammar) = Edge(key, trav, id, grammar)
Item(key::ConsKey, comp, id, grammar) = Constituent(key, comp, id, grammar)
start(item::Item) = item.start
tail(item::Item) = item.tail
length(item::Item) = item.tail - item.start
id(item::Item) = item.id
lastpopscore(item::Item) = item.lastpopscore
insidepopnumber(item::Item) = item.insidepopnumber
outsidepopnumber(item::Item) = item.outsidepopnumber

##################
### Edge class ###
##################

mutable struct Edge{St,S} <: Item
    start              :: Int
    tail               :: Int
    state              :: St
    score              :: S
    isfinished         :: Bool
    traversals         :: Vector{Traversal{S}}
    id                 :: Int # edges have positive ids
    lastpopscore       :: S
    lastoutsidepopscore:: S
    insidepopnumber    :: Int
    outsidepopnumber   :: Int
    outside_finished   :: Bool
    trav_outs_summands :: Dict{Tuple{Int, Int}, LogProb}
    comp_outs_summands :: Dict{Int, LogProb}
end
Edge(start, tail, state, traversals::Vector, id, score) =
    Edge(start, tail, state, score, false, traversals, id, score, score, 0, 0, false, Dict{Tuple{Int, Int}, LogProb}(), Dict{Int, LogProb}())
Edge(start, tail, state, traversal::Traversal{S}, id) where S =
    Edge(start, tail, state, [traversal], id, zero(S))
Edge(key::EdgeKey, traversal::Traversal, id, grammar) =
    Edge(start(key), tail(key), state(key), traversal, id)

score(edge::Edge) = if edge.isfinished
    edge.score
else
    sum(score(trav) for trav in traversals(edge))
end
prob(edge::Edge) = prob(score(edge))
outside(edge::Edge) = sum(values(edge.trav_outs_summands)) + sum(values(edge.comp_outs_summands))
state(edge::Edge) = edge.state
cat_or_state(edge::Edge) = edge.state
key(edge::Edge) = EdgeKey(edge.start, edge.tail, edge.state)
traversals(edge::Edge) = edge.traversals
show(io::IO, edge::Edge) = print(io, state(edge), "(", start(edge), ",", tail(edge), ")")

function add!(edge::Edge, traversal)
    for (i, trav) in enumerate(traversals(edge))
        if edgeid(trav) == edgeid(traversal) && consid(trav) == consid(traversal)
            traversals(edge)[i] = @set traversal.inloop = true
            return edge
        end
    end
    push!(traversals(edge), traversal)
    edge
end

#########################
### Constituent class ###
#########################

mutable struct Constituent{C,T,CR,TR,S} <: Item
    start               :: Int
    tail                :: Int
    cat                 :: C
    score               :: S
    isfinished          :: Bool
    completions         :: Vector{EdgeCompletion{CR,S}}
    terminal_completion :: Nullable{TerminalCompletion{T,TR,S}}
    id                  :: Int # constituents have negative ids
    lastpopscore        :: S
    lastoutsidepopscore :: S
    insidepopnumber     :: Int
    outsidepopnumber    :: Int
    outside_finished    :: Bool
    trav_outs_summands  :: Dict{Tuple{Int, Int}, LogProb}
end
Constituent(st,ta,ca,co,te,id,sc) =
    Constituent(st,ta,ca,sc,false,co,te,id,sc,sc,0,0,false,Dict{Tuple{Int, Int}, LogProb}())
function Constituent(start, tail, cat, comp::EdgeCompletion, id, grammar)
    T, TR, S = terminal_type(grammar), terminal_rule_type(grammar), score_type(grammar)
    Constituent(start, tail, cat, [comp], Nullable{TerminalCompletion{T,TR,S}}(), id, zero(S))
end
function Constituent(start, tail, cat, comp::TerminalCompletion, id, grammar)
    CR, S = category_rule_type(grammar), score_type(grammar)
    Constituent(start, tail, cat, Vector{EdgeCompletion{CR,S}}(), Nullable(comp), id, zero(S))
end
Constituent(key::ConsKey, comp, id, grammar) =
    Constituent(start(key), tail(key), cat(key), comp, id, grammar)

score(cons::Constituent) = if cons.isfinished
    cons.score
else
    if has_terminal(cons)
        if isempty(completions(cons))
            score(get(cons.terminal_completion)) # get(cons.terminal_completion))
        else
            sum(score(comp) for comp in completions(cons)) +
                score(get(cons.terminal_completion)) # get(cons.terminal_completion))
        end
    else
        sum(score(comp) for comp in completions(cons))
    end
end

prob(cons::Constituent) = prob(score(cons))
outside(cons::Constituent) = sum(values(cons.trav_outs_summands))
cat(cons::Constituent) = cons.cat
cat_or_state(cons::Constituent) = cons.cat
key(cons::Constituent) = ConsKey(cons.start, cons.tail, cons.cat)
completions(cons::Constituent) = cons.completions
has_terminal(cons::Constituent) = !isnull(cons.terminal_completion)
terminal_completion(cons::Constituent) = get(cons.terminal_completion)
terminal(cons::Constituent) = terminal(terminal_completion(cons))

function add!(cons::Constituent, completion::EdgeCompletion)
    for (i, comp) in enumerate(completions(cons))
        if edgeid(comp) == edgeid(completion)
            completions(cons)[i] = @set completion.inloop = true
            return cons
        end
    end
    push!(completions(cons), completion)
    nothing
end

function show(io::IO, c::Constituent)
    print(io, cat(c), "(", start(c), ",", tail(c), ")")
    if has_terminal(c)
        print(io, "[", terminal(c) ,"]")
    end
end

###########################
### ParserLogbook class ###
###########################

mutable struct ParserLogbook{C,T,CR,TR,St,S}
    edges :: Vector{Edge{St,S}}
    conss :: Vector{Constituent{C,T,CR,TR,S}}
    edgeids :: Dict{EdgeKey{St}, Int}
    consids :: Dict{ConsKey{C}, Int}
end
ParserLogbook(C, T, CR, TR, St, S) =
    ParserLogbook(Vector{Edge{St,S}}(),
                  Vector{Constituent{C,T,CR,TR,S}}(),
                  Dict{EdgeKey{St}, Int}(),
                  Dict{ConsKey{C}, Int}())

function discover!(logbook, edge::Edge)
    push!(logbook.edges, edge)
    id = length(logbook.edges)
    logbook.edgeids[key(edge)] = id
    edge.id = id
    nothing
end
function discover!(logbook, cons::Constituent)
    push!(logbook.conss, cons)
    id = -length(logbook.conss)
    logbook.consids[key(cons)] = id
    cons.id = id
    nothing
end
isdiscovered(logbook, key::EdgeKey) = haskey(logbook.edgeids, key)
isdiscovered(logbook, key::ConsKey) = haskey(logbook.consids, key)
get_edge(logbook, id::Int) = logbook.edges[id]
get_cons(logbook, id::Int) = logbook.conss[-id]
get_edge(logbook, key::EdgeKey) = get_edge(logbook, logbook.edgeids[key])
get_cons(logbook, key::ConsKey) = get_cons(logbook, logbook.consids[key])
get_item(logbook, key::EdgeKey) = get_edge(logbook, logbook.edgeids[key])
get_item(logbook, key::ConsKey) = get_cons(logbook, logbook.consids[key])

###################
### Chart class ###
###################

# maybe be used in a future implementation, maybe not
"like normal merge, but value vectors with the same key will be concatenated"
function Base.merge(dicts::Dict{K,Vector{V}}...) where {K,V}
    result = Dict{K,Vector{V}}()
    for d in dicts
        for (key,valvec) in d # key and value vector
            if haskey(result, key)
                append!(result[key], valvec)
            else
                result[key] = valvec
            end
        end
    end
    result
end

mutable struct ChartCell{C,St}
    edgeids :: Dict{St, Vector{Int}}
    consids :: Dict{C, Vector{Int}}
end
ChartCell(C, St) =
    ChartCell(Dict{St, Vector{Int}}(), Dict{C, Vector{Int}}())

mutable struct Chart{C,St}
    cells :: Vector{ChartCell{C,St}}
end
edgeids(chart, edge::Edge) = chart.cells[tail(edge)].edgeids
consids(chart, edge::Edge) = chart.cells[tail(edge)].consids
edgeids(chart, cons::Constituent) = chart.cells[start(cons)].edgeids
consids(chart, cons::Constituent) = chart.cells[start(cons)].consids
# dict(chart, edge::Edge) = chart.cells[tail(edge)].edgeids
# dict(chart, cons::Constituent) = chart.cells[start(cons)].consids
function insert!(chart::Chart, edge::Edge)
    dict = edgeids(chart, edge)
    if haskey(dict, state(edge))
        push!(dict[state(edge)], id(edge))
    else
        dict[state(edge)] = [id(edge)]
    end
    nothing
end
function insert!(chart::Chart, cons::Constituent)
    dict = consids(chart, cons)
    if haskey(dict, cat(cons))
        push!(dict[cat(cons)], id(cons))
    else
        dict[cat(cons)] = [id(cons)]
    end
    nothing
end

####################
### Agenda class ###
####################

# fast viterbi parsing does not work yet

mutable struct Agenda{T}
    pqueue :: PriorityQueue{Int, T}
end

function Agenda(parsing_method)
    @assert parsing_method in (:full, :viterbi)
    if parsing_method == :full
        Agenda(PriorityQueue{Int, Int}())
    else
        Agenda(PriorityQueue{Int, Float64}())
    end
end

function enqueue!(agenda::Agenda, item, just_used)
    agenda.pqueue[id(item)] = priority(agenda, item, just_used) ::Int
    nothing
end

dequeue!(agenda::Agenda) = dequeue!(agenda.pqueue)
isempty(agenda::Agenda) = isempty(agenda.pqueue)

priority(agenda::Agenda{Int}, edge::Edge, just_completed::Bool) =
    4 * length(edge) - 2*!just_completed - 1
priority(agenda::Agenda{Int}, cons::Constituent, just_introduced_edge::Bool) =
    4 * length(cons) - 2*!just_introduced_edge

priority(agenda::Agenda{Float64}, item, just_used) =
    - prob(item).value # negative logprob, agenda favours high probs but is min orientated

##########################
### ParseForest class ###
##########################

mutable struct ParseForest{C,T,CR,TR,St,S}
    heads :: Vector{Constituent{C,T,CR,TR,S}}
    logbook :: ParserLogbook{C,T,CR,TR,St,S}
    terminals :: Vector{T}
end

function ParseForest(chart::Chart, logbook::ParserLogbook{C,T,CR,TR,St,S}, terminals, grammar) where {C,T,CR,TR,St,S}
    heads = Constituent{C,T,CR,TR,S}[]
    cell = chart.cells[1]
    for cat in startsymbols(grammar)
        if haskey(cell.consids, cat)
            for consid in cell.consids[cat]
                cons = get_cons(logbook, consid)
                if tail(cons) == length(terminals)+1
                    push!(heads, cons)
                end
            end
        end
    end
    ParseForest(heads, logbook, terminals)
end

heads(f::ParseForest) = f.heads
is_complete(f::ParseForest) = !isempty(f.heads)
score(f::ParseForest) = sum(map(score, f.heads))

### OutsideAgenda ###

mutable struct OutsideAgenda
    pqueue :: PriorityQueue{Int, Int}
end

function OutsideAgenda()
    OutsideAgenda(PriorityQueue{Int, Int}())
end

function enqueue!(agenda::OutsideAgenda, item, just_used)
    agenda.pqueue[id(item)] = priority(agenda, item, just_used)
    agenda
end

dequeue!(agenda::OutsideAgenda) = dequeue!(agenda.pqueue)
isempty(agenda::OutsideAgenda) = isempty(agenda.pqueue)

priority(agenda::OutsideAgenda, edge::Edge, just_used::Bool) =
    - 4 * length(edge) - 2*!just_used -1
priority(agenda::OutsideAgenda, cons::Constituent, just_used::Bool) =
    - 4 * length(cons) - 2*!just_used

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

inside = prob

@inline function acc_or_init!(d::Dict, k, v)
    if haskey(d, k)
        d[k] += v
    else
        d[k] = v
    end
end

function calculate_outside_probs_and_expected_rule_usages!(f::ParseForest{C,T,CR,TR,St,S}, grammar) where {C,T,CR,TR,St,S}
    @assert is_complete(f)
    outside_agenda = OutsideAgenda()
    # mapping categories to rules to (startindex,endindex) to probabilities
    rule_usage_probs = Dict{C, Dict{typejoin(CR,TR), Dict{Tuple{Int,Int}, LogProb}}}()
    for head in heads(f)
        # traversal outside score summands
        head.trav_outs_summands[0,0] = one(LogProb)
        enqueue!(outside_agenda, head, false)
    end
    while !isempty(outside_agenda)
        id = dequeue!(outside_agenda)
        if id > 0
            parent_edge = get_edge(f.logbook, id)
            parent_edge.outsidepopnumber += 1
            # @show parent_edge.id
            # println()
            o = outside(parent_edge)
            if !(outside(parent_edge) ≈ parent_edge.lastoutsidepopscore) &&
                outsidepopnumber(parent_edge) < MAX_POP_NUMBER
                for trav in traversals(parent_edge)
                    if trav.inloop
                        cons = get_cons(f.logbook, consid(trav))
                        if edgeid(trav) == 0
                            acc_or_init!(cons.trav_outs_summands, (id,0), o)
                            # cons.trav_outs_summands[id, 0] = o
                            enqueue!(outside_agenda, cons, false)
                        end
                    end
                end
                parent_edge.lastoutsidepopscore = o
                enqueue!(outside_agenda, parent_edge, true)
            elseif !(parent_edge.outside_finished)
                parent_edge.outside_finished = true
                for trav in traversals(parent_edge)
                    if true # !(trav.inloop)
                        # println("edge")
                        cons = get_cons(f.logbook, consid(trav))
                        if edgeid(trav) == 0
                            # acc_or_init!(cons.trav_outs_summands, (id,0), o)
                            cons.trav_outs_summands[id, 0] = o
                            enqueue!(outside_agenda, cons, false)
                        else
                            edge = get_edge(f.logbook, edgeid(trav))
                            edge.trav_outs_summands[id,cons.id] = o * inside(cons)
                            # acc_or_init!(edge.trav_outs_summands, (id, cons.id), o * inside(cons))
                            enqueue!(outside_agenda, edge, false)
                            cons.trav_outs_summands[id,edge.id] = o * inside(edge)
                            # acc_or_init!(cons.trav_outs_summands, (id, edge.id), o * inside(edge))
                            enqueue!(outside_agenda, cons, false)
                        end
                    end
                end
            end
        else
            parent_cons = get_cons(f.logbook, id)
            parent_cons.outsidepopnumber += 1
            # @show parent_cons
            # println(map(c->c.edgeid, parent_cons.completions))
            # println(exp(score(parent_cons)))
            # println(parent_cons.trav_outs_summands)
            # println(outside(parent_cons) ≈ parent_cons.lastoutsidepopscore)
            # println()
            o = outside(parent_cons)
            if !(outside(parent_cons) ≈ parent_cons.lastoutsidepopscore) &&
                outsidepopnumber(parent_cons) < MAX_POP_NUMBER
                for comp in completions(parent_cons)
                    if comp.inloop
                        edge = get_edge(f.logbook, edgeid(comp))
                        p = prob(grammar, cat(parent_cons), rule(comp))
                        acc_or_init!(edge.comp_outs_summands, id, o * p)
                        # edge.comp_outs_summands[id] = o * p
                        enqueue!(outside_agenda, edge, false)
                    end
                end
                parent_cons.lastoutsidepopscore = o
                enqueue!(outside_agenda, parent_cons, true)
            else#if !(parent_cons.outside_finished)
                parent_cons.outside_finished = true
                for comp in completions(parent_cons)
                    if true # !(comp.inloop)
                        # println("cons")
                        edge = get_edge(f.logbook, edgeid(comp))
                        p = prob(grammar, cat(parent_cons), rule(comp))
                        # acc_or_init!(edge.comp_outs_summands, id, o * p)
                        edge.comp_outs_summands[id] = o * p
                        enqueue!(outside_agenda, edge, false)
                        if haskey(rule_usage_probs, cat(parent_cons))
                            if haskey(rule_usage_probs[cat(parent_cons)], rule(comp))
                                # if haskey(rule_usage_probs[cat(parent_cons)][rule(comp)], (start(parent_cons),tail(parent_cons)))
                                #     rule_usage_probs[cat(parent_cons)][rule(comp)][(start(parent_cons),tail(parent_cons))] += o * p * inside(edge)
                                # else
                                #     rule_usage_probs[cat(parent_cons)][rule(comp)][(start(parent_cons),tail(parent_cons))] = o * p * inside(edge)
                                # end
                                rule_usage_probs[cat(parent_cons)][rule(comp)][(start(parent_cons),tail(parent_cons))] = o * p * inside(edge)
                            else
                                rule_usage_probs[cat(parent_cons)][rule(comp)] = Dict((start(parent_cons),tail(parent_cons)) => o * p * inside(edge))
                            end
                        else
                            rule_usage_probs[cat(parent_cons)] = Dict(rule(comp)=>Dict((start(parent_cons),tail(parent_cons)) => o * p * inside(edge)))
                        end
                    end
                end
                if has_terminal(parent_cons)
                    comp = terminal_completion(parent_cons)
                    p = prob(grammar, cat(parent_cons), rule(comp))
                    if haskey(rule_usage_probs, cat(parent_cons))
                        if haskey(rule_usage_probs[cat(parent_cons)], rule(comp))
                            # if haskey(rule_usage_probs[cat(parent_cons)][rule(comp)], (start(parent_cons),tail(parent_cons)))
                            #     rule_usage_probs[cat(parent_cons)][rule(comp)][(start(parent_cons),tail(parent_cons))] += o * p
                            # else
                            #     rule_usage_probs[cat(parent_cons)][rule(comp)][(start(parent_cons),tail(parent_cons))] = o * p
                            # end
                            rule_usage_probs[cat(parent_cons)][rule(comp)][(start(parent_cons),tail(parent_cons))] = o * p
                        else
                            rule_usage_probs[cat(parent_cons)][rule(comp)] = Dict((start(parent_cons),tail(parent_cons)) => o * p)
                        end
                    else
                        rule_usage_probs[cat(parent_cons)] = Dict(rule(comp)=>Dict((start(parent_cons),tail(parent_cons)) => o * p))
                    end
                end
            end
        end
    end

    sequence_prop = sum(inside(head) for head in f.heads)

    expected_rule_usages = Dict(
        cat => Dict(
            rule => exp(sum(values(rule_usage_probs[cat][rule])).value - sequence_prop.value)
            for rule in keys(rule_usage_probs[cat]))
        for cat in keys(rule_usage_probs)
    )
end

function calculate_outside_probs_and_expected_rule_usages!(sent :: Vector{T}, grammar :: Grammar{C,T,CR,TR,S,Cond}) where {C,T,CR,TR,S,Cond}
    forest = run_chartparser(sent, grammar)
    d = calculate_outside_probs_and_expected_rule_usages!(forest, grammar)
    (forest, d)
end

function best_tree(f::ParseForest)
    @assert is_complete(f)
    best_tree(f, f.heads[findmax(map(prob, f.heads))[2]])
end

function best_tree(f::ParseForest{C,T,CR,TR,St,S}, cons::Constituent) where {C,T,CR,TR,St,S}
    if has_terminal(cons)
        TreeNode((cons,rule(terminal_completion(cons))),
                 Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}})
    else
        comp = completions(cons)[findmax(map(prob, completions(cons)))[2]]
        node = TreeNode((cons, rule(comp)),
                        Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}})
        trees = best_trees(f, get_edge(f.logbook, edgeid(comp))) # children trees
        for t in trees
            insert_child!(node, t)
        end
        node
    end
end

function best_trees(f::ParseForest, edge::Edge)
    best_traversal = traversals(edge)[findmax(map(prob, traversals(edge)))[2]]
    if edgeid(best_traversal) == 0
        (best_tree(f, get_cons(f.logbook, consid(best_traversal))),)
    else
        (best_trees(f, get_edge(f.logbook, edgeid(best_traversal)))...,
         best_tree(f,  get_cons(f.logbook, consid(best_traversal))))
    end
end

function sample_tree(f::ParseForest)
    @assert is_complete(f)
    head = categorical_sample(heads(f), map(prob, heads(f)))
    sample_tree(f, head)
end

function sample_tree(f::ParseForest{C,T,CR,TR,St,S}, cons::Constituent) where{C,T,CR,TR,St,S}
    if has_terminal(cons)
        TreeNode((cons,rule(terminal_completion(cons))),
                 Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}})
    else
        comp = categorical_sample(completions(cons),
                                           map(prob, completions(cons)))
        node = TreeNode((cons, rule(comp)),
                        Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}})
        trees = sample_trees(f, get_edge(f.logbook, edgeid(comp))) # children trees
        for t in trees
            insert_child!(node, t)
        end
        node
    end
end

function sample_trees(f::ParseForest, edge::Edge)
    traversal = categorical_sample(traversals(edge),
                                            map(prob, traversals(edge)))
    if edgeid(traversal) == 0 # then this edge contains the start state
        (sample_tree(f, get_cons(f.logbook, consid(traversal))),)
    else
        (sample_trees(f, get_edge(f.logbook, edgeid(traversal)))...,
         sample_tree(f,  get_cons(f.logbook, consid(traversal))))
    end
end

# # all_trees no yet implemented
# function all_trees{C,T,CR,TR,St,S}(f::ParseForest{C,T,CR,TR,St,S})
#     if is_complete(f)
#         vcat([all_trees(f, cons) for cons in heads(f)]...)
#     else
#         Vector{TreeNode{Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}}}}
#     end
# end
#
# function all_trees{C,T,CR,TR,St,S}(
#     f    :: ParseForest{C,T,CR,TR,St,S},
#     cons :: Constituent
#     )
#     if has_terminal(cons)
#         [TreeNode((cons,rule(terminal_completion(cons))),
#                   Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}})]
#     else
#         for comp in completions(cons)
#             node = TreeNode((cons, rule(comp)),
#                             Tuple{Constituent{C,T,CR,TR,S}, Union{CR,TR}})
#         end
#     end
# end
#
# function all_trees(f::ParseForest, edge::Edge)
#
# end

#######################
### ParseTree class ###
#######################

ParseTree{Cons<:Constituent,R} = TreeNode{Tuple{Cons,R}}

function show(io::IO, tree::ParseTree)
    print(io, "[", tree.data[1])
    for child in tree.children
        print(io, child)
    end
    print(io, "]")
end

cats_and_rules(t::ParseTree) = ((node.data[1].cat, node.data[2]) for node in t)


text:
using Documenter, SignalingDimension

DocMeta.setdocmeta!(SignalingDimension, :DocTestSetup, :(using SignalingDimension); recursive=true)

using BellScenario
using QBase

include("../script/utilities.jl")
using .ScriptUtilities

makedocs(;
    modules=[SignalingDimension, ScriptUtilities],
    authors="Brian Doolittle <brian.d.doolittle@gmail.com>",
    repo="https://github.com/ChitambarLab/SignalingDimension.jl/blob/{commit}{path}#L{line}",
    sitename="SignalingDimension.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ChitambarLab.github.io/SignalingDimension.jl",
        assets=String["assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
        "Background" => [
            "Signaling Correlations" => "background/signaling_correlations.md",
            "Signaling Polytopes" => "background/signaling_polytopes.md",
            "Signaling Dimension" => "background/signaling_dimension.md",
        ],
        "Bell Inequalities" => "bell_inequalities.md",
        "Certifying Signaling Dimension" => "certifying_signaling_dimension.md",
        "Script Utilities" => "script_utilities.md",
    ],
)

deploydocs(;
    repo="github.com/ChitambarLab/SignalingDimension.jl",
    devbranch = "main",
)


text:
module LuaCall

using Lua_jll

export luacall, @lua_str

const LUA_STATE = Ref{Ptr{Cvoid}}(C_NULL)

const LuaInt = let t = ccall((:jl_lua_int_type, liblua), Cint, ())
    t == 1 ? Cint :
    t == 2 ? Clong :
    t == 3 ? Clonglong :
    error("unknown Lua integer type")
end

const LuaFloat = let t = ccall((:jl_lua_float_type, liblua), Cint, ())
    t == 1 ? Cfloat :
    t == 2 ? Cdouble :
    t == 3 ? error("Julia does not natively support long doubles") :
    error("unknown Lua float type")
end

function __init__()
    state = ccall((:luaL_newstate, liblua), Ptr{Cvoid}, ())
    state == C_NULL && error("Failed to initialize Lua")
    ccall((:luaL_openlibs, liblua), Cvoid, (Ptr{Cvoid},), state)
    global LUA_STATE[] = state
    atexit(()->ccall((:lua_close, liblua), Cvoid, (Ptr{Cvoid},), LUA_STATE[]))
    nothing
end

struct Table
    name::Symbol
end

struct LuaError <: Exception
    msg::String

    LuaError(msg) = new(msg)
    function LuaError()
        msg = unsafe_getstack(String, -1)
        ccall((:lua_settop, liblua), Cvoid, (Ptr{Cvoid}, Cint), LUA_STATE[], -2)
        new(msg)
    end
end

Base.showerror(io::IO, ex::LuaError) = print(io, "LuaError: ", ex.msg)

#define LUA_TNONE		(-1)
#define LUA_TNIL		0
#define LUA_TBOOLEAN		1
#define LUA_TLIGHTUSERDATA	2
#define LUA_TNUMBER		3
#define LUA_TSTRING		4
#define LUA_TTABLE		5
#define LUA_TFUNCTION		6
#define LUA_TUSERDATA		7
#define LUA_TTHREAD		8

function checkstack(size::Int=1)
    ok = ccall((:lua_checkstack, liblua), Cint, (Ptr{Cvoid}, Cint), LUA_STATE[], size)
    ok == 1 || throw(LuaError("insufficient space on the Lua stack for $size items"))
    nothing
end

function stacktype(pos::Int=-1)
    t = ccall((:lua_type, liblua), Cint, (Ptr{Cvoid}, Cint), LUA_STATE[], pos)
    if t <= 0
        Nothing  # TODO: Differentiate nil and none?
    elseif t == 1
        Bool
    elseif t == 3
        int = ccall((:lua_isinteger, liblua), Cint, (Ptr{Cvoid}, Cint), LUA_STATE[], pos)
        int == 0 ? LuaFloat : LuaInt
    elseif t == 4
        String
    else
        name = ccall((:lua_typename, liblua), Cstring, (Ptr{Cvoid}, Cint), LUA_STATE[], t)
        error("Lua type ", unsafe_string(name), " is not currently supported")
    end
end

pushstack(value::Integer) =
    ccall((:lua_pushinteger, liblua), Cvoid, (Ptr{Cvoid}, LuaInt), LUA_STATE[], value)

pushstack(value::Bool) =
    ccall((:lua_pushboolean, liblua), Cvoid, (Ptr{Cvoid}, Cint), LUA_STATE[], value)

pushstack(value::Real) =
    ccall((:lua_pushnumber, liblua), Cvoid, (Ptr{Cvoid}, LuaFloat), LUA_STATE[], value)

pushstack(::Nothing) =
    ccall((:lua_pushnil, liblua), Cvoid, (Ptr{Cvoid},), LUA_STATE[])

pushstack(value::AbstractString) =
    ccall((:lua_pushlstring, liblua), Cvoid,
          (Ptr{Cvoid}, Ptr{UInt8}, Csize_t),
          LUA_STATE[], value, length(value))

function getstack(::Type{T}, pos::Int=-1) where T
    checkstack(-pos)
    unsafe_getstack(T, pos)
end

getstack() = getstack(stacktype(-1), -1)

function unsafe_getstack(::Type{T}, pos::Int) where T<:Integer
    isnum = Ref{Cint}(-1)
    x = ccall((:lua_tointegerx, liblua), LuaInt,
              (Ptr{Cvoid}, Cint, Ref{Cint}),
              LUA_STATE[], pos, isnum)
    @assert isnum[] == 1
    x
end

unsafe_getstack(::Type{Bool}, pos::Int) =
    ccall((:lua_toboolean, liblua), Cint, (Ptr{Cvoid}, Cint), LUA_STATE[], pos) != 0

function unsafe_getstack(::Type{T}, pos::Int) where T<:Real
    isnum = Ref{Cint}(-1)
    x = ccall((:lua_tonumberx, liblua), LuaFloat,
              (Ptr{Cvoid}, Cint, Ref{Cint}),
              LUA_STATE[], pos, isnum)
    @assert isnum[] == 1
    x
end

unsafe_getstack(::Type{Nothing}, pos::Int) = nothing

function unsafe_getstack(::Type{String}, pos::Int)
    len = Ref{Csize_t}(0)
    ptr = ccall((:lua_tolstring, liblua), Ptr{UInt8},
                (Ptr{Cvoid}, Cint, Ref{Csize_t}),
                LUA_STATE[], pos, len)
    unsafe_string(ptr, len[])
end

function pcallk(f::Symbol, rt::NTuple{N,DataType}, args...) where N
    checkstack(length(args) + N + 1)
    #stackdump()
    ccall((:lua_getglobal, liblua), Cvoid, (Ptr{Cvoid}, Cstring), LUA_STATE[], f)
    #stackdump()
    for arg in args
        pushstack(arg)
    end
    #stackdump()
    rc = ccall((:lua_pcallk, liblua), Cint,
               (Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}),
               LUA_STATE[], length(args), N, 0, C_NULL, C_NULL)
    rc == 0 || throw(LuaError())
    #stackdump()
    rc
end

function luacall(f::Symbol, rt::DataType, args...)
    rc = pcallk(f, (rt,), args...)
    @assert rc == 0
    result = getstack(rt)
    ccall((:lua_settop, liblua), Cvoid, (Ptr{Cvoid}, Cint), LUA_STATE[], -2)
    #stackdump()
    result
end

function luacall(f::Symbol, rt::NTuple{N,DataType}, args...) where N
    rc = pcallk(f, rt, args...)
    @assert rc == 0
    results = Vector{Any}(undef, N)
    @inbounds for i = N:-1:1
        results[i] = getstack(rt[N-i+1], -i)
    end
    ccall((:lua_settop, liblua), Cvoid, (Ptr{Cvoid}, Cint), LUA_STATE[], -N - 1)
    #stackdump()
    results
end

macro lua_str(code::String)
    rc = ccall((:luaL_loadstring, liblua), Cint, (Ptr{Cvoid}, Cstring), LUA_STATE[], code)
    rc == 0 || throw(LuaError())
    # Do an initial "priming" call, which ensures that Lua tells C everything it needs to
    # know about what has been defined
    rc = ccall((:lua_pcallk, liblua), Cint,
               (Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}),
               LUA_STATE[], 0, 0, 0, C_NULL, C_NULL)
    rc == 0 || throw(LuaError())
end

function stackdump()
    top = ccall((:lua_gettop, liblua), Cint, (Ptr{Cvoid},), LUA_STATE[])
    println("Lua stack (", top, " items):")
    for i = 1:top
        t = ccall((:lua_type, liblua), Cint, (Ptr{Cvoid}, Cint), LUA_STATE[], i)
        n = unsafe_string(ccall((:lua_typename, liblua), Ptr{UInt8},
                                (Ptr{Cvoid}, Cint), LUA_STATE[], t))
        print("  [", i, "]: ", n, " ")
        if t == 4
            len = Ref{Csize_t}(0)
            ptr = ccall((:lua_tolstring, liblua), Ptr{UInt8},
                        (Ptr{Cvoid}, Cint, Ref{Csize_t}),
                        LUA_STATE[], i, len)
            s = unsafe_string(ptr, len[])
            println(repr(s))
        elseif t == 1
            b = ccall((:lua_toboolean, liblua), Cint, (Ptr{Cvoid}, Cint), LUA_STATE[], i)
            println(b != 0)
        elseif t == 3
            isnum = Ref{Cint}(-1)
            x = ccall((:lua_tonumberx, liblua), Cdouble,
                      (Ptr{Cvoid}, Cint, Ref{Cint}),
                      LUA_STATE[], i, isnum)
            @assert isnum[] == 1
            println(x)
        else
            println()
        end
    end
    println()
    nothing
end

end # module


text:
import Base.:*, Base.:>, Base.:<
using AbstractAlgebra: GFElem

export term,
    Field,
    binom,
    gt,
    gte,
    lt,
    |,
    /,
    div,
    *,
    lcm,
    minus,
    LT,
    S,
    mdiv,
    buchberger,
    is_groebner_basis,
    tuples_lt

Field = GF(32003)
/(x::AbstractAlgebra.GFElem{Int64}, y::AbstractAlgebra.GFElem{Int64}) = x * inv(y)

# A term is a pair of coefficient (Float64) and exponent
struct term{N}
    l :: GFElem{Int64}
    a :: NTuple{N, Int64}
end

binom(t :: term) = [t, term(Field(0), ntuple(x->0, length(t.a)))]
binom(t :: Vector{term{N}}) where N = length(t) == 1 ? binom(t[1]) : t

# Grevlex order
#
# a > b if total (sum?) degree of a > b or,
# if equal the last non-zero term of a - b is negative
gt(a::term, b::term) = gt(a.a, b.a)
function gt(a, b)
    if a == b                return false
    elseif (sum(a) > sum(b)) return true
    elseif (sum(a) < sum(b)) return false
    else   return filter(!=(0), a .- b)[end] < 0 end
end

gte(a::term, b::term) = a.a == b.a || gt(a, b)
lt(a::term, b::term) = ! gte(a, b)

>(a :: term, b :: term) = gt(a.a, b.a)
<(a :: term, b :: term) = not(gt(a.a, b.a))

# Term division
(|)(t :: term, r :: term) = all(t.a .<= r.a) # Ja, std <=, ikke grevlex ;)
div = (|)
(/)(t :: term, r :: term) = term(/(t.l, r.l), t.a .- r.a)

# Term lcm
lcm(t :: term, r :: term) = term(one(t.l), max.(t.a, r.a))

# Polynomial-term multiplication
(*)(t :: term, r :: term) = term(t.l * r.l, t.a .+ r.a)
(*)(t :: term, f ) = Ref(t) .* f

# Polynomial subtraction
minus(f, g) = let d = merge(+, Dict(x.a => x.l for x in f),
                              Dict(x.a => -1*x.l for x in g))
    [term(v, k) for (k, v) in d if v != 0]
end
add(f, g) = let d = merge(+, Dict(x.a => x.l for x in f),
                          Dict(x.a => x.l for x in g))
    [term(v, k) for (k, v) in d if v != 0]
end

# Leading term
#
# (lambda, a) at the largest a according to grevlex
LT(f) = reduce((a, b) -> if (gt(a.a, b.a)) a else b end, f)

# S(f, g) = let gamma = lcm(LT(f), LT(g)).a
#     ((term(1, gamma)/LT(f))*f) - ((term(1, gamma)/LT(g))*g)
# end

S(f, g) = let gamma = lcm(LT(f), LT(g))
    minus((gamma/LT(f))*f, (gamma/LT(g))*g)
end

# Multivariate division

function mdiv_count(h, F)
    # r = h
    # count = 0
    # sort!(F, lt = (x, y) -> ! gt(x[1].a, y[1].a))
    # while r != [] && any(LT.(F) .| Ref(LT(r)))
    #     i = findfirst(LT.(F) .| Ref(LT(r)))
    #     f = F[i]
    #     r = minus(r, (LT(r)/LT(f))*f)
    #     count = count + 1
    #     # println(r)
    # end
    # return (r, count)

    r = minus(h, h)

    count = 0
    while length(h) != 0
        if any(LT.(F) .| Ref(LT(h)))
            i = findfirst(LT.(F) .| Ref(LT(h)))
            f = F[i]
            h = minus(h, (LT(h) / LT(f)) * f)
            count = count + 1
        else
            r = add(r, [LT(h)])
            h = minus(h, [LT(h)])
        end
    end
    return (add(r, h), count)

end

mdiv(h, F) = mdiv_count(h, F)[1]

# reduce = mdiv

# Update function; TODO: refine this
# BEWARE: mutation!!!
# update!(P, G, r) = append!(P, [(f, r) for f in G])

function update!(P, G, r)
    m = length(G)
    ltG = LT.(G)
    ltr = LT(r)
    function can_drop(p)
        f, g = p
        gam = lcm(LT(f), LT(g))
        return (ltr | gam) && gam != lcm(LT(f), ltr) && gam != lcm(LT(g), ltr)
    end
    filter!(p -> !can_drop(p), P)

    # println("filtered")
    lcms = Dict{term, Vector{Vector{term}}}()
    for f in G
        if lcm(LT(f), ltr) in keys(lcms)
            push!(lcms[lcm(LT(f), ltr)], f)
        else
            lcms[lcm(LT(f), ltr)] = [f]
        end
    end

    # println("computed lcms")
    # println(length(keys(lcms)))
    min_lcms = []
    P_ = []
    for gam in sort(collect(keys(lcms)), lt = ((x, y) -> ! gt(x.a, y.a)))
    # for gam in keys(lcms)
        if all([! (gam | t) for t in min_lcms])
            push!(min_lcms, gam)
            if ! any([lcm(LT(t), ltr) == LT(t) * ltr for t in lcms[gam] ])
                push!(P_, (lcms[gam][1], r))
            end
        end
    end
    append!(P, P_)
end

select!(P) = pop!(P)

function buchberger(G)
    P = [(f, g) for f in G for g in G]
    t = 0
    while length(P) > 0
        # println("start")
        (f, g) = select!(P)
        # println("selected")
        r = mdiv(S(f, g), G)
        # println("division done")
        if r != []
            update!(P, G, r)
            append!(G, [r])
        end
        t += 1
        # println(t)
        # println(length(P))

    end
    return G
end

function is_groebner_basis(G)
    return all([mdiv(S(G[i], G[j]), G) == [] for i in 1:length(G) for j in i:length(G)])
end

function tuples_lt(n, k)
    if n == 0
        return Vector{Int64}[Int64[]]
    else
        tups = []
        for i in 0:k
            recs = tuples_lt(n - 1, k - i)
            for rec in recs
                push!(rec, i)
            end
            append!(tups, recs)
        end
        return tups
    end
end

# function poly(s, n)
#     ss = split(replace(s, r"\w" => ""), "+")
#     terms = match.(r"^(\d*)(.*)", ss)
#     for term in terms

#     sss = [SubString.(term, findall(r"\d*x\d+(\^\d+)?", term)) for term in ss]

# end


text:
module KdV
using ApproxFun, LinearAlgebra
export ReflectionCoefficient, tFun

struct ReflectionCoefficient{VM, VP} <: Function
    V₋::VM
    V₊::VP
    ReflectionCoefficient{VM, VP}(V₋::VM, V₊::VP) where {VM,VP} = new{VM,VP}(V₋, V₊)
end

function ReflectionCoefficient(V, a=-6.0, x₀=0.0, b=6.0)
    d₋, d₊ = a .. x₀ , x₀ .. b
    V₋,V₊ = Fun(V, d₋), Fun(V, d₊)
    ReflectionCoefficient{typeof(V₋),typeof(V₊)}(V₋,V₊)
end

# ψ'' + (V(x) + k^2) ψ = 0
function (R::ReflectionCoefficient)(k)
    k == 0 && return -one(ComplexF64)

    a,x₀ = endpoints(domain(R.V₋))
    b = rightendpoint(domain(R.V₊))
    D = Derivative()
    V₋,V₊ = R.V₋, R.V₊
    ψ = [ivp(); D^2  + (V₋ + k^2)] \ [exp(-im*k*a), -im*k*(exp(-im*k*a)), 0.0]

    F = qr([rdirichlet(space(V₊)); rneumann(); D^2  + (V₊ + k^2)])
    ϕ₊ = F \ [exp(im*k*b), im*k*(exp(im*k*b)), 0.0]
    ϕ₋ = F \ [exp(-im*k*b), -im*k*(exp(-im*k*b)), 0.0]

    a,b = [ϕ₋(x₀) ϕ₊(x₀);
           ϕ₋'(x₀) ϕ₊'(x₀) ] \ [ψ(x₀); ψ'(x₀)]
    b/a
end

# use multiple threads since reflection coefficient is slow
function tvalues(f, d, n)
    p = points(d, n)
    F = similar(p, ComplexF64)
    Threads.@threads for k=1:length(p)
        F[k] = f(p[k])
    end
    F
end

tFun(f, d::Space, n) = Fun(d, transform(d,tvalues(f,d,n)))
tFun(f, d, n) = tFun(f, Space(d), n)

function quickinv(F)
    A = Array(F)
    V = values.(pad.(A, maximum(ncoefficients.(A))))
    Vi = [inv([V[k,j][p] for k=1:2, j=1:2]) for p=1:length(V[1,1])]
    Fun(Fun([Fun(sp,transform(sp,[Vi[p][k,j] for p=1:length(Vi)])) for k=1:2, j=1:2]), space(F))
end
end # module


text:
function de_propose(rng::AbstractRNG, density, particles::AbstractVector, i::Int)
    γ = 2.38 / sqrt(2 * length(density)) * exp(randn(rng) * 0.1)
    a = b = i
    while a ∈ (i,)
        a = rand(rng, eachindex(particles))
    end
    while b ∈ (a, i)
        b = rand(rng, eachindex(particles))
    end
    W = op(*, op(-, particles[a], particles[b]), γ)
    T = op(
        x -> γ * x / 300 * randn(rng),
        op(
            +,
            op(abs, op(-, particles[a], particles[b])),
            op(abs, op(-, particles[i], particles[b])),
            op(abs, op(-, particles[a], particles[i])),
        ),
    )
    op(+, particles[i], W, T), 0.0
end

function ais_walk_propose(rng::AbstractRNG, density, particles::AbstractVector, i::Int)
    a = b = c = i
    while a ∈ (i,)
        a = rand(rng, eachindex(particles))
    end
    while b ∈ (a, i)
        b = rand(rng, eachindex(particles))
    end
    while c ∈ (b, a, i)
        c = rand(rng, eachindex(particles))
    end
    Xs = op(/, op(+, particles[a], op(+, particles[b], particles[c])), 3)
    W = op(
        +,
        op(*, randn(rng), op(-, particles[a], Xs)),
        op(*, randn(rng), op(-, particles[b], Xs)),
        op(*, randn(rng), op(-, particles[c], Xs)),
    )
    op(+, particles[i], W), 0.0
end

"Inverse cdf of g-pdf, see eq. 10 of Foreman-Mackey et al. 2013."
cdf_g_inv(u, a) = (u * (sqrt(a) - sqrt(1 / a)) + sqrt(1 / a))^2

"Sample from g using inverse transform sampling.  a=2.0 is recommended."
sample_g(rng::AbstractRNG, a) = cdf_g_inv(rand(rng), a)

function stretch_propose(rng::AbstractRNG, density, particles::AbstractVector, i::Int)
    a = i
    while i == a
        a = rand(rng, eachindex(particles))
    end
    Z = sample_g(rng, 3.0)
    W = op(*, op(-, particles[i], particles[a]), Z)
    op(+, particles[a], W), (length(density) - 1) * log(Z)
end

function propose(rng::AbstractRNG, density, particles::AbstractVector, i::Int)
    p = rand(rng, (1, 1, 1, 1, 2, 2, 3))
    pr = (stretch_propose, de_propose, ais_walk_propose)
    pr[p](rng, density, particles, i)
end

function transition!(
    density::AbstractDensity,
    particles::AbstractVector,
    logdensity::AbstractVector,
    particle_index::Int,
    rng::AbstractRNG,
)
    p, ld_correction = propose(rng, density, particles, particle_index)
    ld = loglike(density, push_p(density, p))
    if accept(density, rng, logdensity[particle_index], ld, ld_correction)
        particles[particle_index] = p
        logdensity[particle_index] = ld
        return true
    end
    return false
end


text:
import NodeModule.Node

function euclidean_distance(n1::Node, n2::Node)::Float32
    """Calculates the euclidean distance between given two nodes.
    """
    return sqrt((n1.x - n2.x)^2 + (n1.y - n2.y)^2)
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
using BEAST
using CompScienceMeshes
using Base.Test

x = Vec(1.0,0.0,0.0)
y = Vec(0.0,1.0,0.0)
z = Vec(0.0,0.0,1.0)

g = BEAST.creategaussian(1.0, 3.0)
pw = BEAST.planewave(x, z, g)

t = linspace(0.0, 6.0, 100);
y = [pw([0.0,0.0,0.0],x)[1] for x in t];

Δx = 1.0
Γ = meshrectangle(10.0, Δx, Δx)
X = raviartthomas(Γ)

Δt = 1.0
Nt = 20
U = timebasisc0d1(Float64, Δt, Nt)

b = zeros(Float64, numfunctions(X), numfunctions(U))
store(v, m, k) = (b[m,k] += v)
BEAST.assemble!(pw, X, U, store)


text:
#=
Code for simulating, visualizing and manipulating *scalar* additive and
multiplicative functionals.

@authors: Shunsuke Hori

=#

using QuantEcon

"""
This type and method are written to transform a scalar additive functional
into a linear state space system.
"""
type AMF_LSS_VAR{TR<:Real}
    A::TR
    B::TR
    D::TR
    F::TR
    nu::TR
    lss::LSS
end
function AMF_LSS_VAR(A::Real, B::Real, D::Real, F::Real=0.0, nu::Real=0.0)
    # Construct BIG state space representation
    lss = construct_ss(A,B,D,F,nu)
    return AMF_LSS_VAR(A,B,D,F,nu,lss)
end

"""
This creates the state space representation that can be passed
into the LSS type from QuantEcon.
"""
function construct_ss(A::Real, B::Real, D::Real, F::Real, nu::Real)
    H, g = additive_decomp(A, B, D, F)

    # Build A matrix for LSS
    # Order of states is: [1, t, xt, yt, mt]
    A1 = [1 0 0 0 0]       # Transition for 1
    A2 = [1 1 0 0 0]       # Transition for t
    A3 = [0 0 A 0 0]       # Transition for x_{t+1}
    A4 = [nu 0 D 1 0]      # Transition for y_{t+1}
    A5 = [0 0 0 0 1]       # Transition for m_{t+1}
    Abar = vcat(A1, A2, A3, A4, A5)

    # Build B matrix for LSS
    Bbar = [0, 0, B, F, H]

    # Build G matrix for LSS
    # Order of observation is: [xt, yt, mt, st, tt]
    G1 = [0 0 1 0 0]               # Selector for x_{t}
    G2 = [0 0 0 1 0]               # Selector for y_{t}
    G3 = [0 0 0 0 1]               # Selector for martingale
    G4 = [0 0 -g 0 0]              # Selector for stationary
    G5 = [0 nu 0 0 0]              # Selector for trend
    Gbar = vcat(G1, G2, G3, G4, G5)

    # Build LSS type
    x0 = [0, 0, 0, 0, 0]
    S0 = zeros(5, 5)
    lss = LSS(Abar, Bbar, Gbar, mu_0=x0, Sigma_0=S0)

    return lss
end

"""
Return values for the martingale decomposition (Proposition 4.3.3.)
    - `H`         : coefficient for the (linear) martingale component (kappa_a)
    - `g`         : coefficient for the stationary component g(x)
"""
function additive_decomp(A::Real,B::Real,D::Real,F::Real)
    A_res = 1 / (1 - A)
    g = D * A_res
    H = F + D * A_res * B

    return H, g
end

"""
Return values for the multiplicative decomposition (Example 5.4.4.)
    - `nu_tilde`  : eigenvalue
    - `H`         : vector for the Jensen term
"""
function multiplicative_decomp(A::Real, B::Real, D::Real, F::Real, nu::Real)
    H, g = additive_decomp(A, B, D, F)
    nu_tilde = nu + 0.5*H^2

    return nu_tilde, H, g
end

function loglikelihood_path(amf::AMF_LSS_VAR, x::Vector, y::Vector)
    A, B, D, F = amf.A, amf.B, amf.D, amf.F
    T = length(y)
    FF = F^2
    FFinv = 1/FF
    temp = y[2:end] - y[1:end-1] - D*x[1:end-1]
    obs =  temp .* FFinv .* temp
    obssum = cumsum(obs)
    scalar = (log(FF) + log(2pi))*collect(1:T-1)

    return (-0.5)*(obssum + scalar)
end

function loglikelihood(amf::AMF_LSS_VAR,
                       x::Vector, y::Vector)
    llh = loglikelihood_path(amf, x, y)
    return llh[end]
end


text:
# Here we override broadcasting for banded matrices.
# The design is to to exploit the broadcast machinery so that
# banded matrices that conform to the banded matrix interface but are not
# <: AbstractBandedMatrix can get access to fast copyto!, lmul!, rmul!, axpy!, etc.
# using broadcast variants (B .= A, B .= 2.0 .* A, etc.)

abstract type AbstractBlockBandedStyle <: AbstractArrayStyle{2} end
struct BandedBlockBandedStyle <: AbstractBlockBandedStyle end
abstract type AbstractBlockSkylineStyle <: AbstractBlockBandedStyle end
struct BlockSkylineStyle <: AbstractBlockSkylineStyle end
struct BlockBandedStyle <: AbstractBlockSkylineStyle end
BlockSkylineStyle(::Val{2}) = BlockBandedStyle()
BlockBandedStyle(::Val{2}) = BlockBandedStyle()
BandedBlockBandedStyle(::Val{2}) = BandedBlockBandedStyle()

BroadcastStyle(::BlockSkylineStyle, ::BandedBlockBandedStyle) = BlockSkylineStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::BlockSkylineStyle) = BlockSkylineStyle()
BroadcastStyle(::BlockSkylineStyle, ::BlockBandedStyle) = BlockSkylineStyle()
BroadcastStyle(::BlockBandedStyle, ::BlockSkylineStyle) = BlockSkylineStyle()

BroadcastStyle(::BlockBandedStyle, ::BandedBlockBandedStyle) = BlockBandedStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::BlockBandedStyle) = BlockBandedStyle()

BroadcastStyle(::DefaultArrayStyle{2}, ::BlockSkylineStyle) = BlockSkylineStyle()
BroadcastStyle(::BlockSkylineStyle, ::DefaultArrayStyle{2}) = BlockSkylineStyle()

BroadcastStyle(::DefaultArrayStyle{2}, ::BlockBandedStyle) = BlockBandedStyle()
BroadcastStyle(::BlockBandedStyle, ::DefaultArrayStyle{2}) = BlockBandedStyle()

BroadcastStyle(::DefaultArrayStyle{2}, ::BandedBlockBandedStyle) = BandedBlockBandedStyle()
BroadcastStyle(::BandedBlockBandedStyle, ::DefaultArrayStyle{2}) = BandedBlockBandedStyle()

###
# broadcast blockbandwidths
###
_blockbnds(bc) = blocksize(bc) .- 1

_blockbandwidth_u(u, a, ax) = blockisequal(axes(a,1),ax) ? (blockbandwidth(a,1),u) : (blocksize(ax,1)-1,u)
_blockbandwidth_l(l, a, ax) = blockisequal(axes(a,2),ax) ? (l,blockbandwidth(a,2)) : (l,blocksize(ax,1)-1)

_broadcast_blockbandwidths(bnds, _, _) = bnds
_broadcast_blockbandwidths((l,u), a::AbstractVector, (ax1,ax2)) =
    _blockbandwidth_u(u, a, ax1)

function _broadcast_blockbandwidths((l,u), A::AbstractArray, (ax1,ax2))
    if size(A,2) == 1
        _blockbandwidth_u(u, A, ax1)
    elseif size(A,1) == 1
        _blockbandwidth_l(l, A, ax2)
    elseif blockisequal(axes(A), (ax1,ax2))
        blockbandwidths(A) # need to special case vector broadcasting
    else
        _blockbnds(A)
    end
end

blockbandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(*)}) =
    min.(_broadcast_blockbandwidths.(Ref(_blockbnds(bc)), bc.args, Ref(axes(bc)))...)

blockbandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(/)}) = _broadcast_blockbandwidths(_blockbnds(bc), first(bc.args), axes(bc))
blockbandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(\)}) = _broadcast_blockbandwidths(_blockbnds(bc), last(bc.args), axes(bc))

_subblockbnds(bc) = size(bc) .- 1

_broadcast_subblockbandwidths(bnds) = bnds
_broadcast_subblockbandwidths(bnds, _) = bnds
_broadcast_subblockbandwidths((l,u), a::AbstractVector) = (bandwidth(a,1),u)

function _broadcast_subblockbandwidths((l,u), A::AbstractArray)
    size(A,2) == 1 && return (subblockbandwidth(A,1),u)
    size(A,1) == 1 && return (l, subblockbandwidth(A,2))
    subblockbandwidths(A) # need to special case vector broadcasting
end

subblockbandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(*)}) =
    min.(_broadcast_subblockbandwidths.(Ref(_subblockbnds(bc)), bc.args)...)

subblockbandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(/)}) = _broadcast_subblockbandwidths(_subblockbnds(bc), first(bc.args))
subblockbandwidths(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Any,typeof(\)}) = _broadcast_subblockbandwidths(_subblockbnds(bc), last(bc.args))

function subblockbandwidths(bc::Broadcasted)
    (a,b) = size(bc)
    bnds = (a-1,b-1)
    _isweakzero(bc.f, bc.args...) && return min.(bnds, max.(_broadcast_subblockbandwidths.(Ref(bnds), bc.args)...))
    bnds
end

# zero is preserved. Take the maximum bandwidth
import BandedMatrices: _isweakzero

function blockbandwidths(bc::Broadcasted)
    (a,b) = size(bc)
    bnds = (a-1,b-1)
    _isweakzero(bc.f, bc.args...) && return min.(bnds, max.(_broadcast_blockbandwidths.(Ref(bnds), bc.args, Ref(axes(bc)))...))
    bnds
end

similar(bc::Broadcasted{BlockBandedStyle}, ::Type{T}) where T = BlockBandedMatrix{T}(undef, axes(bc), blockbandwidths(bc))
similar(bc::Broadcasted{BandedBlockBandedStyle}, ::Type{T}) where T = BandedBlockBandedMatrix{T}(undef, axes(bc), blockbandwidths(bc), subblockbandwidths(bc))

####
# Default to standard Array broadcast
#
# This is because, for example, exp.(B) is not banded
####

copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle}) =
   copyto!(dest, Broadcasted{DefaultArrayStyle{2}}(bc.f, bc.args, bc.axes))

similar(bc::Broadcasted{<:AbstractBlockBandedStyle}, ::Type{T}) where T =
    similar(Broadcasted{DefaultArrayStyle{2}}(bc.f, bc.args, bc.axes), T)

##
# copyto!
##

checkblocks(A, B) = blockisequal(axes(A), axes(B)) || throw(DimensionMismatch("*"))

function _blockbanded_copyto!(dest::AbstractMatrix, src::AbstractMatrix)
    @boundscheck checkblocks(dest, src)

    dl, du = colblockbandwidths(dest)
    sl, su = colblockbandwidths(src)
    M,N = blocksize(src)
    # Source matrix must fit within bands of destination matrix
    all(dl .≥ min.(sl,Ref(M-1))) && all(du .≥ min.(su,Ref(N-1))) || throw(BandError(dest))

    for J = 1:N
        for K = max(1,J-du[J]):min(J-su[J]-1,M)
            zero!(view(dest,Block(K),Block(J)))
        end
        for K = max(1,J-su[J]):min(J+sl[J],M)
            copyto!(view(dest,Block(K),Block(J)), view(src,Block(K),Block(J)))
        end
        for K = max(1,J+sl[J]+1):min(J+dl[J],M)
            zero!(view(dest,Block(K),Block(J)))
        end
    end
    dest
end

function blockbanded_copyto!(dest::AbstractMatrix, src::AbstractMatrix)
    if isblockbanded(dest)
        _blockbanded_copyto!(dest, src)
    else
        _blockbanded_copyto!(PseudoBlockArray(dest, axes(src)), src)
    end
    dest
end

_copyto!(_, ::AbstractBlockBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) = blockbanded_copyto!(dest, src)
function _copyto!(_, ::BlockLayout{<:BandedColumns}, dest::AbstractMatrix, src::AbstractMatrix)
    if !blockisequal(axes(dest), axes(src))
        copyto!(PseudoBlockArray(dest, axes(src)), src)
        return dest
    end

    srcB = blocks(src)
    srcD = bandeddata(srcB)

    dl, du = colblockbandwidths(dest)
    sl, su = bandwidths(srcB)
    M,N = size(srcB)
    # Source matrix must fit within bands of destination matrix
    all(dl .≥ min(sl,M-1)) && all(du .≥ min(su,N-1)) || throw(BandError(dest))

    for J = 1:N
        for K = max(1,J-du[J]):min(J-su-1,M)
            zero!(view(dest,Block(K),Block(J)))
        end
        for K = max(1,J-su):min(J+sl,M)
            copyto!(view(dest,Block(K),Block(J)), srcD[K-J+su+1,J])
        end
        for K = max(1,J+sl+1):min(J+dl[J],M)
            zero!(view(dest,Block(K),Block(J)))
        end
    end
    dest
end

function _copyto!(_, ::BlockLayout{<:AbstractBandedLayout}, dest::AbstractMatrix, src::AbstractMatrix)
    if !blockisequal(axes(dest), axes(src))
        copyto!(PseudoBlockArray(dest, axes(src)), src)
        return dest
    end

    srcB = blocks(src)

    dl, du = colblockbandwidths(dest)
    sl, su = bandwidths(srcB)
    M,N = size(srcB)
    # Source matrix must fit within bands of destination matrix
    all(dl .≥ min(sl,M-1)) && all(du .≥ min(su,N-1)) || throw(BandError(dest))

    for J = 1:N
        for K = max(1,J-du[J]):min(J-su-1,M)
            zero!(view(dest,Block(K),Block(J)))
        end
        for K = max(1,J-su):min(J+sl,M)
            copyto!(view(dest,Block(K),Block(J)), inbands_getindex(srcB, K, J))
        end
        for K = max(1,J+sl+1):min(J+dl[J],M)
            zero!(view(dest,Block(K),Block(J)))
        end
    end
    dest
end

function _copyto!(::BandedBlockBandedColumns, ::BandedBlockBandedColumns, dest::AbstractMatrix, src::AbstractMatrix)
    if blockbandwidths(dest) == blockbandwidths(src) && subblockbandwidths(dest) == subblockbandwidths(src)
        copyto!(bandedblockbandeddata(dest), bandedblockbandeddata(src))
        dest
    else # TODO: if subblockbandwidths match we can still do something better
        _copyto!(BandedBlockBandedLayout(), BandedBlockBandedLayout(), dest, src)
    end
end

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(identity)})
    (A,) = bc.args
    dest ≡ A && return dest
    blockbanded_copyto!(dest, A)
end

##
# lmul!/rmul!
##

######
# extra marrix routines
#####

function blockbanded_fill!(B::AbstractMatrix{T}, x) where T
    x == zero(T) || throw(BandError(B))

    for J = blockaxes(B,2), K = blockcolsupport(B,J)
        fill!(view(B,K,J), x)
    end
    B
end

function blockbanded_rmul!(B::AbstractMatrix{T}, x::Number) where T
    x == zero(T) || throw(BandError(B))

    for J = blockaxes(B,2), K = blockcolsupport(B,J)
        rmul!(view(B,K,J), x)
    end
    B
end

function blockbanded_lmul!(x::Number, B::AbstractMatrix{T}) where T
    x == zero(T) || throw(BandError(B))

    for J = blockaxes(B,2), K = blockcolsupport(B,J)
        lmul!(x, view(B,K,J))
    end
    B
end

function blockbanded_fill!(B::AbstractBlockBandedMatrix, x::Number)
    x == zero(T) || throw(BandError(B))

    fill!(B.data, x)
end

function blockbanded_rmul!(A::AbstractBlockBandedMatrix, x::Number)
    rmul!(A.data, x)
    A
end

function blockbanded_lmul!(x::Number, A::AbstractBlockBandedMatrix)
    lmul!(x, A.data)
    A
end

lmul!(x::Number, A::AbstractBlockBandedMatrix) = blockbanded_lmul!(x, A)
rmul!(A::AbstractBlockBandedMatrix, x::Number) = blockbanded_rmul!(A, x)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(-), <:Tuple{<:AbstractMatrix}})
    A, = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_lmul!(-1, dest)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(-), <:Tuple{<:AbstractMatrix}}, ::Type{T}) where T =
    similar(bc.args[1], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:Number,<:AbstractMatrix}})
    α,A = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_lmul!(α, dest)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:Number,<:AbstractMatrix}}, ::Type{T}) where T =
    similar(bc.args[2], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:AbstractMatrix,<:Number}})
    A,α = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_rmul!(dest, α)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(*), <:Tuple{<:AbstractMatrix,<:Number}}, ::Type{T}) where T =
    similar(bc.args[1], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(/), <:Tuple{<:AbstractMatrix,<:Number}})
    A,α = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_rmul!(dest, inv(α))
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(/), <:Tuple{<:AbstractMatrix,<:Number}}, ::Type{T}) where T =
    similar(bc.args[1], T)

function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(\), <:Tuple{<:Number,<:AbstractMatrix}})
    α,A = bc.args
    dest ≡ A || copyto!(dest, A)
    blockbanded_lmul!(inv(α), dest)
end

similar(bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(\), <:Tuple{<:Number,<:AbstractMatrix}}, ::Type{T}) where T =
    similar(bc.args[2], T)

##
# axpy!
##

# these are the routines of the banded interface of other AbstractMatrices
function blockbanded_axpy!(a, X::AbstractMatrix, Y::AbstractMatrix)
    size(X) == size(Y) || throw(DimensionMismatch())

    for J = blockaxes(X,2), K = blockcolsupport(X,J)
        view(Y,K,J) .= a .* view(X,K,J) .+ view(Y,K,J)
    end
    Y
end

for op in (:+, :-)
    @eval function copyto!(C::AbstractArray{T}, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof($op),
                                                                <:Tuple{<:AbstractMatrix,<:AbstractMatrix}}) where T
        bc_axes = Base.Broadcast.combine_axes(bc.args...)# Use combine_axes as `bc` might get axes from `C`
        if !blockisequal(axes(C), bc_axes)
            copyto!(PseudoBlockArray(C, bc_axes), bc)
            return C
        end

        A,B = bc.args
        if !blockisequal(axes(A), axes(B))
            copyto!(C, Base.broadcasted(BlockStyle{2}(), $op, A, B))
            return C
        end

        A_l,A_u = blockbandwidths(A)
        B_l,B_u = blockbandwidths(B)
        C_l,C_u = blockbandwidths(C)

        N,M = blocksize(C)

        for J̃ = 1:M
            J = Block(J̃)
            # joint block range
            for K = Block.(max(1,J̃-min(A_u,B_u)):min(N,J̃+min(A_l,B_l)))
                view(C,K,J) .= $op.(view(A,K,J), view(B,K,J))
            end
            # above, A but not B
            for K = Block.(max(1,J̃-A_u):min(N,J̃-B_u-1))
                view(C,K,J) .= view(A,K,J)
            end
            # above, B but not A
            for K = Block.(max(1,J̃-B_u):min(N,J̃-A_u-1))
                view(C,K,J) .= $op.(view(B,K,J))
            end
            # below, A but not B
            for K = Block.(max(1,J̃+B_l+1):min(N,J̃+A_l))
                view(C,K,J) .= view(A,K,J)
            end
            # below, B but not A
            for K = Block.(max(1,J̃+A_l+1):min(N,J̃+B_l))
                view(C,K,J) .= $op.(view(B,K,J))
            end
            # above, neither
            for K = Block.(max(J̃-C_u,1):min(J̃-max(A_u,B_u)-1,N))
                view(C,K,J) .= zero(T)
            end
            # below, neither
            for K = Block.(max(J̃+max(A_l,B_l)+1,1):min(J̃+C_l,N))
                view(C,K,J) .= zero(T)
            end
        end

        C
    end
end

##
# special cases to hack tests where 2.0 .* A .+ B can be constructed in different ways depending on REPL or in a function.
# Should be made more general.
##

function copyto!(dest::AbstractArray{T}, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(+),
                                                        <:Tuple{<:Broadcasted{<:AbstractBlockBandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                                                        <:AbstractMatrix}}) where T
    αA,B = bc.args
    α,A = αA.args
    dest ≡ B || blockbanded_copyto!(dest, B)
    blockbanded_axpy!(α, A, dest)
end

function similar(bc::Broadcasted{BlockBandedStyle, <:Any, typeof(+),
                        <:Tuple{<:Broadcasted{<:AbstractBlockBandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                        <:AbstractMatrix}}, ::Type{T}) where T
    αA,B = bc.args
    α,A = αA.args
    Al,Au = blockbandwidths(A)
    Bl,Bu = blockbandwidths(B)
    BlockBandedMatrix{T}(undef, axes(bc), (max(Al,Bl), max(Au,Bu)))
end

function similar(bc::Broadcasted{BandedBlockBandedStyle, <:Any, typeof(+),
                        <:Tuple{<:Broadcasted{<:AbstractBlockBandedStyle,<:Any,typeof(*),<:Tuple{<:Number,<:AbstractMatrix}},
                        <:AbstractMatrix}}, ::Type{T}) where T
    αA,B = bc.args
    α,A = αA.args
    Al,Au = blockbandwidths(A)
    Bl,Bu = blockbandwidths(B)
    Aλ,Aμ = subblockbandwidths(A)
    Bλ,Bμ = subblockbandwidths(B)

    BandedBlockBandedMatrix{T}(undef, axes(bc), (max(Al,Bl), max(Au,Bu)), (max(Aλ,Bλ), max(Aμ,Bμ)))
end

####
# Special case for Diagonal multiplicartion
####


text:
using BinaryBuilder

name = "MbedTLS"
version = v"2.24.0"

# Collection of sources required to build MbedTLS
sources = [
    GitSource("https://github.com/ARMmbed/mbedtls.git",
              "523f0554b6cdc7ace5d360885c3f5bbcc73ec0e8"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/mbedtls

atomic_patch -p1 ../patches/fix_incorrect_EOF_check.patch

# llvm-ranlib gets confused, use the binutils one
if [[ "${target}" == *apple* ]]; then
    ln -sf /opt/${target}/bin/${target}-ranlib /opt/bin/ranlib
    ln -sf /opt/${target}/bin/${target}-ranlib /opt/bin/${target}-ranlib
fi

# enable MD4
sed "s|//#define MBEDTLS_MD4_C|#define MBEDTLS_MD4_C|" -i include/mbedtls/config.h

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DCMAKE_C_STANDARD=99 \
    -DUSE_SHARED_MBEDTLS_LIBRARY=On \
    ..
make -j${nproc} && make install

if [[ "${target}" == *mingw* ]]; then
    # For some reason, the build system doesn't set the `.dll` files as
    # executable, which prevents them from being loaded.  Also, we need
    # to explicitly use `${prefix}/lib` here because the build system
    # is a simple one, and blindly uses `/lib`, even on Windows.
    chmod +x ${prefix}/lib/*.dll
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libmbedcrypto", :libmbedcrypto),
    LibraryProduct("libmbedx509", :libmbedx509),
    LibraryProduct("libmbedtls", :libmbedtls),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: vchuravy/Yggdrasil

max_issues_repo_name: vchuravy/Yggdrasil

max_forks_repo_name: vchuravy/Yggdrasil