text:
"""
Abstract atmospheric conditions type. The suptypes of AbstractAtmosphere should describe the
atmospheric conditions for one time-step only, see *e.g.* [`Atmosphere`](@ref)
"""
abstract type AbstractAtmosphere end

"""
Atmosphere structure to hold all values related to the meteorology / atmosphere.

# Arguments

- `date = Dates.now()`: the date of the record.
- `duration = 1.0` (seconds): the duration of the time-step.
- `T` (°C): air temperature
- `Wind` (m s-1): wind speed
- `P = 101.325` (kPa): air pressure. The default value is at 1 atm, *i.e.* the mean sea-level
atmospheric pressure on Earth.
- `Rh = rh_from_vpd(VPD,eₛ)` (0-1): relative humidity
- `Cₐ` (ppm): air CO₂ concentration
- `e = vapor_pressure(T,Rh)` (kPa): vapor pressure
- `eₛ = e_sat(T)` (kPa): saturated vapor pressure
- `VPD = eₛ - e` (kPa): vapor pressure deficit
- `ρ = air_density(T, P, constants.Rd, constants.K₀)` (kg m-3): air density
- `λ = latent_heat_vaporization(T, constants.λ₀)` (J kg-1): latent heat of vaporization
- `γ = psychrometer_constant(P, λ, constants.Cₚ, constants.ε)` (kPa K−1): psychrometer "constant"
- `ε = atmosphere_emissivity(T,e,constants.K₀)` (0-1): atmosphere emissivity
- `Δ = e_sat_slope(meteo.T)` (0-1): slope of the saturation vapor pressure at air temperature
- `clearness::A = 9999.9` (0-1): Sky clearness
- `Ri_SW_f::A = 9999.9` (W m-2): Incoming short wave radiation flux
- `Ri_PAR_f::A = 9999.9` (W m-2): Incoming PAR radiation flux
- `Ri_NIR_f::A = 9999.9` (W m-2): Incoming NIR radiation flux
- `Ri_TIR_f::A = 9999.9` (W m-2): Incoming TIR radiation flux
- `Ri_custom_f::A = 9999.9` (W m-2): Incoming radiation flux for a custom waveband

# Notes

The structure can be built using only `T`, `Rh`, `Wind` and `P`. All other variables are optional
and can be automatically computed using the functions given in `Arguments`.

# Examples

```julia
Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65)
```
"""
struct Atmosphere{A,D1,D2} <: AbstractAtmosphere
    date::D1
    duration::D2
    T::A
    Wind::A
    P::A
    Rh::A
    Cₐ::A
    e::A
    eₛ::A
    VPD::A
    ρ::A
    λ::A
    γ::A
    ε::A
    Δ::A
    clearness::A
    Ri_SW_f::A
    Ri_PAR_f::A
    Ri_NIR_f::A
    Ri_TIR_f::A
    Ri_custom_f::A
end

function Atmosphere(;
    date = Dates.now(), duration = 1.0, T, Wind, P = 101.325, Rh,
    Cₐ = 400.0, e = vapor_pressure(T, Rh), eₛ = e_sat(T), VPD = eₛ - e,
    ρ = air_density(T, P), λ = latent_heat_vaporization(T),
    γ = psychrometer_constant(P, λ), ε = atmosphere_emissivity(T, e),
    Δ = e_sat_slope(T), clearness = 9999.9, Ri_SW_f = 9999.9, Ri_PAR_f = 9999.9,
    Ri_NIR_f = 9999.9, Ri_TIR_f = 9999.9, Ri_custom_f = 9999.9)

        # Checking some values:
    if Wind <= 0
        @warn "Wind should always be > 0, forcing it to 1e-6"
        Wind = 1e-6
    end

    if Rh <= 0
        @warn "Rh should always be > 0, forcing it to 1e-6"
        Rh = 1e-6
    end

    if Rh > 1
        if 1 < Rh < 100
            @warn "Rh should be 0 < Rh < 1, assuming it is given in % and dividing by 100"
            Rh /= 100
        else
            @error "Rh should be 0 < Rh < 1, and its value is $(Rh)"
        end
    end

    if clearness != 9999.9 && (clearness <= 0 || clearness > 1)
        @error "clearness should always be 0 < clearness < 1"
    end

    param_A =
    promote(
        T, Wind, P, Rh, Cₐ, e, eₛ, VPD, ρ, λ, γ, ε, Δ, clearness, Ri_SW_f, Ri_PAR_f,
        Ri_NIR_f, Ri_TIR_f, Ri_custom_f
    )

    Atmosphere(date, duration, param_A...)
end

"""
    Weather(D <: AbstractArray{<:AbstractAtmosphere}[, S])
    Weather(df::DataFrame[, mt])

Defines the weather, *i.e.* the local conditions of the Atmosphere for one or more time-steps.
Each time-step is described using the [`Atmosphere`](@ref) structure.

The simplest way to instantiate a `Weather` is to use a `DataFrame` as input.

The `DataFrame` should be formated such as each row is an observation for a given time-step
and each column is a variable. The column names should match exactly the field names of the
[`Atmosphere`](@ref) structure, `i.e.`:

```@example
fieldnames(Atmosphere)
```

## See also

- the [`Atmosphere`](@ref) structure
- the [`read_weather`](@ref) function to read Archimed-formatted meteorology data.

## Examples

```julia
# Example of weather data defined by hand (cumbersome):
w = Weather(
    [
        Atmosphere(T = 20.0, Wind = 1.0, P = 101.3, Rh = 0.65),
        Atmosphere(T = 23.0, Wind = 1.5, P = 101.3, Rh = 0.60),
        Atmosphere(T = 25.0, Wind = 3.0, P = 101.3, Rh = 0.55)
    ],
    (
        site = "Test site",
        important_metadata = "this is important and will be attached to our weather data"
    )
)

# Example using a DataFrame, that you would usually import from a file:
using CSV, DataFrames
file = joinpath(dirname(dirname(pathof(PlantBiophysics))),"test","inputs","meteo.csv")
df = CSV.read(file, DataFrame; header=5, datarow = 6)
# Select and rename the variables:
select!(df, :date, :VPD, :temperature => :T, :relativeHumidity => :Rh, :wind => :Wind, :atmosphereCO2_ppm => :Cₐ)
df[!,:duration] .= 1800 # Add the time-step duration, 30min

# Make the weather, and add some metadata:
Weather(df, (site = "Aquiares", file = file))
```
"""
struct Weather{D <: AbstractArray{<:AbstractAtmosphere},S <: MutableNamedTuple}
    data::D
    metadata::S
end

function Weather(df::T) where T <: AbstractArray{<:AbstractAtmosphere}
    Weather(df, MutableNamedTuple())
end

function Weather(df::T, mt::S) where {T <: AbstractArray{<:AbstractAtmosphere},S <: NamedTuple}
    Weather(df, MutableNamedTuple(;mt...))
end

function Weather(df::DataFrame, mt::S) where S <: MutableNamedTuple
    Weather([Atmosphere(; i...) for i in eachrow(df)], mt)
end

function Weather(df::DataFrame, mt::S) where S <: NamedTuple
    mt = MutableNamedTuple(;mt...)
    Weather([Atmosphere(; i...) for i in eachrow(df)], mt)
end

function Weather(df::DataFrame, dict::S) where S <: AbstractDict
    # There must be a better way for transforming a Dict into a MutableNamedTuple...
    Weather(df, MutableNamedTuple(; NamedTuple{Tuple(Symbol.(keys(dict)))}(values(dict))...))
end

function Weather(df::DataFrame)
    Weather(df, MutableNamedTuple())
end

function Base.show(io::IO, n::Weather)
    printstyled(io, "Weather data.\n", bold = true, color = :green)
    printstyled(io, "Metadata: `$(NamedTuple(n.metadata))`.\n", color = :cyan)
    printstyled(io, "Data:\n", color = :green)
    # :normal, :default, :bold, :black, :blink, :blue, :cyan, :green, :hidden, :light_black, :light_blue, :light_cyan, :light_green, :light_magenta, :light_red, :light_yellow, :magenta, :nothing, :red,
#   :reverse, :underline, :white, or :yellow
    print(DataFrame(n))
    return nothing
end

Base.getindex(w::Weather, i::Integer) = w.data[i]
Base.getindex(w::Weather, s::Symbol) = [getproperty(i,s) for i in w.data]

"""
    DataFrame(data::Weather)

Transform a Weather type into a DataFrame.

See also [`Weather`](@ref) to make the reverse.
"""
function DataFrame(data::Weather)
return DataFrame(data.data)
end

"""
    vapor_pressure(Tₐ, rh)

Vapor pressure (kPa) at given temperature (°C) and relative hunidity (0-1).
"""
function vapor_pressure(Tₐ, rh)
    rh * e_sat(Tₐ)
end

"""
    e_sat(T)

Saturated water vapour pressure (es, in kPa) at given temperature `T` (°C).
See Jones (1992) p. 110 for the equation.
"""
function e_sat(T)
0.61375 * exp((17.502 * T) / (T + 240.97))
end

"""
    e_sat_slope(T)

Slope of the vapor pressure saturation curve at a given temperature `T` (°C).
"""
function e_sat_slope(T)
  (e_sat(T + 0.1) - e_sat(T)) / 0.1
end

"""
    air_density(Tₐ, P)
    air_density(Tₐ, P, Rd, K₀)

ρ, the air density (kg m-3).

# Arguments

- `Tₐ` (Celsius degree): air temperature
- `P` (kPa): air pressure
- `Rd` (J kg-1 K-1): gas constant of dry air (see Foken p. 245, or R bigleaf package).
- `K₀` (Celsius degree): temperature in Celsius degree at 0 Kelvin

# Note

Rd and K₀ are Taken from [`Constants`](@ref) if not provided.

# References

Foken, T, 2008: Micrometeorology. Springer, Berlin, Germany.
"""
function air_density(Tₐ, P, Rd, K₀)
    (P * 1000) / (Rd * (Tₐ - K₀))
end

function air_density(Tₐ, P)
    constants = Constants()
    air_density(Tₐ, P, constants.Rd, constants.K₀)
end

"""
    psychrometer_constant(P, λ, Cₚ, ε)
    psychrometer_constant(P, λ)

γ, the psychrometer constant, also called psychrometric constant (kPa K−1). See Monteith and
Unsworth (2013), p. 222.

# Arguments

- `P` (kPa): air pressure
- `λ` (``J\\ kg^{-1}``): latent heat of vaporization for water (see [`latent_heat_vaporization`](@ref))
- `Cₚ` (J kg-1 K-1): specific heat of air at constant pressure (``J\\ K^{-1}\\ kg^{-1}``)
- `ε` (Celsius degree): temperature in Celsius degree at 0 Kelvin

# Note

Cₚ, ε and λ₀ are taken from [`Constants`](@ref) if not provided.

```julia
Tₐ = 20.0

λ = latent_heat_vaporization(Tₐ, λ₀)
psychrometer_constant(100.0, λ)
```

# References

Monteith, John L., et Mike H. Unsworth. 2013. « Chapter 13 - Steady-State Heat Balance: (i)
Water Surfaces, Soil, and Vegetation ». In Principles of Environmental Physics (Fourth Edition),
edited by John L. Monteith et Mike H. Unsworth, 217‑47. Boston: Academic Press.

"""
function psychrometer_constant(P, λ, Cₚ, ε)
    γ = (Cₚ * P) / (ε * λ)
    return γ
end

function psychrometer_constant(P, λ)
    constant = Constants()
    γ = (constant.Cₚ * P) / (constant.ε * λ)
    return γ
end

"""
γ_star(γ, a_sh, a_s, rbv, Rsᵥ, Rbₕ)

γ∗, the apparent value of psychrometer constant (kPa K−1).

# Arguments

- `γ` (kPa K−1): psychrometer constant
- `aₛₕ` (1,2): number of faces exchanging heat fluxes (see Schymanski et al., 2017)
- `aₛᵥ` (1,2): number of faces exchanging water fluxes (see Schymanski et al., 2017)
- `Rbᵥ` (s m-1): boundary layer resistance to water vapor
- `Rsᵥ` (s m-1): stomatal resistance to water vapor
- `Rbₕ` (s m-1): boundary layer resistance to heat

# Note

Using the corrigendum from Schymanski et al. (2017) in here so the definition of
[`latent_heat`](@ref) remains generic.

Not to be confused with [`Γ_star`](@ref) in FcVB model, which is the CO₂ compensation point.

# References

Monteith, John L., et Mike H. Unsworth. 2013. « Chapter 13 - Steady-State Heat Balance: (i)
Water Surfaces, Soil, and Vegetation ». In Principles of Environmental Physics (Fourth Edition),
edited by John L. Monteith et Mike H. Unsworth, 217‑47. Boston: Academic Press.

Schymanski, Stanislaus J., et Dani Or. 2017. Leaf-Scale Experiments Reveal an Important
Omission in the Penman–Monteith Equation ». Hydrology and Earth System Sciences 21 (2): 685‑706.
https://doi.org/10.5194/hess-21-685-2017.
"""
function γ_star(γ, aₛₕ, aₛᵥ, Rbᵥ, Rsᵥ, Rbₕ)
    γ * aₛₕ / aₛᵥ * (Rbᵥ + Rsᵥ) / Rbₕ # rv + Rsᵥ= Boundary + stomatal conductance to water vapour
end

"""
    latent_heat_vaporization(Tₐ,λ₀)
    latent_heat_vaporization(Tₐ)

λ, the latent heat of vaporization for water (J kg-1).

# Arguments

- `Tₐ` (°C): air temperature
- `λ₀`: latent heat of vaporization for water at 0 degree Celsius. Taken from `Constants().λ₀`
if not provided (see [`Constants`](@ref)).

"""
function latent_heat_vaporization(Tₐ, λ₀)
  λ₀ - 2.365e3 * Tₐ
end

function latent_heat_vaporization(Tₐ)
    latent_heat_vaporization(Tₐ, Constants().λ₀)
end


text:
@testset "6 Hyperbolic functions" begin
    include("6.1 Hyperbolic sine/.jl")
    include("6.2 Hyperbolic cosine/.jl")
    include("6.3 Hyperbolic tangent/.jl")
    include("6.4 Hyperbolic cotangent/.jl")
    include("6.5 Hyperbolic secant/.jl")
    include("6.6 Hyperbolic cosecant/.jl")
    include("6.7 Miscellaneous/.jl")
end

ext: jl

max_stars_repo_path: test/6 Hyperbolic functions/.jl

max_issues_repo_path: test/6 Hyperbolic functions/.jl

max_forks_repo_path: test/6 Hyperbolic functions/.jl


text:
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np
@pyimport json
@pyimport os
@pyimport pickle
@pyimport pandas as pd
@pyimport collections

@pyimport parameters_julia
@pyimport utils_julia

out_dir = parameters_julia.out_dir;
files_ID = parameters_julia.files_ID;
month_w = parameters_julia.month_w;
month = parameters_julia.month;
year = parameters_julia.year;
instances_1 = parameters_julia.instances_ID;
deg_grid = parameters_julia.deg_grid;
c_grid = parameters_julia.c_grid;
lamb_grid = parameters_julia.lamb_grid;
week_day_Apr_list = parameters_julia.week_day_list;

@pyimport GLS_julia
@pyimport Compute_Jacobian
@pyimport get_observed_flow

numNodes = Compute_Jacobian.numNodes;
numLinks = Compute_Jacobian.numLinks;
numODpairs = Compute_Jacobian.numODpairs;
numZones = Compute_Jacobian.numZones
od_pairs = Compute_Jacobian.od_pairs;
link_list_js = Compute_Jacobian.link_list_js;
link_length_list = Compute_Jacobian.link_length_list;

function read_demand_file(file_path, n_nodes);
    file = open(file_path)
    demands = Dict()
    n = n_nodes  # number of nodes
    for i = 1:n
        demands[(i,i)] = 0.0
    end
    for line in eachline(file)
        OD_demand = split(line, ",")
        key, value = (parse(Int, OD_demand[1]), parse(Int, OD_demand[2])), parse(Float64, split(OD_demand[3], "\n")[1])
        demands[key] = value
    end
    close(file)
    return demands
end

# Sioux Falls network data
# http://www.bgu.ac.il/~bargera/tntp/

#Link travel time = free flow time * ( 1 + B * (flow/capacity)^Power ).
#Link generalized cost = Link travel time + toll_factor * toll + distance_factor * distance

# Traffic Assignment Data structure
type TA_Data
    network_name::String

    number_of_zones::Int64
    number_of_nodes::Int64
    first_thru_node::Int64
    number_of_links::Int64

    start_node::Array
    end_node::Array
    capacity::Array
    link_length::Array
    free_flow_time::Array
    B::Array
    power::Array
    speed_limit::Array
    toll::Array
    link_type::Array

    total_od_flow::Float64

    travel_demand::Array
    od_pairs::Array

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64
end

function load_ta_network_(out_dir, files_ID, month_w, day,  instan)
    toll_factor = 0
    distance_factor = 0
    network_data_file = files_ID * "_net_" * month_w * "_" * string(day) * '_' * instan * ".txt"
    trip_table_file = files_ID * "_trips_" * month_w * "_" * string(day) * '_' * instan * ".txt"
    best_objective = 0

    network_name = files_ID

    network_data_file =  out_dir  * "data_traffic_assignment_uni-class/" * network_data_file
    trip_table_file =   out_dir  * "data_traffic_assignment_uni-class/" * trip_table_file

    number_of_zones = 0
    number_of_links = 0
    number_of_nodes = 0
    first_thru_node = 0

    n = 0
    A = 0
    n = open(network_data_file, "r")
    A = readlines(n)
    for line in A
        #print( line)
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<NUMBER OF NODES>")
            number_of_nodes = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<FIRST THRU NODE>")
            first_thru_node = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<NUMBER OF LINKS>")
            number_of_links = parse( Int, line[ search(line, '>')+1 : end] )
        elseif contains(line, "<END OF METADATA>")
       # println(line)
        end
    end

    @assert number_of_links > 0
    start_node = Array{Int}(number_of_links)
    end_node = Array{Int}(number_of_links)
    capacity = zeros(number_of_links)
    link_length = zeros(number_of_links)
    free_flow_time = zeros(number_of_links)
    B = zeros(number_of_links)
    power = zeros(number_of_links)
    speed_limit = zeros(number_of_links)
    toll = zeros(number_of_links)
    link_type = Array{Int}(number_of_links)

    idx = 1

    apa = 0
    for line in A
        if contains(line, "~")
            continue
        end

        if contains(line, ";")
            line = strip(line, '\n')
            line = strip(line, ';')
            numbers = split(line)
            start_node[idx] = parse(Int, numbers[1])
            end_node[idx] = parse(Int, numbers[2])

            capacity[idx] = parse(Float64, numbers[3])
            link_length[idx] = parse(Float64, numbers[4])

            free_flow_time[idx] = parse(Float64, numbers[5])
            B[idx] = parse(Float64, numbers[6])
            power[idx] = parse(Float64, numbers[7])
            speed_limit[idx] = parse(Float64, numbers[8])
            toll[idx] = parse(Float64, numbers[9])
            link_type[idx] = parse(Float64, numbers[10])

            idx = idx + 1
        end
    end

    ##################################################
    # Trip Table
    ##################################################

    number_of_zones_trip = 0
    total_od_flow = 0

    f = open(trip_table_file, "r")
    fe = readlines(f)

    for line in fe
        #println(line)
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones_trip = parse( Int, line[ search(line, '>')+1 : end ] )
        elseif contains(line, "<TOTAL OD FLOW>")
            total_od_flow = parse( Float64, line[ search(line, '>')+1 : end ] )
        elseif contains(line, "<END OF METADATA>")
            break
        end
    end

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = []
    origin = 0
    for line in fe
        if contains(line, "Origin")
            origin = parse( Int, split(line)[2] )
        elseif contains(line, ";")
            pairs = split(line, ";")
            for i=1:size(pairs)[1]
                if contains(pairs[i], ":")
                    pair = split(pairs[i], ":")
                    destination = parse( Int, strip(pair[1]) )
                    od_flow = parse( Float64, strip(pair[2]) )
                    travel_demand[origin, destination] = od_flow
                    push!(od_pairs, (origin, destination))
                   # println("origin=$origin, destination=$destination, flow=$od_flow")
                end
            end
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        network_name,
        number_of_zones,
        number_of_nodes,
        first_thru_node,
        number_of_links,
        start_node,
        end_node,
        capacity,
        link_length,
        free_flow_time,
        B,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective)

    return ta_data

end # end of load_network

type Arc
    initNode::Int
    termNode::Int
    capacity::Float64
    freeflowtime::Float64
    flow::Float64
end

Arc(initNode::Int, termNode::Int, capacity::Float64, freeflowtime::Float64) =
    Arc(initNode, termNode, capacity, freeflowtime, 0.);

using Graphs

function create_graph(start_node, en_node)
    @assert Base.length(start_node)==Base.length(en_node)

    no_node = max(maximum(start_node), maximum(en_node))
    no_arc = Base.length(start_node)

    graph = simple_inclist(no_node)
    for i=1:no_arc
        add_edge!(graph, start_node[i], en_node[i])
    end
    return graph
end

function get_vector(state, origin, destination, link_dic)
    current = destination
    parent = -1
    x = zeros(Int, maximum(link_dic))

    while parent != origin
        parent = state.parents[current]
        link_idx = link_dic[parent,current]
        if link_idx != 0
            x[link_idx] = 1
        end
        current = parent
    end
    return x
end

function BPR(flowVec, fcoeffs, free_flow_time, capacity)
    bpr = similar(flowVec)
    for a = 1:length(bpr)
        bpr[a] = free_flow_time[a] * sum([fcoeffs[i] * (flowVec[a]/capacity[a])^(i-1) for i = 1:length(fcoeffs)])
    end
    return bpr
end

function all_or_nothing(graph, link_dic, travel_time, demands, start_node, numZones)
    state = []
    path = []
    x = zeros(size(start_node))

    for r=1:numZones
        # for each origin node r, find shortest paths to all destination nodes
        state = dijkstra_shortest_paths(graph, travel_time, r)

        for s=1:numZones
            # for each destination node s, find the shortest-path vector
            # load travel demand
            x = x + demands[(r,s)] * get_vector(state, r, s, link_dic)
        end
    end

    return x
end

function iniDemand(trip_file, numZones, flag=0)
    file = open(trip_file)
    demands = Dict{}()
    for s=1:numZones
        for t=1:numZones
            demands[(s,t)] = 0
        end
    end
    s = 0
    for line in eachline(file)
        if contains(line, "Origin")
            s = Int(parse(Float64,(split(line)[2])))
        elseif contains(line, ";")
            pairs = split(line, ";")
            for pair in pairs[1:end-1]
                if !contains(pair, "\n")
                    pair_vals = split(pair, ":")
                    t, demand = Int(parse(Float64,pair_vals[1])), parse(Float64,pair_vals[2])
                    demands[(s,t)] = demand
                end
            end
        end
    end
    close(file)
    return demands
end

function tapMSA(graph, link_dic, demands, fcoeffs, free_flow_time, capacity, start_node, en_node, numZones, numIter=500, tol=1e-6)
    # Finding a starting feasible solution
    travel_time = BPR(zeros(numLinks), fcoeffs, free_flow_time, capacity)
    xl = all_or_nothing(graph, link_dic, travel_time, demands, start_node, numZones)

    l = 1

    while l < numIter
        l += 1

        xl_old = xl

        # Finding yl
        travel_time = BPR(xl, fcoeffs, free_flow_time, capacity )

        yl = all_or_nothing(graph, link_dic, travel_time, demands, start_node, numZones)

        # assert(yl != xl)

        xl = xl + (yl - xl)/l

        xl_new = xl

        relative_gap = norm(xl_new - xl_old, 1) / norm(xl_new, 1)

        if relative_gap < tol
            break
        end

    end

    tapFlows = Dict()

    for i = 1:length(start_node)
        key = (start_node[i], en_node[i])
        tapFlows[key] = xl[i]
    end

    tapFlowVect = xl

    return tapFlows, tapFlowVect
end

function extract_demandDict(out_dir, files_ID, month_w ,day, instance1, numZones)
    demandsDict = Dict()
    # get ground trueth demands, indexed by 0
    demandsDict[0] = iniDemand(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * string(day) * "_" * instance1 * ".txt", numZones)
    # get initial demands, indexed by 1
    demandsDict[1] = iniDemand(out_dir * "data_traffic_assignment_uni-class/" * files_ID * "_trips_" * month_w * "_" * string(day) * "_" * instance1 * ".txt", numZones, 1)

    return demandsDict
end

function demandsDicToVec(demandsDic, odPairLabel_)
    demandsVec = zeros(length(odPairLabel_))
    for i = 1:length(demandsVec)
        demandsVec[i] = demandsDic[(odPairLabel_["$i"][1], odPairLabel_["$i"][2])]
    end
    return demandsVec
end

function demandsVecToDic(demandsVec, odPairLabel_)
    demandsDic = Dict{}()
    for i = 1:numNodes
        demandsDic[(i, i)] = 0
    end
    for i = 1:length(demandsVec)
        demandsDic[(odPairLabel_["$i"][1], odPairLabel_["$i"][2])] = demandsVec[i]
    end
    return demandsDic
end

function arcData(arc_file)
    arcs = Dict()
    file = open(arc_file)
    inHeader=true
    for line in eachline(file)
        if inHeader
            inHeader = !contains(line, "Init node")
            continue
        end
        vals = split(line, )
        arcs[(parse(Int, vals[1]), parse(Int, vals[2]))] = Arc(parse(Int, vals[1]), parse(Int, vals[2]), parse(Float64, vals[3]), parse(Float64, vals[5]))
    end
    close(file)
    return arcs
end

# add flow data to arcs
function observFlow(arc_file, tapFlowDic)
    arcs = arcData(arc_file)
    ix = 0
    for key in keys(arcs)
        arcs[key].flow = tapFlowDic[key]
    end
    return arcs
end

function objF(graph, link_dic, gamma1, gamma2, demandsVec, demandsVec0, tapFlowVecDict, fcoeffs, free_flow_time, capacity, start_node, en_node, numZones, odPairLabel_)
    demandsDic = demandsVecToDic(demandsVec, odPairLabel_)
    tapFlowVec = tapMSA(graph, link_dic, demandsDic, fcoeffs, free_flow_time, capacity, start_node, en_node, numZones)[2]
    return gamma1 * sum([(demandsVec[i] - demandsVec0[i])^2 for i = 1:length(demandsVec)]) + gamma2 * sum([(tapFlowVec[a] - tapFlowVecDict[a])^2 for a = 1:length(tapFlowVec)])
end

function tapFlowVecToLinkCostDict(tapFlowVec, fcoeffsInvVI, free_flow_time, capacity)
    linkCostVec = BPR(tapFlowVec, fcoeffsInvVI, free_flow_time, capacity)
    temp_dict = Dict{}()
    for i in 1:length(linkCostVec)
        temp_dict["$(i-1)"] = linkCostVec[i]
    end
    return temp_dict
end

# compute the gradient
function gradient_(gamma1, gamma2, demandsVec, demandsVec0, tapFlowVec, observFlowVec, jacob, numODpairs, numLinks)
    gradi = zeros(numODpairs)
    for i = 1:numODpairs
        gradi[i] = 2 * gamma1 * (demandsVec[i] - demandsVec0[i]) + 2 * gamma2 * sum([(tapFlowVec[j] - observFlowVec[j]) * jacob[i, j] for j = 1:numLinks])
    end
    return gradi
end

function descDirec(gamma1, gamma2, demandsVec, demandsVec0, tapFlowVec, observFlowVec, jacob, numODpairs, numLinks)
    gradi = gradient_(gamma1, gamma2, demandsVec, demandsVec0, tapFlowVec, observFlowVec, jacob, numODpairs, numLinks)
    h = similar(gradi)
    for i = 1:length(gradi)
        h[i] = -1 * gradi[i]
    end
    return h
end

# compute a search direction
function searchDirec(demandsVec, descDirect, epsilon_1)
    h = descDirect
    h_ = similar(h)
    for i = 1:length(h)
            if (demandsVec[i] > epsilon_1) || (demandsVec[i] <= epsilon_1 && h[i] > 0)
            h_[i] = h[i]
        else
            h_[i] = 0
        end
    end
    return h_
end

# line search
function thetaMax(demandsVec, searchDirect)
    h_ = searchDirect
    thetaList = Float64[]
    for i = 1:length(h_)
        if h_[i] < 0
            push!(thetaList, - demandsVec[i]/h_[i])
        end
    end
    theta_max = minimum(thetaList)
    return theta_max
end

function armijo(gamma1, gamma2, objFunOld, demandsVecOld, demandsVec0, tapFlowVecDict, fcoeffs, searchDirec, thetaMax, Theta, N, graph, link_dic, free_flow_time, capacity, start_node, en_node, numZones , odPairLabel_)
    demandsVecList = Array{Float64}[]
    objFunList = Float64[]
    push!(demandsVecList, demandsVecOld)
    push!(objFunList, objFunOld)
    for n = 0:N
        demandsVecNew = similar(demandsVecOld)
        for i = 1:length(demandsVecOld)
            demandsVecNew[i] = demandsVecOld[i] + (thetaMax/(Theta^n)) * searchDirec[i]
        end

        objFun_New = objF(graph, link_dic, gamma1, gamma2, demandsVecNew, demandsVec0, tapFlowVecDict, fcoeffs, free_flow_time, capacity, start_node, en_node, numZones , odPairLabel_)
        push!(demandsVecList, demandsVecNew)
        push!(objFunList, objFun_New)
    end
    idx = indmin(objFunList)
    objFunNew = objFunList[idx]
    assert(objFunNew <= objFunOld)
    return demandsVecList[idx], objFunNew
end

function BPRSocial(flowVec, fcoeffs, free_flow_time, capacity)
    bpr = similar(flowVec)
    # refer to [Page 50; Patriksson 1994, 2015]
    for a = 1:length(bpr)
        bpr[a] = free_flow_time[a] * sum([fcoeffs[i] * (flowVec[a]/capacity[a])^(i-1) for i = 1:length(fcoeffs)])
        + free_flow_time[a] * sum([fcoeffs[i] * (i-1) * (flowVec[a]/capacity[a])^(i-1) for i = 2:length(fcoeffs)])
    end
    return bpr
end

function tapMSASocial(demands, fcoeffs, graph, link_dic, start_node, en_node, free_flow_time, capacity, numLinks, numZones, numIter=1000, tol=1e-6)
    # Finding a starting feasible solution
    travel_time = BPRSocial(zeros(numLinks), fcoeffs, free_flow_time, capacity)
    xl = all_or_nothing(graph, link_dic, travel_time, demands, start_node, numZones)

    l = 1

    while l < numIter
        l += 1
        xl_old = xl

        # Finding yl
        travel_time = BPRSocial(xl, fcoeffs, free_flow_time, capacity)
        yl = all_or_nothing(graph, link_dic, travel_time, demands, start_node, numZones)

        # assert(yl != xl)
        xl = xl + (yl - xl)/l
        xl_new = xl

        relative_gap = norm(xl_new - xl_old, 1) / norm(xl_new, 1)
        if relative_gap < tol
            break
        end

    end

    tapFlows = Dict()

    for i = 1:length(start_node)
        key = (start_node[i], en_node[i])
        tapFlows[key] = xl[i]
    end

    tapFlowVect = xl

    return tapFlows, tapFlowVect
end

function socialObj(linkFlowVec, free_flow_time, polyDeg, fcoeffs, capacity, numLinks)
    objVal = sum([sum([free_flow_time[a] * fcoeffs[i] * linkFlowVec[a]^i / capacity[a]^(i-1) for i=1:polyDeg])
        for a = 1:numLinks])
    return objVal
end

using JuMP, Ipopt
function socialOpt(out_dir, files_ID, instance1, demandsVec, polyDeg, free_flow_time, fcoeffs, capacity)

    #load OD pair-route incidence
    odPairRoute = readstring(out_dir * "od_pair_route_incidence_" * instance1 *  files_ID * ".json");
    odPairRoute = JSON.parse(odPairRoute);

    for route in keys(odPairRoute)
        if odPairRoute[route]>0
            odPairRoute[route] = 1
        end
    end

    #load link labels
    linkLabel = readstring(out_dir * "link_label_dict.json");
    linkLabel = JSON.parse(linkLabel);

    #load link-route incidence
    linkRoute = readstring(out_dir * "link_route_incidence_" * instance1 *  files_ID * ".json");
    linkRoute = JSON.parse(linkRoute);

    m = Model(solver=IpoptSolver());
    numLinks = length(linkLabel)
    numRoute = length(linkRoute)
    #numRoute = 400

    numOD = length(demandsVec)

    @variable(m, linkFlow[1:numLinks])
    @variable(m, pathFlow[1:numRoute])

    pathFlowSum = Dict()

    for i=1:numOD
        pathFlowSum[i] = 0
        for j=1:numRoute
            if "$(i)-$(j)" in keys(odPairRoute)
                pathFlowSum[i] += pathFlow[j]
            end
        end
        @constraint(m, pathFlowSum[i] == demandsVec[i])
    end

    pathFlowLinkSum = Dict()

    for a=1:numLinks
        pathFlowLinkSum[a] = 0
        for j=1:numRoute
            if "$(a)-$(j)" in keys(linkRoute)
                pathFlowLinkSum[a] += pathFlow[j];
            end
        end
        @constraint(m, pathFlowLinkSum[a] == linkFlow[a]);
    end

    for j=1:numRoute
        @constraint(m, pathFlow[j] >= 0);
    end

    for a=1:numLinks
        @constraint(m, pathFlowLinkSum[a] <= capacity[a]);
    end

    #@expression(m, f, sum{free_flow_time[a]*linkFlow[a] + .03*free_flow_time[a]*((linkFlow[a])^5)/((capacity[a])^4), a = 1:numLinks} )

#   @NLexpression(m, f, sum{ free_flow_time[a] * fcoeffs[i]  *linkFlow[a]^i / capacity[a]^(i-1) , i = 1:polyDeg , a = 1:numLinks }) ;

    @NLexpression(m, f, sum{free_flow_time[a] * fcoeffs[1] * linkFlow[a] +
            free_flow_time[a] * fcoeffs[2] * linkFlow[a]^2 / capacity[a] +
            free_flow_time[a] * fcoeffs[3] * linkFlow[a]^3 / capacity[a]^2 +
            free_flow_time[a] * fcoeffs[4] * linkFlow[a]^4 / capacity[a]^3 +
            free_flow_time[a] * fcoeffs[5] * linkFlow[a]^5 / capacity[a]^4 +
            free_flow_time[a] * fcoeffs[6] * linkFlow[a]^6 / capacity[a]^5 +
            free_flow_time[a] * fcoeffs[7] * linkFlow[a]^7 / capacity[a]^6 +#, a = 1:numLinks})
            free_flow_time[a] * fcoeffs[8] * linkFlow[a]^8 / capacity[a]^7, a = 1:numLinks})

    @NLobjective(m, Min, f);
    #print(m)

    solve(m);

    return getobjectivevalue(m)
end

objFunDict = Dict()

function adjustingODdemands(instance, day, cnt, out_dir, month_w, month, year, files_ID, ta_data, link_list,numNodes, gamma1, gamma2 )
    ## Load initialization
    flow_observ, y, link_vector = get_observed_flow.get_observed_flow(out_dir, files_ID, instance, month_w,  month, year, week_day_Apr_list);

    # demands from GLS
    g_0 = read_demand_file(out_dir * "OD_demands/OD_demand_matrix_" * month_w * "_" * string(day) * "_weekday_" * instance * files_ID * ".txt", numNodes)

    # cost function parameters
    coeffs_dict_Apr_PM_ = readstring(out_dir * "coeffs_dict_" * month_w *  "_" * instance * ".json");
    coeffs_dict_Apr_PM_ = JSON.parse(coeffs_dict_Apr_PM_);

    # cross-validation best key selection
    best_key = readstring(out_dir * "cross_validation_best_key/cross_validation_best_key_" * month_w * "_" * string(day) * "_" * instance * ".json");
    best_key = JSON.parse(best_key);
    best_key = "(8, 0.5, 1.0, 2)"
    fcoeffs = coeffs_dict_Apr_PM_[best_key];
    #fcoeffs = [1,0,0,0,.5,0,0,0]
    #fcoeffs = [1.0,-0.00302509, 0.0577279, -0.195632, 0.620696, -0.905963, 0.936143, -0.469483, 0.108584];
    polyDeg = length(fcoeffs);

    # load network
    ta_data = load_ta_network_(out_dir, files_ID, month_w, day,  instance);

    # Load OdPairLabel
    odPairLabel_ = readstring(out_dir * "od_pair_label_dict__refined.json");
    odPairLabel_ = JSON.parse(odPairLabel_);

    # Renaming variables
    numNodes = maximum(map(pair->pair[1], keys(g_0)));
    start_node = ta_data.start_node;
    end_node = ta_data.end_node;
    capacity = ta_data.capacity;
    free_flow_time = ta_data.free_flow_time;
    number_of_zones = ta_data.number_of_zones;
    numLinks = size(start_node)[1];

    numODpairs = numNodes * (numNodes - 1);
    graph = create_graph(start_node, end_node);
    link_list = sparse(start_node, end_node, 1:numLinks);

    # Run MSA and adjust demand vector
    tapFlows = Dict();
    tapFlowVect = Dict();
    tapFlowDicDict = Dict();
    tapFlowVecDict = Dict();
    linkCostDicDict = Dict();
    jacobiSpiessDict = Dict();
    jacobDict = Dict();
    descDirecDict = Dict();
    searchDirecDict = Dict();
    thetaMaxDict = Dict();
    demandsVecDict = Dict();
    demandsDiffDict = Dict();
    norObjFunDict = Dict();
    tapSocialFlowDicDict = Dict();
    tapSocialFlowVecDict = Dict();
    tapSocialFlowDicDict_ = Dict();
    tapSocialFlowVecDict_ = Dict();
    user_sol_dict = Dict();
    social_sol_dict = Dict();
    PoA_dict = Dict();

    #demandsVecDict[0] = demandsDicToVec(g_0, odPairLabel_);
    #demandsDiffDict[1] = norm(demandsDicToVec(demandsDict[1], odPairLabel_) - demandsDicToVec(g_0, odPairLabel_))/
    #                     norm(demandsDicToVec(g_0, odPairLabel_));

    if isdir(out_dir * "demandsDict") == false
            mkdir(out_dir * "demandsDict");
    end

    #obj_dict[day] = demandsDictFixed(demandsDict, flow_observ,link_vector, graph, ta_data, link_dic, day, gamma1, gamma2,
    #    out_dir, files_ID, month_w, instance, key_, free_flow_time, capacity, start_node, en_node, numZones, cnt)

    xl = flow_observ[:, cnt]
    tapFlows[0] = Dict()
    for i = 1:length(start_node)
        key = link_vector[i]
        tapFlows[0][key] = xl[i]
    end

    tapFlowVecDict[0] = [];
    for i = 1:length(start_node)
        key = link_vector[i];
        append!(tapFlowVecDict[0], tapFlows[0][key])
    end

    network_data_file = files_ID * "_net_" * month_w * "_full_" * instance * ".txt";
    arcsDict = observFlow(out_dir  * "data_traffic_assignment_uni-class/" * network_data_file, tapFlows[0]);

    demandsVecDict[0] = demandsDicToVec(g_0, odPairLabel_);
    #demandsVecDict[0] = ones(length(demandsVecDict[0]))
    demandsVecDict[1] = demandsDicToVec(g_0, odPairLabel_);
    #demandsVecDict[1] = ones(length(demandsVecDict[1]))

    objFunDict[1] = objF(graph, link_list, gamma1, gamma2, demandsVecDict[0],
        demandsVecDict[0], tapFlowVecDict[0], fcoeffs, free_flow_time, capacity, start_node, end_node, numZones, odPairLabel_);

    demandsDic = demandsVecToDic(demandsVecDict[0], odPairLabel_)

    # get initial flow vector (corresponding to initial demands)
    tapFlows[1], tapFlowVecDict[1] = tapMSA(graph, link_list, g_0, fcoeffs, free_flow_time, capacity, start_node, end_node, number_of_zones);

    # Computation of a Descent direction
    linkCostDicDict[1] = tapFlowVecToLinkCostDict(tapFlowVecDict[1], fcoeffs, free_flow_time, capacity);
    #linkCostDicDict[1]["0"], link_length_list[1]
    jacobiSpiessDict[1] = Compute_Jacobian.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs, link_list_js, [linkCostDicDict[1]["$(i)"] for i=0:numLinks-1]);

    demandsDict = extract_demandDict(out_dir, files_ID, month_w ,day, instance, numZones);

    # maximum number of iterations
    N = 100;

    # Armijo rule parameters
    rho = 2;
    M = 10;

    # search direction parameter
    epsilon_1 = 0;

    # stop criterion parameter
    epsilon_2 = 1e-20;
    obj = 9999999999999

    las = 0

    for l = 1:N
        jacobDict[l] = jacobiSpiessDict[l]

        descDirecDict[l] = descDirec(gamma1, gamma2, demandsVecDict[l], demandsVecDict[0],  tapFlowVecDict[l],
                tapFlowVecDict[0], jacobDict[l], numODpairs, numLinks);

        demandsVecDict[l] = demandsDicToVec(demandsDict[l], odPairLabel_);

        searchDirecDict[l] = searchDirec(demandsVecDict[l], descDirecDict[l], epsilon_1);

        thetaMaxDict[l] = thetaMax(demandsVecDict[l], searchDirecDict[l]);

        demandsVecDict[l+1] = similar(demandsVecDict[0]);

        demandsVecDict[l+1], objFunDict[l+1] = armijo(gamma1, gamma2, objFunDict[l], demandsVecDict[l], demandsVecDict[0], tapFlowVecDict[0], fcoeffs,
                searchDirecDict[l], thetaMaxDict[l], rho, M, graph, link_list, free_flow_time, capacity, start_node, end_node, numZones, odPairLabel_);

        demandsDict[l+1] = demandsVecToDic(demandsVecDict[l+1], odPairLabel_);

        tapFlows[l+1], tapFlowVecDict[l+1] = tapMSA(graph, link_list, demandsDict[l+1], fcoeffs, free_flow_time, capacity, start_node, end_node, numZones);

        arcsDict[l+1] = observFlow(out_dir  * "data_traffic_assignment_uni-class/" * network_data_file, tapFlows[l+1]);

        linkCostDicDict[l+1] = tapFlowVecToLinkCostDict(tapFlowVecDict[l+1], fcoeffs, free_flow_time, capacity);

        jacobiSpiessDict[l+1] = Compute_Jacobian.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                                      link_list_js, [linkCostDicDict[l+1]["$(i)"] for i=0:numLinks-1]);

        demandsDiffDict[l+1] = norm(demandsVecDict[l+1] - demandsVecDict[0]) / norm(demandsVecDict[0]);

        obj =  objFunDict[1] - objFunDict[l]

        las = l
            # stopping criterion
        if (objFunDict[l] - objFunDict[l+1]) / objFunDict[1] < epsilon_2
            break
        end

        println("iteration $(l) finished...")

    end

    # normalize objective function value
    for a = 1:(length(objFunDict))
        norObjFunDict[a] = objFunDict[a] / objFunDict[1];
    end

    # write files
    outfile = open(out_dir * "demandsDict/demandsVecDict$(day)_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, demandsVecDict)
    close(outfile)

    outfile = open(out_dir * "demandsDict/demandsDict$(day)_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, demandsDict)
    close(outfile)

    outfile = open(out_dir * "demandsDict/tapFlowDicDict$(day)_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, tapFlows)
    close(outfile)

    outfile = open(out_dir * "demandsDict/tapFlowVecDict$(day)_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, tapFlowVecDict)
    close(outfile)

    outfile = open(out_dir * "demandsDict/jacobi$(day)_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, jacobiSpiessDict[las])
    close(outfile)

    demandsDict[length(demandsDict)-1]
    demandsDict_ = Dict{}()

    for key in keys(demandsDict[length(demandsDict)-1])
        demandsDict_[key] = demandsDict[length(demandsDict)-1][key]
    end

    outfile = open(out_dir * "demandsDict/demandsDictFixed$(day)_"* month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, demandsDict_)
    close(outfile)

    # Calculating POA
    tapSocialFlowDicDict_[day], tapSocialFlowVecDict_[day] =  tapMSASocial(demandsDict[las], fcoeffs, graph, link_list, start_node, end_node, free_flow_time, capacity, numLinks, numZones);
    #tapSocialFlowDicDict[day] = socialOpt(out_dir, files_ID, instance, demandsVecDict[las], polyDeg, free_flow_time, fcoeffs, capacity)
    #print(tapSocialFlowDicDict[day])

    #user_sol_dict[day] = socialObj(flow_observ[:, cnt], free_flow_time, polyDeg, fcoeffs, capacity, numLinks) ;
    user_sol_dict[day] = socialObj(tapFlowVecDict[las], free_flow_time, polyDeg, fcoeffs, capacity, numLinks) ;
    tapSocialFlowDicDict[day] = socialObj(tapSocialFlowVecDict_[day], free_flow_time, polyDeg, fcoeffs, capacity, numLinks);

    PoA_dict[day] = user_sol_dict[day] / tapSocialFlowDicDict[day];

    println("-----------------------------------------")
    println("day $(day) finished...")
    println("PoA for day $(day) and instance " * instance * " is " * string(PoA_dict[day]))
    println("-----------------------------------------")

    return PoA_dict[day], user_sol_dict[day], obj
end

PoA_dict = Dict()
day = 3
cnt = 2
instance = "PM"
gamma1 = 0
gamma2 = 1
ta_data = load_ta_network_(out_dir, files_ID, month_w, day,  instance);
link_list = sparse(ta_data.start_node, ta_data.end_node, 1:numLinks);
instances_1 = ["PM"]
cnt = 0
for instance in instances_1
    PoA_dict = Dict()
    cong = Dict()
    obj = Dict()
    for day in week_day_Apr_list
        cnt += 1
        PoA_dict[day], cong[day], obj[day] = adjustingODdemands(instance, day, cnt, out_dir, month_w, month, year, files_ID, ta_data, link_list, numNodes, gamma1, gamma2 )
    end

    outfile =  open(out_dir * "cong_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, cong)
    close(outfile)

    outfile =  open(out_dir * "PoA_dict_" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, PoA_dict)
    close(outfile)

    outfile =  open(out_dir * "obj_dict" * month_w * "_" * instance * ".json", "w")
    JSON.print(outfile, obj)
    close(outfile)
end


text:
# Magic numbers
@constants MAGICS "" begin
    const MH_MAGIC = 0xfeedface
    const MH_CIGAM = bswap(MH_MAGIC)

    const MH_MAGIC_64 = 0xfeedfacf
    const MH_CIGAM_64 = bswap(MH_MAGIC_64)

    const FAT_MAGIC = 0xCAFEBABE
    const FAT_CIGAM = bswap(FAT_MAGIC)
end

const CPU_ARCH_MASK  = 0xff000000
const CPU_ARCH_ABI64 = 0x01000000

@constants CPUTYPES "CPU_TYPE_" begin
    const CPU_TYPE_ANY          = reinterpret(UInt32,Int32(-1))
    const CPU_TYPE_VAX          = 1
    # skip
    # skip
    # skip
    # skip
    const CPU_TYPE_MC680x0      = 6
    const CPU_TYPE_X86          = 7
    const CPU_TYPE_X86_64       = CPU_TYPE_X86 | CPU_ARCH_ABI64
    const CPU_TYPE_MIPS         = 8
    # skip
    const CPU_TYPE_MC98000      = 10
    const CPU_TYPE_HPPA         = 11
    const CPU_TYPE_ARM          = 12
    const CPU_TYPE_MC88000      = 13
    const CPU_TYPE_SPARC        = 14
    const CPU_TYPE_I860         = 15
    const CPU_TYPE_ALPHA        = 16
    # skip
    const CPU_TYPE_POWERPC      = 18
    const CPU_TYPE_POWERPC64    = CPU_TYPE_POWERPC | CPU_ARCH_ABI64
end

# TODO subtype constants
@constants NLISTTYPES "N_" begin
    const N_STAB         = 0xe0
    const N_PEXT         = 0x10
    const N_TYPE         = 0x0e
    const N_EXT          = 0x01
    const N_UNDF         = 0x00
    const N_ABS          = 0x02
    const N_SECT         = 0x0e
    const N_PBUD         = 0x0c
    const N_INDR         = 0x0a
end

@constants NDESC "N_" begin
    const N_WEAK_REF      = 0x0040
    const N_WEAK_DEF      = 0x0080
    const N_ARM_THUMB_DEF = 0x0008
end

const N_REF_TO_WEAK = N_WEAK_DEF

const NO_SECT  = 0x0
const MAX_SECT = 0xff

@constants FILETYPES "MH_" begin
    const MH_OBJECT      = 0x1     # relocatable object file
    const MH_EXECUTE     = 0x2     # demand paged executable file
    const MH_FVMLIB      = 0x3     # fixed VM shared library file
    const MH_CORE        = 0x4     # core file
    const MH_PRELOAD     = 0x5     # preloaded executable file
    const MH_DYLIB       = 0x6     # dynamically bound shared library
    const MH_DYLINKER    = 0x7     # dynamic link editor
    const MH_BUNDLE      = 0x8     # dynamically bound bundle file
    const MH_DYLIB_STUB  = 0x9     # shared library stub for static linking only
    const MH_DSYM        = 0xa     # companion file with only debug sections
    const MH_KEXT_BUNDLE = 0xb     # x86_64 kexts
end

# Constants for the flags field of the mach_header

@constants FLAGS "MH_" begin
    # the object file has no undefined references
    const MH_NOUNDEFS                   = 0x1
    # the object file is the output of an incremental link
    # against a base file and can't be link edited again
    const MH_INCRLINK                   = 0x2
    # the object file is input for the dynamic linker
    # and can't be staticly link edited again
    const MH_DYLDLINK                   = 0x4
    # the object file's undefined references are bound
    # by the dynamic linker when loaded
    const MH_BINDATLOAD                 = 0x8
    # the file has its dynamic undefined references prebound
    const MH_PREBOUND                   = 0x10
    # the file has its read-only and read-write segments split
    const MH_SPLIT_SEGS                 = 0x20
    # the shared library init routine is to be run
    # lazily via catching memory faults to its writeable segments
    # (obsolete)
    const MH_LAZY_INIT                  = 0x40
    # the image is using two-level name space bindings
    const MH_TWOLEVEL                   = 0x80
    # the executable is forcing all image to use flat
    # name space bindings
    const MH_FORCE_FLAT                 = 0x100
    # this umbrella guarantees no multiple defintions of
    # symbols in its sub-images so the two-level namespace
    # hints can always be used
    const MH_NOMULTIDEFS                = 0x100
    # do not have dyld notify the prebinding agent about this
    # executable
    const MH_NOFIXPREBINDING            = 0x200
    # the binary is not prebound but can have its prebinding
    # redone. only used when MH_PREBOUND is not set.
    const MH_PREBINDABLE                = 0x400
    # indicates that this binary binds to all two-level
    # namespace modules of its dependent libraries.
    # only used when MH_PREBINDABLE and MH_TWOLEVEL are both set
    const MH_ALLMODSBOUND               = 0x800
    # safe to divide up the sections into sub-sections
    # via symbols for dead code stripping
    const MH_SUBSECTIONS_VIA_SYMBOLS    = 0x1000
    const MH_CANONICAL                  = 0x2000
    const MH_WEAK_DEFINES               = 0x4000
    const MH_BINDS_TO_WEAK              = 0x8000
    const MH_ALLOW_STACK_EXECUTION      = 0x10000
    const MH_ROOT_SAFE                  = 0x20000
    const MH_SETUID_SAFE                = 0x40000
    const MH_NO_REEXPORTED_DYLIBS       = 0x80000
    const MH_PIE                        = 0x100000
    const MH_DEAD_STRIPPABLE_DYLIB      = 0x200000
    const MH_HAS_TLV_DESCRIPTORS        = 0x400000
    const MH_NO_HEAP_EXECUTION          = 0x800000
end

# Load command types.  Note that we have dropped the LC_REQ_DYLD flag OR'ed
# into many of these types to aid in easy usage of these constants.  Since
# the meaning of LC_REQ_DYLD is to denote to the linker which load commands
# cannot be safely ignored when loading, this does not affect us much, and
# we disregard it entirely.

const LC_REQ_DYLD = 0x80000000
@constants LCTYPES "LC_" begin
    const LC_SEGMENT            = 0x1   # segment of this file to be mapped
    const LC_SYMTAB             = 0x2   # link-edit stab symbol table info
    const LC_SYMSEG             = 0x3   # link-edit gdb symbol table info (obsolete)
    const LC_THREAD             = 0x4   # thread
    const LC_UNIXTHREAD         = 0x5   # unix thread (includes a stack)
    const LC_LOADFVMLIB         = 0x6   # load a specified fixed VM shared library
    const LC_IDFVMLIB           = 0x7   # fixed VM shared library identification
    const LC_IDENT              = 0x8   # object identification info (obsolete)
    const LC_FVMFILE            = 0x9   # fixed VM file inclusion (internal use)
    const LC_PREPAGE            = 0xa   # prepage command (internal use)
    const LC_DYSYMTAB           = 0xb   # dynamic link-edit symbol table info
    const LC_LOAD_DYLIB         = 0xc   # load a dynamically linked shared library
    const LC_ID_DYLIB           = 0xd   # dynamically linked shared lib ident
    const LC_LOAD_DYLINKER      = 0xe   # load a dynamic linker
    const LC_ID_DYLINKER        = 0xf   # dynamic linker identification
    const LC_PREBOUND_DYLIB     = 0x10  # modules prebound for a dynamically linked shared library
    const LC_ROUTINES           = 0x11  # image routines
    const LC_SUB_FRAMEWORK      = 0x12  # sub framework
    const LC_SUB_UMBRELLA       = 0x13  # sub umbrella
    const LC_SUB_CLIENT         = 0x14  # sub client
    const LC_SUB_LIBRARY        = 0x15  # sub library
    const LC_TWOLEVEL_HINTS     = 0x16  # two-level namespace lookup hints
    const LC_PREBIND_CKSUM      = 0x17  # prebind checksum
    #
    # load a dynamically linked shared library that is allowed to be missing
    # (all symbols are weak imported).
    #
    const LC_LOAD_WEAK_DYLIB        = 0x18
    const LC_SEGMENT_64             = 0x19 # 64-bit segment of this file to be mapped
    const LC_ROUTINES_64            = 0x1a # 64-bit image routines
    const LC_UUID                   = 0x1b # the uuid
    const LC_RPATH                  = 0x1c # runpath additions
    const LC_CODE_SIGNATURE         = 0x1d # local of code signature
    const LC_SEGMENT_SPLIT_INFO     = 0x1e # local of info to split segments
    const LC_REEXPORT_DYLIB         = 0x1f # load and re-export dylib
    const LC_LAZY_LOAD_DYLIB        = 0x20 # delay load of dylib until first use
    const LC_ENCRYPTION_INFO        = 0x21 # encrypted segment information
    const LC_DYLD_INFO              = 0x22 # compressed dyld information
    const LC_DYLD_INFO_ONLY         = 0x22 # compressed dyld information only
    const LC_LOAD_UPWARD_DYLIB      = 0x23 # load upward dylib
    const LC_VERSION_MIN_MACOSX     = 0x24 # build for MacOSX min OS version
    const LC_VERSION_MIN_IPHONEOS   = 0x25 # build for iPhoneOS min OS version
    const LC_FUNCTION_STARTS        = 0x26 # compressed table of function start addresses
    const LC_DYLD_ENVIRONMENT       = 0x27 # string for dyld to treat like environment variable
    const LC_MAIN                   = 0x28 # replacement for LC_UNIXTHREAD
    const LC_DATA_IN_CODE           = 0x29 # table of non-instructions in __text
    const LC_SOURCE_VERSION         = 0x2A # source version used to build binary
    const LC_DYLIB_CODE_SIGN_DRS    = 0x2B # Code signing DRs copied from linked dylibs
    const LC_ENCRYPTION_INFO_64     = 0x2C # 64-bit encrypted segment information
    const LC_LINKER_OPTION          = 0x2D # linker options in MH_OBJECT files

end

const SECTION_TYPE          = 0x000000ff # 256 section types
const SECTION_ATTRIBUTES    = 0xffffff00 # 24 section attributes

@constants SECTYPES "" begin
    # regular section
    const S_REGULAR                                 = 0x0
    # zero fill on demand section
    const S_ZEROFILL                                = 0x1
    # section with only literal C string
    const S_CSTRING_LITERALS                        = 0x2
    # section with only 4 byte literals
    const S_4BYTE_LITERALS                          = 0x3
    # section with only 8 byte literals
    const S_8BYTE_LITERALS                          = 0x4
    # section with only pointers to literals
    const S_LITERAL_POINTERS                        = 0x5
#
#  For the two types of symbol pointers sections and the symbol stubs section
#  they have indirect symbol table entries.  For each of the entries in the
#  section the indirect symbol table entries, in corresponding order in the
#  indirect symbol table, start at the index stored in the reserved1 field
#  of the section structure.  Since the indirect symbol table entries
#  correspond to the entries in the section the number of indirect symbol table
#  entries is inferred from the size of the section divided by the size of the
#  entries in the section.  For symbol pointers sections the size of the entries
#  in the section is 4 bytes and for symbol stubs sections the byte size of the
#  stubs is stored in the reserved2 field of the section structure.

    # section with only non-lazy symbol pointers
    const S_NON_LAZY_SYMBOL_POINTERS                = 0x6
    # section with only lazy symbol
    const S_LAZY_SYMBOL_POINTERS                    = 0x7
    # section with only symbol stubs, byte
    # size of stub in the reserved2 field
    const S_SYMBOL_STUBS                            = 0x8
    # section with only function pointers
    # for initialization
    const S_MOD_INIT_FUNC_POINTERS                  = 0x9
    # section with only function pointers
    # for initialization
    const S_MOD_TERM_FUNC_POINTERS                  = 0xa
    # section contains symbols that
    # are to be coalesced
    const S_COALESCED                               = 0xb
    # zero fill on demand section
    # (that can be larger than 4 gigabytes)
    const S_GB_ZEROFILL                             = 0xc
    # section with only pairs of function
    # pointers for interposing
    const S_INTERPOSING                             = 0xd
    # section with only 16 byte literals
    const S_16BYTE_LITERALS                         = 0xe
    # section contains DTrace Object Format
    const S_DTRACE_DOF                              = 0xf
    # section with only lazy symbol pointers
    # to lazy loaded dylibs
    const S_LAZY_DYLIB_SYMBOL_POINTERS              = 0x10
    # template of initial values for TLVs
    const S_THREAD_LOCAL_REGULAR                    = 0x11
    # template of initial values for TLVs
    const S_THREAD_LOCAL_ZEROFILL                   = 0x12
    # TLV descriptors
    const S_THREAD_LOCAL_VARIABLES                  = 0x13
    # pointers to TLV descriptors
    const S_THREAD_LOCAL_VARIABLE_POINTERS          = 0x14
    # functions to call to initialize TLV values
    const S_THREAD_LOCAL_INIT_FUNCTION_POINTERS     = 0x15
end

const SECTION_ATTRIBUTES_USR = 0xff000000
const SECTION_ATTRIBUTES_SYS = 0x00ffff00

@constants SECATTRS "S_ATTR_" begin
    # section contains only true machine instructions
    const S_ATTR_PURE_INSTRUCTIONS      = 0x80000000
    # section contains coalesced symbols that are not to be
    # in a ranlib table of contents
    const S_ATTR_NO_TOC                 = 0x40000000
    # ok to strip static symbols in this section in files
    # with the MH_DYLDLINK flag
    const S_ATTR_STRIP_STATIC_SYMS      = 0x20000000
    # no dead stripping
    const S_ATTR_NO_DEAD_STRIP          = 0x10000000
    # blocks are live if they reference live blocks
    const S_ATTR_LIVE_SUPPORT           = 0x08000000
    # Used with i386 code stubs written on by dyld
    const S_ATTR_SELF_MODIFYING_CODE    = 0x04000000
    # a debug section
    const S_ATTR_DEBUG                  = 0x02000000
    # section contains some machine instructions
    const S_ATTR_SOME_INSTRUCTIONS      = 0x00000400
    # section has external relocation entries
    const S_ATTR_EXT_RELOC              = 0x00000200
    # section has local relocation entries
    const S_ATTR_LOC_RELOC              = 0x00000100
end

@constants UNWIND_FLAGS "UNWIND_" begin
    const UNWIND_IS_NOT_FUNCTION_START  = 0x80000000
    const UNWIND_HAS_LSDA               = 0x40000000
    const UNWIND_PERSONALITY_MASK       = 0x30000000
end

@constants UNWIND_X86 "UNWIND_X86_" begin
    const UNWIND_X86_MODE_MASK                          = 0x0F000000
    const UNWIND_X86_MODE_EBP_FRAME                     = 0x01000000
    const UNWIND_X86_MODE_STACK_IMMD                    = 0x02000000
    const UNWIND_X86_MODE_STACK_IND                     = 0x03000000
    const UNWIND_X86_MODE_DWARF                         = 0x04000000

    const UNWIND_X86_EBP_FRAME_REGISTERS                = 0x00007FFF
    const UNWIND_X86_EBP_FRAME_OFFSET                   = 0x00FF0000

    const UNWIND_X86_FRAMELESS_STACK_SIZE               = 0x00FF0000
    const UNWIND_X86_FRAMELESS_STACK_ADJUST             = 0x0000E000
    const UNWIND_X86_FRAMELESS_STACK_REG_COUNT          = 0x00001C00
    const UNWIND_X86_FRAMELESS_STACK_REG_PERMUTATION    = 0x000003FF

    const UNWIND_X86_DWARF_SECTION_OFFSET               = 0x00FFFFFF
end

@constants UNWIND_X86_REGS "UNWIND_X86_REG_" begin
    const UNWIND_X86_REG_NONE   = 0
    const UNWIND_X86_REG_EBX    = 1
    const UNWIND_X86_REG_ECX    = 2
    const UNWIND_X86_REG_EDX    = 3
    const UNWIND_X86_REG_EDI    = 4
    const UNWIND_X86_REG_ESI    = 5
    const UNWIND_X86_REG_EBP    = 6
end

@constants UNWIND_X86_64 "UNWIND_X86_64_" begin
    const UNWIND_X86_64_MODE_MASK                          = 0x0F000000
    const UNWIND_X86_64_MODE_EBP_FRAME                     = 0x01000000
    const UNWIND_X86_64_MODE_STACK_IMMD                    = 0x02000000
    const UNWIND_X86_64_MODE_STACK_IND                     = 0x03000000
    const UNWIND_X86_64_MODE_DWARF                         = 0x04000000

    const UNWIND_X86_64_EBP_FRAME_REGISTERS                = 0x00007FFF
    const UNWIND_X86_64_EBP_FRAME_OFFSET                   = 0x00FF0000

    const UNWIND_X86_64_FRAMELESS_STACK_SIZE               = 0x00FF0000
    const UNWIND_X86_64_FRAMELESS_STACK_ADJUST             = 0x0000E000
    const UNWIND_X86_64_FRAMELESS_STACK_REG_COUNT          = 0x00001C00
    const UNWIND_X86_64_FRAMELESS_STACK_REG_PERMUTATION    = 0x000003FF

    const UNWIND_X86_64_DWARF_SECTION_OFFSET               = 0x00FFFFFF
end

@constants UNWIND_X86_64_REGS "UNWIND_X86_64_REG_" begin
    const UNWIND_X86_64_REG_NONE    = 0
    const UNWIND_X86_64_REG_RBX     = 1
    const UNWIND_X86_64_REG_R12     = 2
    const UNWIND_X86_64_REG_R13     = 3
    const UNWIND_X86_64_REG_R14     = 4
    const UNWIND_X86_64_REG_R15     = 5
    const UNWIND_X86_64_REG_RBP     = 6
end

@constants UNWIND_SECOND_LEVEL "UNWIND_SECOND_LEVEL_" begin
    const UNWIND_SECOND_LEVEL_REGULAR       = 2
    const UNWIND_SECOND_LEVEL_COMPRESSED    = 3
end

@constants X86_64_RELOC "X86_64_RELOC_" begin
    const X86_64_RELOC_UNSIGNED        = 0
    const X86_64_RELOC_SIGNED          = 1
    const X86_64_RELOC_BRANCH          = 2
    const X86_64_RELOC_GOT_LOAD        = 3
    const X86_64_RELOC_GOT             = 4
    const X86_64_RELOC_SUBTRACTOR      = 5
    const X86_64_RELOC_SIGNED_1        = 6
    const X86_64_RELOC_SIGNED_2        = 7
    const X86_64_RELOC_SIGNED_4        = 8
    const X86_64_RELOC_TLV             = 9
end

@constants X86_THREAD_FLAVORS "" begin
    const x86_THREAD_STATE32    = 1
    const x86_FLOAT_STATE32     = 2
    const x86_EXCEPTION_STATE32 = 3
    const x86_THREAD_STATE64    = 4
    const x86_FLOAT_STATE64     = 5
    const x86_EXCEPTION_STATE64 = 6
    const x86_THREAD_STATE      = 7
    const x86_FLOAT_STATE       = 8
    const x86_EXCEPTION_STATE   = 9
    const x86_DEBUG_STATE32     = 10
    const x86_DEBUG_STATE64     = 11
    const x86_DEBUG_STATE       = 12
end


text:
# ==================================
# Motor
#
# system type: continuous LTI system
# state dimension: 8
# input dimension: 2
# ==================================
using ReachabilityBenchmarks, MathematicalSystems, LazySets, SparseArrays

function motor_model()
    # system matrix
    I = [1, 2, 2, 3, 3, 3, 3, 4, 5, 6, 6, 7, 7, 7, 7, 8]
    J = [2, 3, 2, 1, 2, 3, 4, 1, 6, 7, 6, 5, 6, 7, 8, 5]
    V = [1, 8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0, 1.0,
         8487.2, -1.0865, -2592.1, -21.119, -698.91, -141399.0, 1.0]
    A = sparse(I, J, V)

    # input matrix
    B = sparse([4, 8], [1, 2], [-1.0, -1.0])

    # state domain
    X = Universe(8)

    # input domain
    U = Hyperrectangle([0.23, 0.3], [0.07, 0.1])

    # continuous LTI system
    S = @system(x' = Ax + Bu, x ∈ X, u ∈ U)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
using Plots
using SortFilters
using Random

Random.seed!(5)
len = 200
x = trunc.(Int, 100*randn(len))
x = cumsum(x)
w = 20

# when p represents probability rational number in the range of [0..1]
probability = 0.5
# output is converted to rational type of p
moving_median = movsort(x, w, probability)

# when p represents integer percent number in the range of [0..100]
percent = 50
# output preserves input type and is calculated as simplified quantiles
# from input elements at nearest index
moving_median_int = movsort(x, w, percent)

# also you can get several quantiles at once:
p = (0.25, 0.75)
moving_quartiles = movsort(x, w, p)
# and then calculate interquantile range
moving_iqr = map( x->x[2]-x[1], moving_quartiles)
moving_q25 = map( x->x[1], moving_quartiles)
moving_q75 = map( x->x[2], moving_quartiles)

# notice that quantiles are calculated with a delay of window/2
# start boundary condition is set as a repeated first point
plot(x, label = "x", legend = :topleft)
plot!(moving_median, label = "exact moving median")
plot!(moving_median_int, label = "simplified moving median")
plot!(moving_q25, label = "25th quartile")
plot!(moving_q75, label = "75th quartile")

# png("plot1.png")


text:
## Scratch pad

using Polygons

P = [ .8 .2; .2 .8 ]
U = [ .5 -.5 ; -.5 .5 ]
bet = 0.9
vlow = ( 1 - bet ) * ( ( eye(2) - bet * P ) \ U )

N=40
dirs = zeros( N, 2 )
dirs[:,1] = [ cos(2 * pi * (i-1)/N) for i in 1:N ]
dirs[:,2] = [ sin(2 * pi * (i-1)/N) for i in 1:N ]

# Search directions
U1(s) = [ 2 -1; 3+s s ]
U2(s) = [ 2 3-s; -1 -s ]
# The period payoff functions

function T_unc( W::Array{ Polygon, 1 }, s::Float64, a::Vector, outer = true)
    # s is the value for the state and a is a vector containing the action indices
    u = vec( [ U1(s)[a[1],a[2]] ; U2(s)[a[1],a[2]] ] )
        # Period payoff
    s_idx = ( s == .5 ) ? 1 : 2
        # The index of s.  Required for selecting the right transition probability row
    return ( ( 1 - bet ) * u ) + ( bet * ( weightedSum( W, vec( P[ s_idx, :] ), dirs, outer ) ) )
end

function T_a( W::Array{ Polygon, 1 }, s::Float64, a::Vector, outer = true)
    # calcultes the incentive-compatible value set conditional on the action a
    unc = T_unc( W, s, a, outer )
        # The unconstrained set
    pddev = [ maximum( U1(s)[ :, a[2] ] ), maximum( U2(s)[ a[1], : ] ) ]
        # The deviating period payoffs
    s_idx = ( s == .5 ) ? 1 : 2
        # The index of s.  Required for selecting the right transition probability row
    dev = ( 1 - bet ) * pddev + bet * vlow[ s_idx, : ]'
        # The deviating payoff
    ic = crop( crop( unc, 1, dev[1] ), 2, dev[2] )
        # The incentive compatible set of payoffs\
    return ic
end

A_idx = [ 1 1; 1 2 ; 2 1; 2 2 ]
    # The matrix of action combinations

function T_operator( W::Array{ Polygon, 1 }, states::Vector, outer = true)
    # Forms the T operator
    out = [ union( [ T_a( W, s, vec(A_idx[j,:]))::Polygon for j in 1:size(A_idx)[1] ] )::Polygon for s in states ]
    return out
end

# polyPlot( [ winit1, T_operator( winit, [ .5 , -.5 ] ) ] )

winit1 = grahamScan( [ 2 2; -1 3.5; 3.5 -1; -0.5 0.5 ; 0.5 -0.5 ] )
winit = [ winit1, winit1 ]
aa = [ T_a( winit, -.5, vec(A_idx[j,:]))::Polygon for j in 1:size(A_idx)[1] ]
ff = union( aa )
polyPlot( [ aa, ff ])
polyPlot( aa[2] )

s = -.5
a = vec(A_idx[2,:])
nn=T_a( winit, s, a)::Polygon
unc = T_unc( winit, s, a)
unc2 = grahamScan(unc)
pddev = [ maximum( U1(s)[ :, a[2] ] ), maximum( U2(s)[ a[1], : ] ) ]
    # The deviating period payoffs
s_idx = ( s == .5 ) ? 1 : 2
    # The index of s.  Required for selecting the right transition probability row
dev = ( 1 - bet ) * pddev + bet * vlow[ s_idx, : ]'
    # The deviating payoff
ic = crop( crop( unc, 1, dev[1] ), 2, dev[2] )

polyPlot([unc, nn])


text:
function _parse_to_vars(
    model::Optimizer,
    f::MOI.VectorOfVariables,
)
    return [_info(model, v) for v in f.variables]
end

function _parse_to_coeffs_vars(
    model::Optimizer,
    f::MOI.ScalarAffineFunction{T},
) where {T <: Real}
    f = MOI.Utilities.canonical(f)
    coeffs, vars = _parse_to_coeffs_vars(model, f.terms)
    return coeffs, vars, f.constant
end

function _parse_to_coeffs_vars(
    model::Optimizer,
    terms::Vector{MOI.ScalarAffineTerm{T}},
) where {T <: Integer}
    coeffs = Int32[t.coefficient for t in terms]
    vars = Variable[_info(model, t.variable) for t in terms]
    return coeffs, vars
end


text:
#!/usr/bin/env julia

using Luxor

using Test

using Colors

using Random
Random.seed!(42)

function draw_color_bars(x, y, n, colmap, label)
    setcolor(1, 1, 1)
    fontsize(8)

    number_of_colours = n
    text(label, x, y - 10)

    tilesize = (currentheight - 40 * 2) / number_of_colours

    for c in colmap
        setcolor(c)
        rect(x, y, tilesize, tilesize, :fill)
        setcolor((1, 1, 1))
        text("$(c.r), $(c.g), $(c.b)", x + tilesize + 5, y + tilesize/2)
        y +=tilesize
    end
end

function palette_test(fname)
    global currentwidth = 1000 # pts
    global currentheight = 1000 # pts
    Drawing(currentwidth, currentheight, fname)

    # background
    setcolor("grey50")
    rect(0,0,currentwidth, currentheight, :fill)

    # sequential_palette(h, [N::Int=100; c=0.88, s=0.6, b=0.75, w=0.15, d=0.0, wcolor=RGB(1,1,0), dcolor=RGB(0,0,1), logscale=false])
    # diverging_palette(h1, h2 [, N::Int=100; mid=0.5,c=0.88, s=0.6, b=0.75, w=0.15, d1=0.0, d2=0.0, wcolor=RGB(1,1,0), dcolor1=RGB(1,0,0), dcolor2=RGB(0,0,1), logscale=false])
    # colormap(cname::AbstractString [, N::Int=100; mid=0.5, logscale=false, kvs...])

    draw_color_bars(20, 20, 50, diverging_palette(0, 360, 50, mid=0.1), "diverging 0 360 mid 0.1")
    draw_color_bars(120, 20, 50, sequential_palette(0, 50), "sequential 0 #50")
    draw_color_bars(220, 20, 50, sequential_palette(128, 50), "sequential 128 #50")
    draw_color_bars(320, 20, 50, sequential_palette(200, 50), "sequential 200 #50")
    draw_color_bars(420, 20, 50, diverging_palette(120, 140, 50, mid=0.1), "diverging 120 140 #50 ")
    draw_color_bars(520, 20, 50, distinguishable_colors(50), "distinguishable")
    draw_color_bars(620,  20, 50, diverging_palette(120, 140, 50, b=1.0), "diverging 120 140 b 0.5")
    draw_color_bars(720,  20, 50, diverging_palette(120, 140, 50, mid=1.0), "diverging 120 140 mid 1.0")

    @test finish() == true
    println("...finished test: output in $(fname)")
end

palette_test("color-palette.pdf")


text:
module RandomMatrixDistributions

#import RandomMatrices

using Random, Distributions
using LinearAlgebra
using BandedMatrices
using KrylovKit

import LinearAlgebra: eigmax

import Distributions: minimum, maximum, quantile

import Distributions: ContinuousUnivariateDistribution,
    ContinuousMatrixDistribution,
    cdf, pdf, entropy, insupport, mean, median, modes, kurtosis, skewness, std, var, moment

export randeigvals, randeigstat,
    minimum, maximum, quantile
    cdf, pdf, entropy, insupport, mean, median, modes, kurtosis, skewness, std, var, moment,
    eigmax

# Generic eigenvalue sampling functions

function randeigvals(d)
    eigvals(randreduced(d))
end

function randeigstat(d, eigstat, nsims::Int)
    statvals = Array{Float64, 1}(undef, nsims)
    for i in 1:nsims
        statvals[i] = eigstat(randreduced(d))
    end
    statvals
end

include("BandedHelpers.jl")

include("SpikedWigner.jl")
include("SpikedWishart.jl")

include("Jacobi.jl")

include("densities/MarchenkoPastur.jl")
include("densities/TracyWidom.jl")
include("densities/Wachter.jl")

end


text:
module TestClustering

# using Revise
using MLJ
using Test
using Random:seed!
using LinearAlgebra:norm
using Distances: evaluate

# load code to be tested:
import MLJModels
import Clustering # MLJModels.Clustering_ now available for loading
using MLJModels.Clustering_

seed!(132442)

task = load_crabs()

X, y = X_and_y(task)

####
#### KMEANS
####

barekm = KMeans()

fitresult, cache, report = MLJ.fit(barekm, 1, X)

R = MLJ.transform(barekm, fitresult, X)

X_array = convert(Matrix{Float64}, X)

# distance from first point to second center
@test R[1, 2] ≈ norm(view(X_array, 1, :) .- view(fitresult, :, 2))^2
@test R[10, 3] ≈ norm(view(X_array, 10, :) .- view(fitresult, :, 3))^2

p = MLJ.predict(barekm, fitresult, X)

R_mat = MLJ.matrix(R)

@test argmin(R_mat[1, :]) == p[1]
@test argmin(R_mat[10, :]) == p[10]

km = machine(barekm, X)

fit!(km)

info(barekm)

####
#### KMEDOIDS
####

barekm = KMedoids()

fitresult, cache, report = MLJ.fit(barekm, 1, X)

R = MLJ.matrix(MLJ.transform(barekm, fitresult, X))

@test R[1, 2] ≈ evaluate(barekm.metric, view(X_array, 1, :), view(fitresult, :, 2))
@test R[10, 3] ≈ evaluate(barekm.metric, view(X_array, 10, :), view(fitresult, :, 3))

p = MLJ.predict(barekm, fitresult, X)

@test all(report[:assignments] .== p)

km = machine(barekm, X)

fit!(km)

info(barekm)

end # module
true


text:
Base.@kwdef mutable struct SARSA{O, T<:AbstractTraceUpdate} <: LearningUpdate
    lambda::Float64
    opt::O
    trace::T
    e::IdDict = IdDict()
    prev_discounts::IdDict = IdDict()
end

function update!(lu::SARSA,
                 learner::QLearner{M, LU},
                 demons,
                 obs,
                 next_obs,
                 state,
                 action,
                 next_state,
                 next_action,
                 is_terminal,
                 behaviour_pi_func,
                 reward) where {M<:AbstractMatrix, LU<:SARSA}

     weights = learner.model
     λ = lu.lambda
     e = get!(()->zero(weights), lu.e, weights)::typeof(weights)

     C, next_discounts, target_pis, next_target_pis = get_demon_parameters(lu, learner, demons, obs, state, action, next_obs, next_state, next_action, reward)
     b_πs = behaviour_pi_func(state, obs)
     discounts = get!(()->zero(next_discounts), lu.prev_discounts, learner)::typeof(next_discounts)

     #NOTE: Off policy SARSA(λ) I think is the same as EASRSA except the backup is sample and not expected.
     ρ = if b_πs[action] == 0
         zeros(length(target_pis[:,action]))
     else
         target_pis[:,action] ./ b_πs[action]
     end

     inds = get_action_inds(action, learner.num_actions, learner.num_demons)
     state_action_row_ind = inds

     update_trace!(lu.trace,
                   e,
                   state,
                   λ,
                   repeat(discounts, inner = learner.num_actions),
                   repeat(ρ, inner = learner.num_actions),
                   inds)

     next_preds = learner(next_state, next_action)
     pred = learner(state, action)

     td_err = C .+ next_discounts .* next_preds .- pred
     td_err_across_demons = repeat(vec(td_err), inner=learner.num_actions)

     if lu.opt isa Auto
         next_state_action_row_ind = get_action_inds(next_action, learner.num_actions, learner.num_demons)
         state_discount = zero(e)
         state_discount[state_action_row_ind,:] .+= state'
         state_discount[next_state_action_row_ind,:] .-= next_discounts * next_state'
         abs_phi = abs.(e)
         update!(lu.opt, weights, e, td_err_across_demons, abs_phi .* max.(state_discount, abs_phi))
     else
         Flux.Optimise.update!(lu.opt, weights,  -(e .* td_err_across_demons))
     end
     discounts .= next_discounts
end

function zero_all_nonfinites!(itm)
    if eltype(itm) <: Number
        itm[findall(.!isfinite.(itm))] .= 0.0

    elseif eltype(itm) <: Vector
        [zero_all_nonfinites(lst) for lst in itm]
    else
        throw(ArgumentError("Not Implemented"))
    end
end

function update!(lu::SARSA,
                 learner::LSTDLearner,
                 demons,
                 obs,
                 next_obs,
                 state,
                 action,
                 next_state,
                 next_action,
                 is_terminal,
                 behaviour_pi_func,
                 env_reward)

     C, γ_tp1, π_t, π_tp1 =
         get_demon_parameters(learner,
                              demons,
                              obs,
                              state,
                              action,
                              next_obs,
                              next_state,
                              next_action)

     na = learner.num_actions
     fs = learner.feature_size
     λ = lu.lambda
     t = learner.t
     μ_t = behaviour_pi_func(state, obs)[action]
     μ_tp1 = behaviour_pi_func(next_state, next_obs)

     ρ_t = π_t ./ μ_t
     ρ_tp1 = π_tp1 ./ μ_tp1'

     zero_all_nonfinites!(ρ_t)
     zero_all_nonfinites!(ρ_tp1)

     γ_t = get!(()->zero(γ_tp1), lu.prev_discounts, learner)::typeof(γ_tp1)

     ϕ_t = state
     ϕ_tp1 = next_state

     x_t = get_active_action_state_vector(
         ϕ_t, action, learner.feature_size, learner.num_actions)
     x_tp1 = get_active_action_state_vector(
         ϕ_tp1, next_action, learner.feature_size, learner.num_actions)

    all_gvf_e = get!(()->zero.(learner.w_real), lu.e, learner.w_real)::typeof(learner.w_real)

     for gvf ∈ 1:learner.num_demons
         A_inv = learner.A_inv[gvf]
         e = all_gvf_e[gvf]
         b = learner.b[gvf]
         c = C[gvf]

         # NOTE TD Update
         e .= γ_t[gvf]*λ*ρ_t[gvf] * e + x_t
         #NOTE TB Update
         # e .= γ_t[gvf]*λ*π_t[gvf] * e + x_t
         b .+= (c*e - b)/(t+1)

         # NOTE TD UPDATE
         u = sum(ρ_tp1[gvf, a] .* get_active_action_state_vector(ϕ_tp1, a, fs, na) for a ∈ 1:na)
         # NOTE TB Update
         # u = sum(π_tp1[gvf, a] .* get_active_action_state_vector(ϕ_tp1, a, fs, na) for a ∈ 1:na)
         v = transpose(transpose(x_t - γ_tp1[gvf]*u) * A_inv)

         if t > 0
             scale = (t+1)/t
             vz = dot(v, e)
             Ainv_zvt = (A_inv * e) * transpose(v)
             A_inv .= scale * (A_inv - Ainv_zvt./(t + vz))
         else
             Aev = A_inv*(e*v')
             ve = dot(v, e)
             A_inv .-= Aev/(1 + ve)
         end
         learner.w_real[gvf] .= A_inv * b
     end
     γ_t .= γ_tp1
     learner.t += 1
end

function zero_eligibility_traces!(lu::SARSA)
    for (k, v) ∈ lu.e
        if eltype(v) <: Integer
            lu.e[k] = Int[]
        elseif eltype(v) <: Vector #Used for LSTD
            lu.e[k] = zero.(v)
        else
            v .= 0
        end
    end
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
"""
Load data from disk into a vector of tuples of the form (X2, T1, X1, T1, Y)
and a vector of frequency counts. This compressed representation of the data
allows approximately a 10x reduction in processing time because redundant
computations are avoided.
"""
function load()
    # Load raw data from CSV file.
    df = DataFrames.readtable("data/choices.csv")

    # Drop nulls.
    df = df[!DataFrames.isna(df[:LaterOptionChosen]), :]

    # Convert to non-nullable arrays.
    X = hcat(
        convert(Array{Float64}, df[:X2]),
        convert(Array{Float64}, df[:T2]),
        convert(Array{Float64}, df[:X1]),
        convert(Array{Float64}, df[:T1]),
    )
    y = convert(Array{Float64}, df[:LaterOptionChosen])

    # Normalize money measurements.
    X[:, 1] = X[:, 1] ./ maximum(X[:, 1])
    X[:, 3] = X[:, 3] ./ maximum(X[:, 3])

    # Translate data into a tuple format to maximize cache locality.
    n = size(X, 1)
    all_tuples = Array(NTuple{5, Float64}, n)
    for i in 1:n
        all_tuples[i] = (X[i, 1], X[i, 2], X[i, 3], X[i, 4], y[i])
    end

    # Work with the unique values in the input and their frequencies to
    # compress data for faster computations.
    unique_tuples_dict = StatsBase.countmap(all_tuples)
    unique_tuples = collect(keys(unique_tuples_dict))
    w = collect(values(unique_tuples_dict))

    return unique_tuples, w
end


text:
#include("libavdevice_h.jl")

include("avdevice.jl")


text:
using Printf

function do_broadcast()
    comm = MPI.COMM_WORLD

    if MPI.Comm_rank(comm) == 0
        println(repeat("-",78))
        println(" Running on $(MPI.Comm_size(comm)) processes")
        println(repeat("-",78))
    end

    MPI.Barrier(comm)

    N = 5
    root = 0

    if MPI.Comm_rank(comm) == root
        A = [1:N;] * (1.0 + im*2.0)
    else
        A = Array{ComplexF64}(undef, N)
    end

    MPI.Bcast!(A,length(A), root, comm)

    @printf("[%02d] A:%s\n", MPI.Comm_rank(comm), A)

    if MPI.Comm_rank(comm) == root
        B = Dict("foo" => "bar")
    else
        B = nothing
    end

    B = MPI.bcast(B, root, comm)
    @printf("[%02d] B:%s\n", MPI.Comm_rank(comm), B)

    if MPI.Comm_rank(comm) == root
        f = x -> x^2 + 2x - 1
    else
        f = nothing
    end

    f = MPI.bcast(f, root, comm)
    @printf("[%02d] f(3):%d\n", MPI.Comm_rank(comm), f(3))
end


text:
using PkgSegments

# Provide the directory to generate the segments from as the argument to
# this function.
PkgSegments.genfromsegfile!(ARGS[1])


text:
using GlobalSensitivity, Test, Distributions

samples = 2000000

f1(x) = x[1] + 2*x[2] + 6.00*x[3]
dist1 = [Uniform(4,10),Normal(4,23),Beta(2,3)]
b =  gsa(f1,DGSM(),dist1,samples=samples)
@test [b.a[1],b.absa[1],b.asq[1]] ≈ [1.0,1.0,1.0] atol=10e-1
@test [b.a[2],b.absa[2],b.asq[2]] ≈ [2.0,2.0,4.0] atol=10e-1
@test [b.a[3],b.absa[3],b.asq[3]] ≈ [6.0,6.0,36.0] atol=10e-1

#Test on classical Ishigami function
#Reference: https://hal.archives-ouvertes.fr/hal-01164215/document
#As provided in paper theoretical values of asq for the variables are
#asq for x[1] = 7.7
#asq for x[2] = 24.5
#asq for x[3] = 11.0

f2(x) = sin(x[1]) + 7sin(x[2])^2 + 0.1* (x[3]^4) *sin(x[1])
dist2 = [Uniform(-pi,pi),Uniform(-pi,pi),Uniform(-pi,pi)]
c =  gsa(f2,DGSM(),dist2,samples=samples)
@test [c.asq[1]] ≈ [7.7] atol= 10e-1
@test [c.asq[2]] ≈ [24.5] atol= 10e-1
@test [c.asq[3]] ≈ [11.0] atol= 10e-1
@test [c.sigma[1]] ≈ [-14.7] atol= 10e-1
@test [c.sigma[2]] ≈ [-38.7] atol= 10e-1
@test [c.sigma[3]] ≈ [-42.0] atol= 10e-1
@test [c.tao[1]] ≈ [6.13] atol= 10e-1
@test [c.tao[2]] ≈ [16.94] atol= 10e-1
@test [c.tao[3]] ≈ [15.91] atol= 10e-1

#Some test functions from a paper
#Reference:New methods for the sensitivity analysis of
#black-box functions with an application to
#sheet metal forming A doctoral Thesis by Jana Fruth

f3(x) = 0.5*sin(50*x[1])+0.5
dist3 = [Uniform(0,1)]
d =  gsa(f3,DGSM(),dist3,samples=samples)
@test [d.asq[1]] ≈ [310.9] atol= 10e-1

#Test on classical Ishigami function
#Reference: Roustant, Olivier, Jana Fruth, Bertrand Iooss and Sonja Kuhnt. “Derivative-Based Sensitivity Measures for Interactions.” (2013).
#As provided in paper theoretical values of asq for the variables are
#crossedsq for x1:x2 = 0.0
#crossedsq for x2:x3 = 0.0
#crossedsq for x1:x3 = 12.686

f2(x) = sin(x[1]) + 7*sin(x[2])^2 + 0.1* (x[3]^4) *sin(x[1])
dist2 = [Uniform(-pi,pi),Uniform(-pi,pi),Uniform(-pi,pi)]
c =  gsa(f2,DGSM(crossed=true),dist2,samples=samples)
@test [c.crossedsq[2]] ≈ [0.0] atol= 20e-1
@test [c.crossedsq[6]] ≈ [0.0] atol= 20e-1
@test [c.crossedsq[3]] ≈ [12.686] atol= 20e-1


text:
export ConnluSentence
export load_connlu_file

using TextAnalysis

struct ConnluSentence
  token_doc::TokenDocument
  string_doc::StringDocument
  pos_tags::Vector{String}
  gold_tree::DependencyTree
end

function load_connlu_file(filename::String)
  sentences = Vector{ConnluSentence}()

  open(filename) do file
    sentence_lines = []

    while !eof(file)
      line = readline(file)

      if !occursin(r"^\s*$", line)
        push!(sentence_lines, line)
      else
        convert(sentence_lines) |> sentence_data -> ConnluSentence(sentence_data...) |> sentence -> push!(sentences, sentence)
        sentence_lines = []
      end
    end
  end

  sentences
end

function convert(lines)
  root = TreeNode(0)
  words = filter(word -> (length(word) > 1) && !(occursin(r"^(# |(\d+\-\d+)|(\d+\.\d+))", word)), lines) .|>
      (word -> split(word, r"\t")) |>
      (words -> map(word -> [word[1], word[2], word[4], word[7], word[8]], words))

  nodes::Vector{TreeNode} = map(words) do word
    dependent_word_number = parse(Int32, word[1])
    head_word_number = parse(Int32, word[4])
    dependency = uppercase(word[5])

    TreeNode(dependent_word_number, String(dependency), head_word_number)
  end

  tree = DependencyTree(root, nodes, length(nodes))
  token_doc = TokenDocument(map(word -> String(word[2]), words))
  string_doc = findfirst(line -> occursin(r"# text = ", line), lines) |>
    index -> StringDocument(replace(lines[index], r"# text = " => ""))
  pos_tags = map(word -> String(word[3]), words)

  [token_doc, string_doc, pos_tags, tree]
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
struct Simulation{DT <: Number, ET <: Equation{DT}, IT <: Integrator{DT}}
    equ::ET
    int::IT
    x::OffsetArray{DT,3,Array{DT,3}}

    function Simulation(equ::ET, int::IT, nt::Int) where {DT, ET <: Equation{DT}, IT <: Integrator{DT}}
        x = OffsetArray(zeros(DT, size(equ.x₀)..., nt+1), axes(equ.x₀)..., 0:nt)
        x[:,:,0] .= equ.x₀
        new{DT,ET,IT}(equ, int, x)
    end
end

ndims(sim::Simulation) = ndims(sim.equ)
nsamples(sim::Simulation) = nsamples(sim.equ)
ntimesteps(sim::Simulation) = lastindex(sim.x,3)

eachsample(sim::Simulation) = axes(sim.x,2)
eachtimestep(sim::Simulation) = axes(sim.x,3)[1:end]

function run!(sim::Simulation{DT}) where {DT}
    for n in eachtimestep(sim)
        for i in eachsample(sim)
            # x₀ = @view sim.x[:,i,n-1]
            # x₁ = @view sim.x[:,i,n]
            # integrate_step!(sim.int, sim.equ, x₀, x₁)

            # integrate_step!(sim.int, sim.equ, view(sim.x, :, i ,n-1), view(sim.x, :, i, n))

            @views integrate_step!(sim.int, sim.equ, sim.x[:,i,n-1], sim.x[:,i,n])
        end
    end
    return sim.x
end

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
# Initialization

export device!, device_reset!

const initialized = Ref{Bool}(false)
const device_contexts = Dict{CuDevice,CuContext}()

# FIXME: support for flags (see `cudaSetDeviceFlags`)

# API calls that are allowed without lazily initializing the CUDA library
#
# this list isn't meant to be complete (ie. many other API calls are actually allowed
# without setting-up a context), and only serves to make multi-device applications possible.
#
# feel free to open a PR adding additional API calls, if you have a specific use for them.
const preinit_apicalls = Set{Symbol}([
    :cuDriverGetVersion,
    # device calls, commonly used to determine the most appropriate device
    :cuDeviceGet,
    :cuDeviceGetAttribute,
    :cuDeviceGetCount,
    :cuDeviceGetName,
    :cuDeviceTotalMem,
    # context calls, for testing
    :cuCtxGetCurrent
])

function maybe_initialize(apicall)
    initialized[] && return
    apicall in preinit_apicalls && return
    @debug "Initializing CUDA after call to $apicall"
    initialize()
end

function initialize(dev = CuDevice(0))
    # NOTE: we could do something smarter here,
    #       eg. select the most powerful device,
    #       or skip devices without free memory
    device!(dev)
end

const device!_listeners = Set{Function}()

"""
    device!(dev)

Sets `dev` as the current active device for the calling host thread. Devices can be
specified by integer id, or as a `CuDevice`. This is intended to be a low-cost operation,
only performing significant work when calling it for the first time for each device.

If your library or code needs to perform an action when the active device changes, add a
callback of the signature `(::CuDevice, ::CuContext)` to the `device!_listeners` set.
"""
function device!(dev::CuDevice)
    if !initialized[]
        initialized[] = true
        CUDAdrv.apicall_hook[] = nothing
    end

    ctx = get!(device_contexts, dev) do
        pctx = CuPrimaryContext(dev)
        CuContext(pctx)
        pop!(CuContext) # we don't maintain a stack...
    end
    activate(ctx)       # replace the top of the stack

    for listener in device!_listeners
        listener(dev, device_contexts[dev])
    end
end
device!(dev::Integer) = device!(CuDevice(dev))

"""
    device_reset!(dev::CuDevice=device())

Reset the CUDA state associated with a device. This call with release the underlying
context, at which point any objects allocated in that context will be invalidated.
"""
function device_reset!(dev::CuDevice=device())
    active_ctx = CuCurrentContext()
    if active_ctx == nothing
        active_dev = nothing
    else
        active_dev = device(active_ctx)
    end

    if haskey(device_contexts, dev)
        # invalidate compiled kernels
        ctx = device_contexts[dev]
        for id in collect(keys(compilecache))
            kernel = compilecache[id]
            if kernel.ctx == ctx
                delete!(compilecache, id)
            end
        end

        delete!(device_contexts, dev)
    end

    # unconditionally reset the primary context,
    # as there might be users outside of CUDAnative.jl
    pctx = CuPrimaryContext(dev)
    unsafe_reset!(pctx)

    if dev == active_dev
        # unless the user switches devices, new API calls should trigger initialization
        CUDAdrv.apicall_hook[] = maybe_initialize
        initialized[] = false
    end

    return
end

"""
    device!(f, dev)

Sets the active device for the duration of `f`.
"""
function device!(f::Function, dev::CuDevice)
    # FIXME: should use Push/Pop
    old_ctx = CuCurrentContext()
    try
        device!(dev)
        f()
    finally
        if old_ctx != nothing
            activate(old_ctx)
        end
    end
end
device!(f::Function, dev::Integer) = device!(f, CuDevice(dev))

function __init__()
    __init_compiler__()
    configured || return

    if CUDAdrv.version() != cuda_driver_version
        error("Your set-up has changed. Please run Pkg.build(\"CUDAnative\") and restart Julia.")
    end

    CUDAdrv.apicall_hook[] = maybe_initialize
end


text:
infile = "day19/input2.txt"
rules, text = split(read(infile, String), "\n\n")
rules = split(rules, "\n")
text = split(text, "\n")

r_dict = Dict()
for rule in rules
    number, matchlist = match(r"(\d+): (.+)", rule).captures
    number = parse(Int, number)
    if contains(matchlist, "\"")
        r_dict[number] = matchlist[2]
    elseif contains(matchlist, "|")
        matchlist1, matchlist2 = split.(strip.(split(matchlist, "|")))
        matchlist1 = parse.(Int, matchlist1)
        matchlist2 = parse.(Int, matchlist2)
        r_dict[number] = (matchlist1, matchlist2)
    else
        r_dict[number] = parse.(Int, split(matchlist))
    end
end

function matches_rule(rule, text)
    global r_dict
    if text === nothing || text == ""
        return []
    end
    if isa(rule, Tuple)
        return vcat([matches_rule(rule[1], text), matches_rule(rule[2], text)]...)
    elseif isa(rule, Array)
        rem_text = [text]
        for r in rule
            rem_text = vcat(map(t -> matches_rule(r_dict[r], t), rem_text)...)
        end
        return rem_text
    elseif isa(rule, Char)
        if text[1] == rule
            return [text[2:end]]
        else
            return []
        end
    end
end

results = map(t -> matches_rule(r_dict[0], t), text)
println(count(r -> in("", r), results))


text:
using Turing, Test
using Turing: reconstruct, invlink, step
using Turing.VarReplay
using Turing.VarReplay: uid, cuid, getvals, getidcs, set_retained_vns_del_by_spl!, is_flagged, unset_flag!

randr(vi::VarInfo, vn::VarName, dist::Distribution, spl::Turing.Sampler, count::Bool) = begin
  if ~haskey(vi, vn)
    r = rand(dist)
    Turing.push!(vi, vn, r, dist, spl.alg.gid)
    r
  elseif is_flagged(vi, vn, "del")
    unset_flag!(vi, vn, "del")
    r = rand(dist)
    Turing.setval!(vi, Turing.vectorize(dist, r), vn)
    r
  else
    if count Turing.checkindex(vn, vi, spl) end
    Turing.updategid!(vi, vn, spl)
    vi[vn]
  end
end

# Test for uid() (= string())
csym = gensym()
vn1 = VarName(csym, :x, "[1]", 1)
@test string(vn1) == "{$csym,x[1]}:1"
# println(string(vn1))

vn2 = VarName(csym, :x, "[1]", 2)
vn11 = VarName(csym, :x, "[1]", 1)

@test cuid(vn1) == cuid(vn2)
@test vn11 == vn1

vi = VarInfo()
dists = [Normal(0, 1), MvNormal([0; 0], [1.0 0; 0 1.0]), Wishart(7, [1 0.5; 0.5 1])]

alg = PG(PG(5,5),2)
spl2 = Turing.Sampler(alg)
vn_w = VarName(gensym(), :w, "", 1)
randr(vi, vn_w, dists[1], spl2, true)

vn_x = VarName(gensym(), :x, "", 1)
vn_y = VarName(gensym(), :y, "", 1)
vn_z = VarName(gensym(), :z, "", 1)
vns = [vn_x, vn_y, vn_z]

alg = PG(PG(5,5),1)
spl1 = Turing.Sampler(alg)
for i = 1:3
  r = randr(vi, vns[i], dists[i], spl1, false)
  val = vi[vns[i]]
  @test sum(val - r) <= 1e-9
end

# println(vi)

@test length(getvals(vi, spl1)) == 3
@test length(getvals(vi, spl2)) == 1

vn_u = VarName(gensym(), :u, "", 1)
randr(vi, vn_u, dists[1], spl2, true)

# println(vi)
vi.num_produce = 1
set_retained_vns_del_by_spl!(vi, spl2)

# println(vi)

vals_of_1 = collect(getvals(vi, spl1))
# println(vals_of_1)
filter!(v -> ~any(map(x -> isnan.(x), v)), vals_of_1)
@test length(vals_of_1) == 3

vals_of_2 = collect(getvals(vi, spl2))
# println(vals_of_2)
filter!(v -> ~any(map(x -> isnan.(x), v)), vals_of_2)
@test length(vals_of_2) == 1

@model gdemo() = begin
  x ~ InverseGamma(2,3)
  y ~ InverseGamma(2,3)
  z ~ InverseGamma(2,3)
  w ~ InverseGamma(2,3)
  u ~ InverseGamma(2,3)
end

# Test the update of group IDs
g_demo_f = gdemo()
g = Turing.Sampler(Gibbs(1000, PG(10, 2, :x, :y, :z), HMC(1, 0.4, 8, :w, :u)))

pg, hmc = g.info[:samplers]

vi = Turing.VarInfo()
g_demo_f(vi, nothing)
vi = step(g_demo_f, pg, vi)
@test vi.gids == [1,1,1,0,0]

g_demo_f(vi, hmc)
@test vi.gids == [1,1,1,2,2]


text:
function setup_RHS(cellxmax, cellymax, cellzmax, E_adv_hat, F_adv_hat, G_adv_hat, E_vis_hat, F_vis_hat, G_vis_hat)
    RHS = zeros(cellxmax, cellymax, cellzmax, 5)
    for i in 2:cellxmax-1
        for j in 2:cellymax-1
            for k in 2:cellzmax-1
                for l in 1:5
                    RHS[i,j,k,l] = -((E_adv_hat[i+1,j,k,l]-E_vis_hat[i+1,j,k,l]) - (E_adv_hat[i,j,k,l]-E_vis_hat[i,j,k,l]) +
                                     (F_adv_hat[i,j+1,k,l]-F_vis_hat[i,j+1,k,l]) - (F_adv_hat[i,j,k,l]-F_vis_hat[i,j,k,l]) +
                                     (G_adv_hat[i,j,k+1,l]-G_vis_hat[i,j,k+1,l]) - (G_adv_hat[i,j,k,l]-G_vis_hat[i,j,k,l]))
                end
            end
        end
    end
    return RHS
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
#!/usr/bin/env julia

using ConcurrentCalls

function countwords_serial(files)
    function countwords(file)
        counts = Dict{String,Int}()
        for word in eachmatch(r"\w+", readstring(file))
            counts[word.match] = get!(counts, word.match, 0) + 1
        end
        return counts
    end
    return [countwords(file) for file in files]
end

function countwords_concurrent(files)
    function countwords(file)
        counts = Dict{String,Int}()
        for word in eachmatch(r"\w+", readstring(file))
            counts[word.match] = get!(counts, word.match, 0) + 1
        end
        return counts
    end
    return @cc [countwords(file) for file in files]
end

files = readlines(STDIN)
if isempty(files)
    error("no input files")
end

println("--- statistics ---")
filesizes = map(filesize, files)
println("number of files:     $(length(files))")
println("mean of filesizes:   $(@sprintf("%.1f", mean(filesizes)))")
println("median of filesizes: $(@sprintf("%.1f", median(filesizes)))")
println("std. of filesizes:   $(@sprintf("%.1f", std(filesizes)))")
println()

println("--- countwords_serial ---")
countwords_serial(files[1:1])
for _ in 1:3
    @time countwords_serial(files)
end
println()

println("--- countwords_concurrent ---")
countwords_concurrent(files[1:1])
for _ in 1:3
    @time countwords_concurrent(files)
end

#= 2017/08/29 tostoint01
-bash-4.1$ find ~/anaconda3/ -type f -name "*.py" | head -3000 | julia -p4 examples/countwords.jl
--- statistics ---
number of files:     3000
mean of filesizes:   10472.5
median of filesizes: 3743.5
std. of filesizes:   28214.8

--- countwords_serial ---
  7.094024 seconds (28.20 M allocations: 1.252 GiB, 11.39% gc time)
  7.005876 seconds (28.20 M allocations: 1.252 GiB, 10.69% gc time)
  7.027203 seconds (28.20 M allocations: 1.252 GiB, 10.74% gc time)

--- countwords_concurrent ---
  3.208697 seconds (915.67 k allocations: 90.871 MiB, 2.09% gc time)
  2.668285 seconds (900.71 k allocations: 89.754 MiB, 4.32% gc time)
  2.911208 seconds (900.89 k allocations: 89.006 MiB, 3.81% gc time)
=#


text:
#!/usr/bin/env julia

using
  FourierFlows,
  Test,
  Statistics,
  Random,
  FFTW

import # use 'import' rather than 'using' for submodules to keep namespace clean
  PassiveTracerFlows.TracerAdvDiff

const rtol_traceradvdiff = 1e-12 # tolerance for rtol_traceradvdiff tests

# Run tests
testtime = @elapsed begin

@testset "TracerAdvDiff" begin
  include("test_traceradvdiff.jl")

  stepper = "RK4"
  dt, nsteps  = 1e-2, 40
  @test test_constvel(stepper, dt, nsteps)
  dt, tfinal  = 0.002, 0.1
  @test test_timedependentvel(stepper, dt, tfinal)
  dt, tfinal  = 0.005, 0.1
  @test test_diffusion(stepper, dt, tfinal; steadyflow=true)
  dt, tfinal  = 0.005, 0.1
  @test test_diffusion(stepper, dt, tfinal; steadyflow=false)
  dt, tfinal  = 0.005, 0.1
  @test test_hyperdiffusion(stepper, dt, tfinal)
end

end
println("Total test time: ", testtime)


text:
module LiquidsDynamics

### Imports

# Load submodules code
include("Projections.jl")
include("LiquidsStructureGrids.jl")

# Load the actual modules
using Reexport
using Parameters
using LinearAlgebra: I

@reexport using .Projections
@reexport using .LiquidsStructureGrids

### Exports
export asymptotics, asymptotics!, dynamics, dynamics!, initialize_asymptotics,
       initialize_dynamics

### Implementation
abstract type SCGLE end
abstract type Dynamics    <: SCGLE end
abstract type Asymptotics <: SCGLE end

include("utils.jl")
include("dynamics/dynamics.jl")
include("dynamics/init.jl")
include("dynamics/shorttimes.jl")
include("dynamics/intermediatetimes.jl")
include("dynamics/utils.jl")
include("asymptotics/asymptotics.jl")
include("asymptotics/init.jl")
include("asymptotics/utils.jl")

end # module


text:
module VOTables

export VOTable

using EzXML
using EzXML: Document, Node
using Tables

struct VOTable{X<:Document,R,T}
    doc::X
    resources::R
    tables::T
end

VOTable(source::AbstractString) = VOTable(parsexml(source))

function VOTable(xml::Document)
    # validate input is VOTable
    if nodename(root(xml)) != "VOTABLE"
        throw(ArgumentError("input XML does not appear to be a VOTable"))
    end
    ns = ["x" => namespace(root(xml))]
    resources = findall("x:RESOURCE", root(xml), ns)
    tables = map(VOTableData, findall("x:RESOURCE/x:TABLE", root(xml), ns))

    return VOTable(xml, resources, tables)
end

function Base.show(io::IO, vo::VOTable)
    n = length(vo.tables)
    print(io, "VOTable($n)")
end

Base.getindex(vo::VOTable, i) = VOTableData(vo.tables[i].node)

struct VOTableData{N,F,S,U,T,D}
    node::N
    fields::F
    names::S
    units::U
    types::T
    data::D
end

Base.names(tbl::VOTableData) = tbl.names
units(tbl::VOTableData) = tbl.units

struct VOTableDatum{T,U,S,D}
    name::S
    unit::U
    type::T
    datum::D
end

Tables.columnnames(row::AbstractVector{<:VOTableDatum}) = map(d -> Symbol(d.name), row)
function Tables.getcolumn(row::AbstractVector{<:VOTableDatum}, nm::Symbol)
    i = findfirst(==(nm), Tables.columnnames(row))
    return Tables.getcolumn(row, i)
end
Tables.getcolumn(row::AbstractVector{<:VOTableDatum}, i::Int) = row[i].datum

function VOTableData(node::Node)
    ns = ["x" => namespace(node)]
    fields = findall("x:FIELD", node, ns)
    names = map(f -> f["name"], fields)
    units = map(fields) do field
        atts = map(f -> nodename(f), attributes(field))
        return "unit" ∈ atts ? field["unit"] : missing
    end
    types = map(fields) do field
        atts = map(f -> nodename(f), attributes(field))
        datatype = field["datatype"]
        if "arraysize" ∈ atts
            if datatype == "char"
                return String
            else
                return Array{TYPE_MAP[datatype]}
            end
        end
        return TYPE_MAP[datatype]
    end
    data_nodes = findfirst("x:DATA/x:TABLEDATA", node, ns)
    rows = map(eachelement(data_nodes)) do row
        map(names, units, types, eachelement(row)) do name, unit, T, td
            VOTableDatum(name, unit, T, rowdata(T, nodecontent(td)))
        end
    end
    return VOTableData(node, fields, names, units, types, rows)
end

rowdata(::Type{<:AbstractArray{T}}, row) where {T} = safeparse.(eltype(T), split(row))
rowdata(::Type{T}, row) where {T} = safeparse(T, row)

function safeparse(T, val)
    out = tryparse(T, val)
    isnothing(out) ? missing : out
end

safeparse(::Type{String}, val) = isempty(val) ? missing : val

TYPE_MAP = Dict(
    "boolean" => Bool,
    "unsignedByte" => UInt8,
    "short" => Int16,
    "int" => Int32,
    "long" => Int64,
    "float" => Float32,
    "double" => Float64,
    "floatComplex" => ComplexF32,
    "doubleComplex" => ComplexF64,
    # "bit" => ,
    "char" => UInt8,
    "unicodeChar" => Char
)

Tables.schema(table::VOTableData) = Tables.Schema(table.names, table.types)

# Tables.jl interface
Tables.istable(::Type{<:VOTableData}) = true
Tables.rowaccess(::Type{<:VOTableData}) = true
Tables.rows(table::VOTableData) = getfield(table, :data)

end # module


text:
using Clustering, Combinatorics, DataArrays, DataFrames
println("Used Clustering, Combinatorics, DataArrays, DataFrames")
using DataStructures, Dierckx, HDF5, Ipopt, JLD, JSON, JuMP
println("Used DataStructures, Dierckx, HDF5, Ipopt, JLD, JSON, JuMP")
using MathProgBase, Nullables, SpecialFunctions, FactCheck, FunctionalData, Match, Syslogs, Memento, HttpCommon, HttpParser, IniFile, HTTP
println("Used MathProgBase, Nullables, SpecialFunctions, FactCheck,FunctionalData, Match, Syslogs, Memento, HttpCommon, HttpParser, IniFile, HTTP")

println("All packages used")

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr