text:
#=
Author : Guy Farmer
generates a distribution of ambulances for various ambulance counts and stores to a jld file
=#
using AmbulanceDeployment
#using DataFrames, Winston, JLD, CSV, Gurobi, JuMP, GLPK

function generate_deployment()
    hourly_calls = CSV.File(PROJECT_ROOT * "/test/austin-data/Full_WeekdayCalls.csv") |> DataFrame
    # weekend_hourly_calls = CSV.File("data/processed/2-weekend_calls.csv") |> DataFrame
    adjacent_nbhd = CSV.File(PROJECT_ROOT * "/test/austin-data/adjacent_nbhd.csv") |> DataFrame
    coverage = CSV.read(PROJECT_ROOT * "/test/austin-data/coverage_real.csv", DataFrame, header=false)
    coverage = convert(Array{Bool, 2}, coverage[:, :])
    incidents = CSV.File(PROJECT_ROOT * "/test/austin-data/austin_incidents.csv") |> DataFrame

    regions = Int[parse(Int,string(x)) for x in names(hourly_calls[:,6:ncol(hourly_calls)])]
    locations = collect(1:size(coverage,2))
    adjacent = convert(Array, adjacent_nbhd[:,2:ncol(adjacent_nbhd)])[regions,regions] .> 0.5
    demand = convert(Array,hourly_calls[:,6:end]);

    incidents = incidents[.~ismissing.(incidents[!,:stn1_min]), :] # drop 44 calls that were "unreachable" (because all stations are reachable from each other)
    incidents[!,:interarrival_seconds] = [0; incidents[!,:arrival_seconds][2:end] - incidents[!,:arrival_seconds][1:end-1]]
    incidents[!,:isweekday] .= true
    incidents[!,:isweekday][incidents[!,:dow] .== "Sun"] .= false
    incidents[!,:isweekday][incidents[!,:dow] .== "Sat"] .= false;
    incidents = incidents[incidents[!,:isweekday],:]

    regions2index = Dict{Int,Int}(regions[i]=>i for i in 1:length(regions))

#= previously there was an error that didnt account for all the regions because it
it only included regions for incidents. the regions without incidents are changed
to 0=#
    for x in incidents[!,:neighborhood]
        if(!haskey(regions2index,x))
            regions2index[x] = 0
        end
    end
    incidents[!,:neighborhood] = [regions2index[x] for x in incidents[!,:neighborhood]];

    calls = incidents[:,[[:hour,:dow,:month,:year,:neighborhood,:interarrival_seconds];
    Symbol[Symbol("stn$(i)_min") for i in locations]]]

    DataFrames.first(calls, 6)

# We focus on emergency calls during the "peak period" (8AM - 8PM),
# with the emergency calls from the first 3 month as our training set,
# and the subsequent emergency calls from the remaining months as our test set
    peak_period = (hourly_calls[!,:hour] .>= 8) .* (hourly_calls[!,:hour] .<= 20)
    indices = 1:DataFrames.nrow(hourly_calls);
    train_filter = (hourly_calls[!,:year] .== 2019) .* (hourly_calls[!,:month] .<= 3)
    test_filter  = .~train_filter;
    train_indices = indices[train_filter]
    test_indices = indices[test_filter];

# Same as for the hourly calls; but this is for individual emergency calls
    inc_peak_period = (calls[!,:hour] .>= 8) .* (calls[!,:hour] .<= 20)
    inc_indices = 1:DataFrames.nrow(calls);

    inc_train_filter = (calls[!,:year] .== 2019) .* (calls[!,:month] .<= 3)
    inc_test_filter  = .~inc_train_filter

    inc_train_indices = inc_indices[inc_train_filter]
    inc_test_indices = inc_indices[inc_test_filter];

# we distinguish between peak and offpeak hours
    train_peak = indices[peak_period .* train_filter]
    train_offpeak = indices[.~peak_period .* train_filter]

    test_peak = indices[peak_period .* test_filter]
    test_offpeak = indices[.~peak_period .* test_filter]

    p = DeploymentProblem(30, length(locations), length(regions), demand, train_indices,
          test_indices, coverage[regions,:], Array{Bool,2}(adjacent));

    test_inc_peak = inc_peak_period .* inc_test_filter
    test_inc_offpeak = .~inc_peak_period .* inc_test_filter;

## make sure to include @time begin

scenarios = Dict{Symbol, Dict{Int, Vector{Vector{Int}}}}()
generated_deployment = Dict{Symbol, Dict{Int, Vector{Vector{Int}}}}()
upperbounds = Dict{Symbol, Dict{Int, Vector{Float64}}}()
lowerbounds = Dict{Symbol, Dict{Int, Vector{Float64}}}()
upptiming = Dict{Symbol, Dict{Int, Vector{Float64}}}()
lowtiming = Dict{Symbol, Dict{Int, Vector{Float64}}}()
amb_deployment = Dict{Symbol, Dict{Int, Vector{Int}}}()
for(deployment_model, name) in ((dp -> RobustDeployment(dp, α=0.1), :Robust01),
                                  (dp -> RobustDeployment(dp, α=0.05), :Robust005),
                                  (dp -> RobustDeployment(dp, α=0.01), :Robust001),
                                  (dp -> RobustDeployment(dp, α=0.001), :Robust0001),
                                  (dp -> RobustDeployment(dp, α=0.0001), :Robust00001))
    println("$name: ")
    amb_deployment[name] = Dict{Int, Vector{Int}}()
    scenarios[name] = Dict{Int, Vector{Vector{Int}}}()
    generated_deployment[name] = Dict{Int, Vector{Vector{Int}}}()
    upperbounds[name] = Dict{Int, Vector{Float64}}()
    lowerbounds[name] = Dict{Int, Vector{Float64}}()
    upptiming[name] = Dict{Int, Vector{Float64}}()
    lowtiming[name] = Dict{Int, Vector{Float64}}()
    for namb in 30:5:50
        println("$namb ")
        p.nambulances = namb
        model = deployment_model(p)
        set_optimizer(model.m, Gurobi.Optimizer)
        # solve(model, p)
        @time AmbulanceDeployment.optimize!(model, p)
        #print("($(toq())) ")
        amb_deployment[name][namb] = deployment(model)

        # for tracking purposes
        scenarios[name][namb] = model.scenarios
        generated_deployment[name][namb] = model.deployment
        upperbounds[name][namb] = model.upperbounds
        lowerbounds[name][namb] = model.lowerbounds
        upptiming[name][namb] = model.upptiming
        lowtiming[name][namb] = model.lowtiming
    end
    println
end

<<<<<<< Updated upstream
    # for (next_deployment_model, name) in ((next_dp -> StochasticDeployment(next_dp, nperiods=500), :Stochastic),
    #                              (next_dp -> MEXCLPDeployment(next_dp, 0.654), :MEXCLP),
    #                              (next_dp -> MALPDeployment(next_dp, 0.654), :MALP))
    # println("$name: ")
    # amb_deployment[name] = Dict{Int, Vector{Int}}()
    #     for namb in 25:5:50
    #         println("$namb ")
    #         p.nambulances = namb
    #         next_model = next_deployment_model(p)
    #         set_optimizer(next_model.m, Gurobi.Optimizer)
    #         @time AmbulanceDeployment.optimize!(next_model, p)
    #         amb_deployment[name][namb] = deployment(next_model)
    #     end
    #     println()
    # end
    solver_stats = Dict{String, Any}()
    push!(solver_stats, "amb_deployment" => amb_deployment)
    push!(solver_stats, "scenarios" => scenarios)
    push!(solver_stats, "generated_deployment" => generated_deployment)
    push!(solver_stats, "upperbounds" => upperbounds)
    push!(solver_stats, "lowerbounds" => lowerbounds)
    push!(solver_stats, "upptiming" => upptiming)
    push!(solver_stats, "lowtiming" => lowtiming)
    json_string = JSON.json(solver_stats)

    open("src/outputs/solver_stats.json","w") do f
        write(f, json_string)
    end

end
=======
         for (next_deployment_model, name) in ((next_dp -> StochasticDeployment(next_dp, nperiods=500), :Stochastic),
                                         (next_dp -> MEXCLPDeployment(next_dp, 0.654), :MEXCLP),
                                         (next_dp -> MALPDeployment(next_dp, 0.654), :MALP))
            println("$name: ")
            amb_deployment[name] = Dict{Int, Vector{Int}}()
            for namb in 25:5:50
                println("$namb ")
                p.nambulances = namb
                next_model = next_deployment_model(p)
                set_optimizer(next_model.m, Gurobi.Optimizer)
                @time AmbulanceDeployment.optimize!(next_model, p)
                amb_deployment[name][namb] = deployment(next_model)
            end
            println()
        end
solver_stats = Dict{String, Any}()
push!(solver_stats, "amb_deployment" => amb_deployment)
push!(solver_stats, "scenarios" => scenarios)
push!(solver_stats, "generated_deployment" => generated_deployment)
push!(solver_stats, "upperbounds" => upperbounds)
push!(solver_stats, "lowerbounds" => lowerbounds)
push!(solver_stats, "upptiming" => upptiming)
push!(solver_stats, "lowtiming" => lowtiming)
json_string = JSON.json(solver_stats)

open("src/outputs/solver_stats.json","w") do f
                   write(f, json_string)
               end
>>>>>>> Stashed changes

max_stars_repo_name: michaelhilborn/AmbulanceDeployement

max_issues_repo_name: joshua-ong/AmbulanceDeployment

max_forks_repo_name: joshua-ong/AmbulanceDeployment


text:
"""

"""
abstract type ParameterSampler end

"""

"""
function sample(ps::ParameterSampler, parameters::NamedTuple)
    sample(ps, values(parameters)...)
end

"""

"""
struct CartesianParameterSampler <: ParameterSampler end

function sample(::CartesianParameterSampler, parameters::Vararg{Parameter,N}) where {N}
    # make sure all parameters have a sample vector
    for p in parameters
        @assert hassamples(p)
    end

    # get all parameter index combinations
    inds = CartesianIndices(zeros([length(p) for p in parameters]...))[:]

    # generate sample matrix
    smps = [parameters[i].samples[inds[j][i]] for j in eachindex(inds), i in 1:N]

    sinds = Tuple(p.name for p in parameters)
    svals = Tuple(smps[:,j] for j in axes(smps,2))

    Table(; NamedTuple{sinds}(svals)...)
end

ext: jl


text:
mutable struct Table{A,B} # A!=B
	f::Dict{A,B}
	inv::Dict{B,A}
end
function Table{A,B}()where {A,B}
	return Table{A,B}(Dict{A,B}(),Dict{B,A}())
end
function insert!(t::Table{A,B},a::A,b::B)where {A,B}
	t.f[a]=b
	t.inv[b]=a
end
function Base.haskey(t::Table{A,B},a::A)where {A,B}
	return haskey(t.f,a)
end
function Base.haskey(t::Table{A,B},b::B)where {A,B}
	return haskey(t.inv,b)
end
function Base.getindex(t::Table{A,B},a::A)where {A,B}
	return t.f[a]
end
function Base.getindex(t::Table{A,B},b::B)where {A,B}
	return t.inv[b]
end


text:
using Distributions

function mc_sample_path(P; init = 1, sample_size = 1000)
    @assert size(P)[1] == size(P)[2] # square required
    N = size(P)[1]; # should be square

    # create vector of discrete RVs for each row of stochastic matrix P
    dists = [Categorical(P[i, :]) for i in 1:N]

    # set up Simulation
    X = fill(0, sample_size); # allocate memory, or zeros(Int64, sample_size)
    X[1] = init # set the initial state

    for t in 2:sample_size
        dist = dists[X[t - 1]];
        X[t] = rand(dist);
    end

    return X
end


text:
for o in [:Scatter, :Bar, :Scattergl, :Pie, :Heatmap, :Image, :Contour, :Table, :Box, :Violin,
            :Histogram, :Histogram2dContour, :Ohlc, :Candlestick, :Waterfall, :Funnel, :Funnelarea,
            :Indicator, :Scatter3d, :Surface, :Mesh3d, :Cone, :Volume, :Isosurface, :Scattergeo,
            :Choropleth, :Scattermapbox, :Choroplethmapbox, :Densitymapbox,
            :Scatterpolar, :Scatterpolargl, :Barpolar, :Scatterternary, :Sunburst, :Treemap, :Sankey,
            :Splom, :Parcats, :Parcoords, :Carpet, :Scattercarpet, :Contourcarpet, :Layout]
    @eval begin
        export $o
        function $o(args...;kwargs...)
            plotly.graph_objects.$o(args...;kwargs...)
        end
    end
end

export Plot
"""
    Plot(rows::Int64 = 1, cols::Int64 = 1, args...; kwargs...)

Makes a figure consists of `rows × cols` subplots.

# Example
```julia
fig = Plot(3,1)
x = LinRange(0,1,100)
y1 = sin.(x)
y2 = sin.(2*x)
y3 = sin.(3*x)
fig.add_trace(Scatter(x=x, y=y1, name = "Line 1"), row = 1, col = 1)
fig.add_trace(Scatter(x=x, y=y2, name = "Line 2"), row = 2, col = 1)
fig.add_trace(Scatter(x=x, y=y3, name = "Line 3"), row = 3, col = 1)
fig.show()
```
"""
function Plot(rows::Int64 = 1, cols::Int64 = 1, args...; kwargs...)
    pkgs = read(`$PIP list`, String)
    if !occursin("plotly", pkgs)
        run(`$PIP install plotly==4.14.3`)
    end
    copy!(plotly, "plotly")
    spt = pyimport("plotly.subplots")
    spt.make_subplots(rows = rows, cols = cols, args...;kwargs...)
end

export to_html
"""
    to_html(fig::PyObject, filename::String, args...;
    include_plotlyjs = "cnd",
    include_mathjax = "cnd",
    full_html = true,
    kwargs...)

Exports the figure `fig` to an HTML file.
"""
function to_html(fig::PyObject, filename::String, args...;
        include_plotlyjs = "cdn",
        include_mathjax = "cdn",
        full_html = true,
        kwargs...)
    open(filename, "w") do io
        cnt = fig.to_html(
            include_plotlyjs = include_plotlyjs,
            include_mathjax = include_mathjax,
            full_html = full_html
        )
        write(io, cnt)
    end
end


text:
# Autogenerated wrapper script for x264_jll for x86_64-w64-mingw32
export libx264, x264

JLLWrappers.@generate_wrapper_header("x264")
JLLWrappers.@declare_library_product(libx264, "libx264-163.dll")
JLLWrappers.@declare_executable_product(x264)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libx264,
        "bin\\libx264-163.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        x264,
        "bin\\x264.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module DBGraphs

# package code goes here

using LightGraphs
using JuliaDB
import JuliaDB: DTable
import LightGraphs: AbstractGraph, vertices, edges, neighbors, indegree, outdegree,
       in_neighbors, out_neighbors, nv, ne, is_directed, eltype, is_directed

export DBGraph

abstract type AbstractDBGraph <: LightGraphs.AbstractGraph end

mutable struct DBGraph{T} <: AbstractDBGraph
  table::T
  src::Symbol
  dst::Symbol
  directed::Bool
end

function DBGraph(table::DTable)
  return DBGraph(table, :src, :dst, true)
end

#is_directed(g::DBGraph) = g.directed
#is_directed(::DBGraph) = true
#is_directed(::Type{DBGraph}) = true
is_directed(::Type{DBGraph{T}}) where T = false
eltype(::DBGraph) = Int

function vertices(g::DBGraph)
  return 1:nv(g) #select(g.table, g.dst, agg=first)
end

colmax(table, column) = reduce(max,map(x->getfield(x, column), table))
nv(g::DBGraph) = max(colmax(g.table, g.src), colmax(g.table, g.dst))
ne(g::DBGraph) = size(g.table,1)

# function vertices(g::DBGraph)
#   return select(g.table, g.dst, agg=first)
# end

function edges(g::DBGraph)
  return select(g.table, [g.src,g.dst])
end

# function out_neighbors(g::DBGraph, v)
#   return g.table[v,:]
# end
# function in_neighbors(g::DBGraph, v)
#   return g.table[:,v]
# end
# function neighbors(g::DBGraph, v)
#   return merge(out_neighbors(g,v), in_neighbors(g, v))
# end

function out_neighbors(g::DBGraph, v::Integer)
  return collect(map(x->x.dst, filter(x->x[1] == v, g.table)))
end
function in_neighbors(g::DBGraph, v::Integer)
  return collect(map(x->x.src, filter(x->x[2] == v, g.table)))
end
function neighbors(g::DBGraph, v::Integer)
  return union(out_neighbors(g,v), in_neighbors(g, v))
end

indegree(g::DBGraph, v::Integer) = length(collect(in_neighbors(g, v)))
outdegree(g::DBGraph, v::Integer) = length(collect(out_neighbors(g, v)))

end # module


text:
# Overload the norm
LinearAlgebra.norm(x, p, λ) = norm(x[abs.(x) .> λ], p)

# Dispatch on vector
function (opt::AbstractOptimizer{T} where T)(
        X::AbstractArray, A::AbstractArray, Y::AbstractVector{V} where V,
        args...; kwargs...
    )
    Y = reshape(Y, 1, length(Y))
    return opt(X, A, Y, args...; kwargs...)
end

function (opt::AbstractOptimizer{T} where T)(
        X::AbstractArray, A::AbstractArray, Y::Adjoint{V, AbstractVector{V}} where V,
        args...; kwargs...
    )
    Y = reshape(Y, 1, length(Y))
    return opt(X, A, Y, args...; kwargs...)
end

"""
$(SIGNATURES)

Set the threshold(s) of an optimizer to (a) specific value(s).
"""
function set_threshold!(opt::AbstractOptimizer{T}, threshold::T) where T
    @assert all(threshold .> zero(T)) "Threshold must be positive definite"
    opt.λ .= threshold
end

"""
$(SIGNATURES)

Get the threshold(s) of an optimizer.
"""
get_threshold(opt::AbstractOptimizer) = opt.λ

"""
$(SIGNATURES)

Initialize the optimizer with the least square solution for explicit and `zeros` for implicit optimization.
"""
CommonSolve.init(o::AbstractOptimizer, A, Y) = A \ Y

CommonSolve.init(o::AbstractSubspaceOptimizer, A, Y) = zeros(eltype(A), size(A, 2), size(Y, 2))

"""
$(SIGNATURES)

Initialize the optimizer with the least square solution for explicit and `zeros` for implicit optimization in place.
"""
init!(X::AbstractArray, o::AbstractOptimizer, A::AbstractArray, Y::AbstractArray) =  begin
    @static if VERSION < v"1.7.0"
        ldiv!(X, qr(A, Val(true)), Y)
    else
        ldiv!(X, qr(A, ColumnNorm()), Y)
    end
end

include("./utils.jl")

include("./proximals.jl")

# Remove the trace right now

include("./stlsq.jl")
include("./admm.jl")
include("./sr3.jl")

#Nullspace for implicit sindy
include("./implicit.jl")

# Init the progressmeters
# For a general optimizer
default_prg_msg(o::AbstractOptimizer) = summary(o)

function init_progress(opt::AbstractOptimizer, X, A, Y, maxiters, start)
    Progress(
        maxiters, default_prg_msg(opt), start
    )
end

include("./sparseregression.jl")


text:
module OCCAHelpersTest
using Base.Test
using OCCAHelpers

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

function _create_binary_switch(ids, exprs)
    if length(exprs) <= 3
        out = Expr(:if, Expr(:call, :(==), :id, ids[1]), exprs[1])
        if length(exprs) > 1
            push!(out.args, _create_binary_switch(ids[2:end], exprs[2:end]))
        end
        return out
    else
        mid = length(exprs) >>> 1
        return Expr(
            :if,
            Expr(:call, :(<=), :id, ids[mid]),
            _create_binary_switch(ids[1:mid], exprs[1:mid]),
            _create_binary_switch(ids[mid+1:end], exprs[mid+1:end]),
        )
    end
end

# We use a let block here for `expr` to create a local variable that does not
# persist in the scope of the module. All we care about is the _eval_univariate
# function that is eval'd as a result.
let exprs = map(SYMBOLIC_UNIVARIATE_EXPRESSIONS) do arg
        return :(return $(arg[1])(x), $(arg[2]))
    end
    @eval @inline function _eval_univariate(id, x::T) where {T}
        $(_create_binary_switch(1:length(exprs), exprs))
        return error("Invalid operator_id")
    end
end

# We use a let block here for `expr` to create a local variable that does not
# persist in the scope of the module. All we care about is the
# _eval_univariate_2nd_deriv function that is eval'd as a result.
let exprs = map(SYMBOLIC_UNIVARIATE_EXPRESSIONS) do arg
        if arg === :(nothing)  # f''(x) isn't defined
            :(error("Invalid operator_id"))
        else
            :(return $(arg[3]))
        end
    end
    @eval @inline function _eval_univariate_2nd_deriv(id, x::T) where {T}
        $(_create_binary_switch(1:length(exprs), exprs))
        return error("Invalid operator_id")
    end
end

struct _UnivariateOperator{F,F′,F′′}
    f::F
    f′::F′
    f′′::F′′
end

struct _MultivariateOperator{F,F′,F′′}
    N::Int
    f::F
    ∇f::F′
    ∇²f::F′′
    function _MultivariateOperator{N}(
        f::Function,
        ∇f::Function,
        ∇²f::Union{Nothing,Function} = nothing,
    ) where {N}
        return new{typeof(f),typeof(∇f),typeof(∇²f)}(N, f, ∇f, ∇²f)
    end
end

"""
    DEFAULT_UNIVARIATE_OPERATORS

The list of univariate operators that are supported by default.
"""
const DEFAULT_UNIVARIATE_OPERATORS = first.(SYMBOLIC_UNIVARIATE_EXPRESSIONS)

"""
    DEFAULT_MULTIVARIATE_OPERATORS

The list of multivariate operators that are supported by default.
"""
const DEFAULT_MULTIVARIATE_OPERATORS = [:+, :-, :*, :^, :/, :ifelse]

"""
    OperatorRegistry()

Create a new `OperatorRegistry` to store and evaluate univariate and
multivariate operators.
"""
struct OperatorRegistry
    # NODE_CALL_UNIVARIATE
    univariate_operators::Vector{Symbol}
    univariate_operator_to_id::Dict{Symbol,Int}
    univariate_user_operator_start::Int
    registered_univariate_operators::Vector{_UnivariateOperator}
    # NODE_CALL_MULTIVARIATE
    multivariate_operators::Vector{Symbol}
    multivariate_operator_to_id::Dict{Symbol,Int}
    multivariate_user_operator_start::Int
    registered_multivariate_operators::Vector{_MultivariateOperator}
    # NODE_LOGIC
    logic_operators::Vector{Symbol}
    logic_operator_to_id::Dict{Symbol,Int}
    # NODE_COMPARISON
    comparison_operators::Vector{Symbol}
    comparison_operator_to_id::Dict{Symbol,Int}
    function OperatorRegistry()
        univariate_operators = copy(DEFAULT_UNIVARIATE_OPERATORS)
        multivariate_operators = copy(DEFAULT_MULTIVARIATE_OPERATORS)
        logic_operators = [:&&, :||]
        comparison_operators = [:<=, :(==), :>=, :<, :>]
        return new(
            # NODE_CALL_UNIVARIATE
            univariate_operators,
            Dict{Symbol,Int}(
                op => i for (i, op) in enumerate(univariate_operators)
            ),
            length(univariate_operators),
            _UnivariateOperator[],
            # NODE_CALL
            multivariate_operators,
            Dict{Symbol,Int}(
                op => i for (i, op) in enumerate(multivariate_operators)
            ),
            length(multivariate_operators),
            _MultivariateOperator[],
            # NODE_LOGIC
            logic_operators,
            Dict{Symbol,Int}(op => i for (i, op) in enumerate(logic_operators)),
            # NODE_COMPARISON
            comparison_operators,
            Dict{Symbol,Int}(
                op => i for (i, op) in enumerate(comparison_operators)
            ),
        )
    end
end

const _FORWARD_DIFF_METHOD_ERROR_HELPER = raw"""
Common reasons for this include:

 * the function assumes `Float64` will be passed as input, it must work for any
   generic `Real` type.
 * the function allocates temporary storage using `zeros(3)` or similar. This
   defaults to `Float64`, so use `zeros(T, 3)` instead.

As an example, instead of:
```julia
function my_function(x::Float64...)
    y = zeros(length(x))
    for i in 1:length(x)
        y[i] = x[i]^2
    end
    return sum(y)
end
```
use:
```julia
function my_function(x::T...) where {T<:Real}
    y = zeros(T, length(x))
    for i in 1:length(x)
        y[i] = x[i]^2
    end
    return sum(y)
end
```

Review the stacktrace below for more information, but it can often be hard to
understand why and where your function is failing. You can also debug this
outside of JuMP as follows:
```julia
import ForwardDiff

# If the input dimension is 1
x = 1.0
my_function(a) = a^2
ForwardDiff.derivative(my_function, x)

# If the input dimension is more than 1
x = [1.0, 2.0]
my_function(a, b) = a^2 + b^2
ForwardDiff.gradient(x -> my_function(x...), x)
```
"""

_intercept_ForwardDiff_MethodError(err, ::Symbol) = rethrow(err)

function _intercept_ForwardDiff_MethodError(::MethodError, op::Symbol)
    return error(
        "JuMP's autodiff of the user-defined function $(op) failed with a " *
        "MethodError.\n\n$(_FORWARD_DIFF_METHOD_ERROR_HELPER)",
    )
end

function _checked_derivative(f::F, op::Symbol) where {F}
    return function (x)
        try
            return ForwardDiff.derivative(f, x)
        catch err
            _intercept_ForwardDiff_MethodError(err, op)
        end
    end
end

"""
    _validate_register_assumptions(
        f::Function,
        name::Symbol,
        dimension::Integer,
    )

A function that attempts to check if `f` is suitable for registration via
[`register`](@ref) and throws an informative error if it is not.

Because we don't know the domain of `f`, this function may encounter false
negatives. But it should catch the majority of cases in which users supply
non-differentiable functions that rely on `::Float64` assumptions.
"""
function _validate_register_assumptions(
    f::Function,
    name::Symbol,
    dimension::Integer,
)
    # Assumption 1: check that `f` can be called with `Float64` arguments.
    y = 0.0
    try
        if dimension == 1
            y = f(0.0)
        else
            y = f(zeros(dimension))
        end
    catch
        # We hit some other error, perhaps we called a function like log(0).
        # Ignore for now, and hope that a useful error is shown to the user
        # during the solve.
    end
    if !(y isa Real)
        error(
            "Expected return type of `Float64` from the user-defined " *
            "function :$(name), but got `$(typeof(y))`.",
        )
    end
    # Assumption 2: check that `f` can be differentiated using `ForwardDiff`.
    try
        if dimension == 1
            ForwardDiff.derivative(f, 0.0)
        else
            ForwardDiff.gradient(x -> f(x...), zeros(dimension))
        end
    catch err
        if err isa MethodError
            error(
                "Unable to register the function :$name because it does not " *
                "support differentiation via ForwardDiff.\n\n" *
                _FORWARD_DIFF_METHOD_ERROR_HELPER,
            )
        end
        # We hit some other error, perhaps we called a function like log(0).
        # Ignore for now, and hope that a useful error is shown to the user
        # during the solve.
    end
    return
end

function _UnivariateOperator(op::Symbol, f::Function)
    _validate_register_assumptions(f, op, 1)
    f′ = _checked_derivative(f, op)
    f′′ = _checked_derivative(f′, op)
    return _UnivariateOperator(f, f′, f′′)
end

function _UnivariateOperator(op::Symbol, f::Function, f′::Function)
    _validate_register_assumptions(f′, op, 1)
    f′′ = _checked_derivative(f′, op)
    return _UnivariateOperator(f, f′, f′′)
end

function _UnivariateOperator(::Symbol, f::Function, f′::Function, f′′::Function)
    return _UnivariateOperator(f, f′, f′′)
end

function _MultivariateOperator{N}(op::Symbol, f::Function) where {N}
    _validate_register_assumptions(f, op, N)
    g = x -> f(x...)
    ∇f = function (ret, x)
        try
            ForwardDiff.gradient!(ret, g, x)
        catch err
            _intercept_ForwardDiff_MethodError(err, op)
        end
        return
    end
    return _MultivariateOperator{N}(g, ∇f)
end

function _MultivariateOperator{N}(::Symbol, f::Function, ∇f::Function) where {N}
    return _MultivariateOperator{N}(x -> f(x...), (g, x) -> ∇f(g, x...))
end

function _MultivariateOperator{N}(
    ::Symbol,
    f::Function,
    ∇f::Function,
    ∇²f::Function,
) where {N}
    return _MultivariateOperator{N}(
        x -> f(x...),
        (g, x) -> ∇f(g, x...),
        (H, x) -> ∇²f(H, x...),
    )
end

function register_operator(
    registry::OperatorRegistry,
    op::Symbol,
    nargs::Int,
    f::Function...,
)
    if nargs == 1
        if haskey(registry.univariate_operator_to_id, op)
            error("Operator $op is already registered.")
        elseif haskey(registry.multivariate_operator_to_id, op)
            error("Operator $op is already registered.")
        end
        operator = _UnivariateOperator(op, f...)
        push!(registry.univariate_operators, op)
        push!(registry.registered_univariate_operators, operator)
        registry.univariate_operator_to_id[op] =
            length(registry.univariate_operators)
    else
        if haskey(registry.multivariate_operator_to_id, op)
            error("Operator $op is already registered.")
        elseif haskey(registry.univariate_operator_to_id, op)
            error("Operator $op is already registered.")
        end
        operator = _MultivariateOperator{nargs}(op, f...)
        push!(registry.multivariate_operators, op)
        push!(registry.registered_multivariate_operators, operator)
        registry.multivariate_operator_to_id[op] =
            length(registry.multivariate_operators)
    end
    return
end

function _is_registered(registry::OperatorRegistry, op::Symbol, nargs::Int)
    if op in (:<=, :>=, :(==), :<, :>, :&&, :||)
        return true
    end
    if nargs == 1 && haskey(registry.univariate_operator_to_id, op)
        return true
    end
    return haskey(registry.multivariate_operator_to_id, op)
end

function _warn_auto_register(op::Symbol, nargs::Int)
    @warn("""Function $op automatically registered with $nargs arguments.

    Calling the function with a different number of arguments will result in an
    error.

    While you can safely ignore this warning, we recommend that you manually
    register the function as follows:
    ```Julia
    model = Model()
    register(model, :$op, $nargs, $op; autodiff = true)
    ```""")
    return
end

"""
    register_operator_if_needed(
        registry::OperatorRegistry,
        op::Symbol,
        nargs::Int,
        f::Function;
    )

Similar to [`register_operator`](@ref), but this function warns if the function
is not registered, and skips silently if it already is.
"""
function register_operator_if_needed(
    registry::OperatorRegistry,
    op::Symbol,
    nargs::Int,
    f::Function,
)
    if !_is_registered(registry, op, nargs)
        register_operator(registry, op, nargs, f)
        _warn_auto_register(op, nargs)
    end
    return
end

"""
    assert_registered(registry::OperatorRegistry, op::Symbol, nargs::Int)

Throw an error if `op` is not registered in `registry` with `nargs` arguments.
"""
function assert_registered(registry::OperatorRegistry, op::Symbol, nargs::Int)
    if !_is_registered(registry, op, nargs)
        msg = """
        Unrecognized function \"$(op)\" used in nonlinear expression.

        You must register it as a user-defined function before building
        the model. For example, replacing `N` with the appropriate number
        of arguments, do:
        ```julia
        model = Model()
        register(model, :$(op), N, $(op), autodiff=true)
        # ... variables and constraints ...
        ```
        """
        error(msg)
    end
    return
end

"""
    check_return_type(::Type{T}, ret::S) where {T,S}

Overload this method for new types `S` to throw an informative error if a
user-defined function returns the type `S` instead of `T`.
"""
check_return_type(::Type{T}, ret::T) where {T} = nothing

function check_return_type(::Type{T}, ret) where {T}
    return error(
        "Expected return type of $T from a user-defined function, but got " *
        "$(typeof(ret)).",
    )
end

"""
    eval_univariate_function(
        registry::OperatorRegistry,
        op::Symbol,
        x::T,
    ) where {T}

Evaluate the operator `op(x)::T`, where `op` is a univariate function in
`registry`.
"""
function eval_univariate_function(
    registry::OperatorRegistry,
    op::Symbol,
    x::T,
) where {T}
    id = registry.univariate_operator_to_id[op]
    if id <= registry.univariate_user_operator_start
        f, _ = _eval_univariate(id, x)
        return f::T
    end
    offset = id - registry.univariate_user_operator_start
    operator = registry.registered_univariate_operators[offset]
    ret = operator.f(x)
    check_return_type(T, ret)
    return ret::T
end

"""
    eval_univariate_gradient(
        registry::OperatorRegistry,
        op::Symbol,
        x::T,
    ) where {T}

Evaluate the first-derivative of the operator `op(x)::T`, where `op` is a
univariate function in `registry`.
"""
function eval_univariate_gradient(
    registry::OperatorRegistry,
    op::Symbol,
    x::T,
) where {T}
    id = registry.univariate_operator_to_id[op]
    if id <= registry.univariate_user_operator_start
        _, f′ = _eval_univariate(id, x)
        return f′::T
    end
    offset = id - registry.univariate_user_operator_start
    operator = registry.registered_univariate_operators[offset]
    ret = operator.f′(x)
    check_return_type(T, ret)
    return ret::T
end

"""
    eval_univariate_hessian(
        registry::OperatorRegistry,
        op::Symbol,
        x::T,
    ) where {T}

Evaluate the second-derivative of the operator `op(x)::T`, where `op` is a
univariate function in `registry`.
"""
function eval_univariate_hessian(
    registry::OperatorRegistry,
    op::Symbol,
    x::T,
) where {T}
    id = registry.univariate_operator_to_id[op]
    if id <= registry.univariate_user_operator_start
        return _eval_univariate_2nd_deriv(id, x)
    end
    offset = id - registry.univariate_user_operator_start
    operator = registry.registered_univariate_operators[offset]
    ret = operator.f′′(x)
    check_return_type(T, ret)
    return ret::T
end

"""
    eval_multivariate_function(
        registry::OperatorRegistry,
        op::Symbol,
        x::AbstractVector{T},
    ) where {T}

Evaluate the operator `op(x)::T`, where `op` is a multivariate function in
`registry`.
"""
function eval_multivariate_function(
    registry::OperatorRegistry,
    op::Symbol,
    x::AbstractVector{T},
)::T where {T}
    if op == :+
        return sum(x; init = zero(T))
    elseif op == :-
        @assert length(x) == 2
        return x[1] - x[2]
    elseif op == :*
        return prod(x; init = one(T))
    elseif op == :^
        @assert length(x) == 2
        return x[1]^x[2]
    elseif op == :/
        @assert length(x) == 2
        return x[1] / x[2]
    elseif op == :ifelse
        @assert length(x) == 3
        return ifelse(Bool(x[1]), x[2], x[3])
    end
    id = registry.multivariate_operator_to_id[op]
    offset = id - registry.multivariate_user_operator_start
    operator = registry.registered_multivariate_operators[offset]
    @assert length(x) == operator.N
    ret = operator.f(x)
    check_return_type(T, ret)
    return ret::T
end

"""
    eval_multivariate_gradient(
        registry::OperatorRegistry,
        op::Symbol,
        g::AbstractVector{T},
        x::AbstractVector{T},
    ) where {T}

Evaluate the gradient of operator `g .= ∇op(x)`, where `op` is a multivariate
function in `registry`.
"""
function eval_multivariate_gradient(
    registry::OperatorRegistry,
    op::Symbol,
    g::AbstractVector{T},
    x::AbstractVector{T},
) where {T}
    @assert length(g) == length(x)
    if op == :+
        fill!(g, one(T))
    elseif op == :-
        g[1] = one(T)
        g[2] = -one(T)
    elseif op == :*
        # Special case performance optimizations for common cases.
        if length(x) == 1
            g[1] = one(T)
        elseif length(x) == 2
            g[1] = x[2]
            g[2] = x[1]
        else
            total = prod(x)
            if iszero(total)
                for i in 1:length(x)
                    g[i] = prod(x[j] for j in 1:length(x) if i != j)
                end
            else
                for i in 1:length(x)
                    g[i] = total / x[i]
                end
            end
        end
    elseif op == :^
        @assert length(x) == 2
        if x[2] == one(T)
            g[1] = one(T)
        elseif x[2] == T(2)
            g[1] = T(2) * x[1]
        else
            g[1] = x[2] * x[1]^(x[2] - one(T))
        end
        if x[1] > zero(T)
            g[2] = x[1]^x[2] * log(x[1])
        else
            g[2] = T(NaN)
        end
    elseif op == :/
        @assert length(x) == 2
        g[1] = one(T) / x[2]
        g[2] = -x[1] / x[2]^2
    elseif op == :ifelse
        @assert length(x) == 3
        g[1] = zero(T)  # It doesn't matter what this is.
        g[2] = x[1] == one(T)
        g[3] = x[1] == zero(T)
    else
        id = registry.multivariate_operator_to_id[op]
        offset = id - registry.multivariate_user_operator_start
        operator = registry.registered_multivariate_operators[offset]
        @assert length(x) == operator.N
        operator.∇f(g, x)
    end
    return
end

_nan_to_zero(x) = isnan(x) ? 0.0 : x

"""
    eval_multivariate_hessian(
        registry::OperatorRegistry,
        op::Symbol,
        H::AbstractMatrix,
        x::AbstractVector{T},
    ) where {T}

Evaluate the Hessian of operator `∇²op(x)`, where `op` is a multivariate
function in `registry`.

The Hessian is stored in the lower-triangular part of the matrix `H`.

!!! note
    Implementations of the Hessian operators will not fill structural zeros.
    Therefore, before calling this function you should pre-populate the matrix
    `H` with `0`.
"""
function eval_multivariate_hessian(
    registry::OperatorRegistry,
    op::Symbol,
    H,
    x::AbstractVector{T},
) where {T}
    if op in (:+, :-, :ifelse)
        return false
    end
    if op == :*
        # f(x)    = *(x[i] for i in 1:N)
        #
        # ∇fᵢ(x)  = *(x[j] for j in 1:N if i != j)
        #
        # ∇fᵢⱼ(x) = *(x[k] for k in 1:N if i != k & j != k)
        N = length(x)
        if N == 1
            # Hessian is zero
        elseif N == 2
            H[2, 1] = one(T)
        else
            for i in 1:N, j in (i+1):N
                H[j, i] =
                    prod(x[k] for k in 1:N if k != i && k != j; init = one(T))
            end
        end
    elseif op == :^
        # f(x)   = x[1]^x[2]
        #
        # ∇f(x)  = x[2]*x[1]^(x[2]-1)
        #          x[1]^x[2]*log(x[1])
        #
        # ∇²f(x) = x[2]*(x[2]-1)*x[1]^(x[2]-2)
        #          x[1]^(x[2]-1)*(x[2]*log(x[1])+1) x[1]^x[2]*log(x[1])^2
        ln = x[1] > 0 ? log(x[1]) : NaN
        if x[2] == one(T)
            H[2, 1] = _nan_to_zero(ln + one(T))
            H[2, 2] = _nan_to_zero(x[1] * ln^2)
        elseif x[2] == T(2)
            H[1, 1] = T(2)
            H[2, 1] = _nan_to_zero(x[1] * (T(2) * ln + one(T)))
            H[2, 2] = _nan_to_zero(ln^2 * x[1]^2)
        else
            H[1, 1] = x[2] * (x[2] - 1) * x[1]^(x[2] - 2)
            H[2, 1] = _nan_to_zero(x[1]^(x[2] - 1) * (x[2] * ln + 1))
            H[2, 2] = _nan_to_zero(ln^2 * x[1]^x[2])
        end
    elseif op == :/
        # f(x)  = x[1]/x[2]
        #
        # ∇f(x) = 1/x[2]
        #         -x[1]/x[2]^2
        #
        # ∇²(x) = 0.0
        #         -1/x[2]^2 2x[1]/x[2]^3
        d = 1 / x[2]^2
        H[2, 1] = -d
        H[2, 2] = 2 * x[1] * d / x[2]
    else
        id = registry.multivariate_operator_to_id[op]
        offset = id - registry.multivariate_user_operator_start
        operator = registry.registered_multivariate_operators[offset]
        if operator.∇²f === nothing
            error("Hessian is not defined for operator $op")
        end
        @assert length(x) == operator.N
        operator.∇²f(H, x)
    end
    return true
end

"""
    eval_logic_function(
        registry::OperatorRegistry,
        op::Symbol,
        lhs::T,
        rhs::T,
    )::Bool where {T}

Evaluate `(lhs op rhs)::Bool`, where `op` is a logic operator in `registry`.
"""
function eval_logic_function(
    ::OperatorRegistry,
    op::Symbol,
    lhs::T,
    rhs::T,
)::Bool where {T}
    if op == :&&
        return lhs && rhs
    else
        @assert op == :||
        return lhs || rhs
    end
end

"""
    eval_comparison_function(
        registry::OperatorRegistry,
        op::Symbol,
        lhs::T,
        rhs::T,
    )::Bool where {T}

Evaluate `(lhs op rhs)::Bool`, where `op` is a comparison operator in
`registry`.
"""
function eval_comparison_function(
    ::OperatorRegistry,
    op::Symbol,
    lhs::T,
    rhs::T,
)::Bool where {T}
    if op == :<=
        return lhs <= rhs
    elseif op == :>=
        return lhs >= rhs
    elseif op == :(==)
        return lhs == rhs
    elseif op == :<
        return lhs < rhs
    else
        @assert op == :>
        return lhs > rhs
    end
end


text:
export solve!

export LinearSystem

include("Utils.jl")
include("Partial.jl")
include("Curl.jl")
include("CurlCurl.jl")

mutable struct LinearSystem
∇ₛxμⁱ∇ₛx::SparseMatrixCSC{CFloat}
ϵᵣ::Vector{CFloat}
A::SparseMatrixCSC{CFloat}
Q::SparseMatrixCSC{CFloat}
b::Vector{CFloat}
end

function LinearSystem(sim::Simulation)
     ∇ₛxμⁱ∇ₛx_M =  sim(∇ₛxμⁱ∇ₛx)
     ϵᵣ =    sim(convert_to_vector,:ϵᵣ)
     k₀²ϵᵣ = spdiagm(0 => (2pi/sim.λ₀)^2 * ϵᵣ)
     A =     ∇ₛxμⁱ∇ₛx_M - k₀²ϵᵣ
     Q   =   sim(convert_to_diagonal_matrix,:Q)
     src  =   sim(convert_to_vector,:J)
     b =  sim.activate_tfsf ?   (Q*A - A*Q) * src   : src
     LinearSystem(∇ₛxμⁱ∇ₛx_M,ϵᵣ,A,Q,b)
end

function solve! end
function (sim::Simulation)(::typeof(solve!); solver = LU(), args...)
    ls =  LinearSystem(sim)
    linsolve!(sim.E,ls.A,ls.b,solver;args...)
end


text:
module SampleChainsWeighted

using SampleChains, MeasureTheory

using Reexport
@reexport using SampleChains
using NestedTuples
using ElasticArrays
using StructArrays
using ConcreteStructs
using MappedArrays
using Random

export WeightedChain

@concrete struct  WeightedChain{T} <: AbstractChain{T}
    samples     # :: AbstractVector{T}
    logq        # log-density for distribution the sample was drawn from
    logw
    p
    q
    qargs
    rng
end

# TODO: Fix this to use weights
SampleChains.summarize(ch::WeightedChain) = summarize(samples(ch))

function SampleChains.pushsample!(chain::WeightedChain{T}, x::
    push!(samples(chain), transform(gettransform(chain), Q.q))
    push!(logp(chain), Q.ℓq)
    push!(info(chain), tree_stats)
end

function SampleChains.step!(chain::WeightedChain)

end

function SampleChains.initialize!(rng::Random.AbstractRNG, ::Type{WeightedChain}, p, q, qargs)
    qa = q(qargs...)
    x = rand(rng, qa)
    logq = logdensity(qa, x)
    logw = logdensity(p, x) - logq

    chain = WeightedChain(x, logq, logw, p, q, qargs, rng)
end

function SampleChains.initialize!(::Type{WeightedChain}, ℓ, tr, ad_backend=Val(:ForwardDiff))
    rng = Random.GLOBAL_RNG
    return initialize!(rng, WeightedChain, p, q, qargs)
end

function WeightedChain(x::T, logq, logw, p, q, qargs, rng) where {T}
    samples = chainvec(x)
    logq = chainvec(logq)
    logw = chainvec(logw)
    qargs = chainvec(qargs)
    state = chainvec(rng)

    return WeightedChain{T}(samples, logp, logw, p, q, qargs, state)
end

function SampleChains.drawsamples!(chain::WeightedChain, n::Int=1000)
    @cleanbreak for j in 1:n
        Q, tree_stats = step!(chain)
        pushsample!(chain, Q, tree_stats)
    end
    return chain
end

end

max_stars_repo_name: cscherrer/SampleChainsWeighted

max_issues_repo_name: cscherrer/SampleChainsWeighted

max_forks_repo_name: cscherrer/SampleChainsWeighted


text:
using AST
using AtomicLevels

@bsr_exe cfile
@bsr_exe order_c
@bsr_exe w_bsw
@bsr_exe bound_bsw

# The mixing files (.l/.j) are sorted according to energy, but the
# lowest eigenvector does not necessarily have index 1. For each value
# of J, find the indices of the eigenvectors in increasing energy
# order.
function load_mix_file(filename)
    j_pat = r"J =([ 0-9]+)"
    state_pat = filename[end] == 'l' ? r"([A-Z])$" : r"([0-9])([A-Z])([0-9]*)$"
    states = []
    open(filename) do file
        name,Z,nel,ncfg = split(readline(file))[[1,4,7,10]]
        j = 0
        for line in eachline(file)
            if ismatch(j_pat, line)
                j = parse(Int,match(j_pat, line)[1])
            end
            if ismatch(state_pat, line)
                n = parse(Int,match(r"([0-9]+)", line)[1])
                push!(states, (n,j))
            end
        end
    end
    states
end

function bsr_prep_state(orig_wfn::AbstractString,
                        name::AbstractString,
                        filetype::AbstractString,
                        filename::AbstractString,
                        directory::AbstractString;
                        n_eiv::Integer = 1,
                        J2::Integer = 0,
                        cutoff::Real = 0)
    orig_dir = dirname(orig_wfn)
    wfn_filename = basename(orig_wfn)
    cpf(orig_wfn, "$directory/$(name).w")
    for ending in "c$(filetype)"
        cpf("$(orig_dir)/$(name).$(ending)", "$directory/$(name).$(ending)")
    end
    dir_run(directory) do
        mix_file = "$(name).$(filetype)"
        # Read out set index (?) corresponding to the n:th lowest
        # eigenvector
        nn = filter(n -> n[2] == J2, load_mix_file(mix_file))[n_eiv][1]
        cfile(mix_file, nn, J2, "$(filename).c", cutoff)
        w_bsw("$(name).w")
        name != filename && mv("$(name).bsw", "$(filename).bsw", remove_destination = true)
        cpf("$(name).w", "$(filename).w")
    end
end

function bsr_prep_state(state::AST.LevelVertex,
                        filename::AbstractString,
                        directory::AbstractString,
                        cutoff::Real = 0)
    bsr_prep_state(state.wfn, string(state.term), "l",
                   filename, directory; cutoff = cutoff)
end

function add_state!(states::StateDict,
                    term::Term,
                    state::AbstractString)
    if length(get(states, term, [])) == 0
        states[term] = AbstractString[]
    end
    push!(states[term], state)
end

function bsr_prep_states_ls!(directory::AbstractString,
                             wfn::AbstractString,
                             ts::Vector{Term},
                             cfgs::Vector{Config},
                             states::StateDict;
                             filetype::AbstractString = "l",
                             cutoff = 0,
                             all_eivs::Bool = true)
    ts_eiv = all_eivs ? count_common_terms(cfgs, ts) : ones(Int, length(ts))
    for (j,term) in enumerate(ts)
        for i = 1:ts_eiv[j]
            name = "$(string(term))_$i"
            bsr_prep_state(wfn, string(term), filetype,
                           name, directory;
                           n_eiv = i,
                           cutoff = cutoff)
            add_state!(states, term, name)
        end
    end
end

export bsr_prep_state, state_dict, bsr_prep_states_ls!


text:
module Discontinuities

    using ..Quadratures

    export Discontinuity, PathIntegral, PathIntegralLinear, PathIntegralTrigonometric

    export evaluate_l, evaluate_r, derivative_l, derivative_r

    include("discontinuities/path_integral.jl")
    include("discontinuities/path_integral_linear.jl")
    include("discontinuities/path_integral_trigonometric.jl")

    include("discontinuities/discontinuity.jl")

end


text:
# (3) Unemployment and part time emp for econ reason
unempus = get_data(api_fred, "UNEMPLOY"; observation_start = "2008-01-01")
ptempus = get_data(api_fred, "LNS12032194"; observation_start = "2008-01-01")

#Adjust scale to millions of Persons
unempus.data[:value] = unempus.data[:value]./1000;
ptempus.data[:value] = ptempus.data[:value]./1000;

plot_unempus = layer(unempus.data, x = :date, y = :value, Geom.line, Theme(default_color=colorant"red"),
    xintercept = [Date("2020-02-08")], Geom.vline(color = "black", style=[[1mm,1mm]])
    );
plot_ptempus = layer(ptempus.data, x = :date, y = :value, Geom.line, Theme(default_color=colorant"deepskyblue"),
    xintercept = [Date("2008-09-01")], Geom.vline(color = "grey", size=[10.5mm])
    );
plt_labour1 = plot(plot_unempus, plot_ptempus,
        Guide.manual_color_key("Labor market",
                                ["Unemployment level ($(maximum(unempus.data[:date])): $(round(last(unempus.data[:value]), digits = 1)))",
                                 "Employment: Part-time, economic reasons ($(maximum(ptempus.data[:date])): $(round(last(ptempus.data[:value]), digits = 1)))"],
                                ["red", "deepskyblue"],
                                pos = [0.13w,0.38h]
                               ),
         Guide.xlabel("Time"), Guide.ylabel("Millions of Persons"),
        Coord.cartesian(xmin=minimum(weius.data[:date]), xmax=maxvalue=maximum(weius.data[:date])),
        Scale.x_continuous(minvalue=minimum(weius.data[:date]), maxvalue=maximum(weius.data[:date]))
    )
pngout = SVG("plt_labour1.svg", 5inch, 4inch)
draw(pngout, plt_labour1)


text:
function Jav(x::Vector{Float64},Vs::Float64,vs::Float64,ts::Float64,zs::Float64,
    Vlv::Float64,zeta::Float64,Q::Float64,rho::Float64,beta::Float64,W2root::Float64,
    c0::Float64,Kvo::Float64,Kvc::Float64,Ks::Float64,E::Float64,V0::Float64,
    A0::Float64,Aann::Float64,leff::Float64,h::Float64,state::String)

    if length(x) == 3
        J = zeros(3,3);
        J[1,1] = (-2.*ts.*vs.^5./Vs.*(rho./beta).^2.*(((x[1] .- W2root./vs)./8 .+
            c0./vs).^3.*((x[1] .+ W2root./vs)./2) .+ ((x[1] .- W2root./vs)./8 .+
            c0./vs).^4));
        J[1,2] = -1.*ts./h;

        if state == "opening"
            J[2,1] = (0.5.*h.*(1./zs-x[3]).*rho.*vs.^2.*Kvo.*((x[1].-W2root./vs)./8 .+ c0./vs));
            J[2,2] = Vs.*(-h.*(1./zs.-x[3]).*Kvo.*E);
            J[2,3] = (1 .+ h.*Kvo.*(E.*Vs.*(x[2] .- V0./Vs).-beta.*((2.*rho./beta).*vs.^2
                .*((x[1] .- W2root./vs)./8 .+ c0./vs).^2 .- sqrt.(A0))));
        else
            J[2,1] = (0.5.*h*(x[3]).*rho.*vs.^2.*Kvc.*((x[1] .- W2root./vs)./8 .+ c0./vs));
            J[2,2] = Vs.*(-h.*(x[3]).*Kvc.*E);
            J[2,3] = (1 .+ h.*Kvc.*(E.*Vs.*(x[2] .- V0./Vs) .- beta.*((2.*rho./beta).*vs.^2.*
                ((x[1].- W2root./vs)./8 .+ c0./vs).^2 .- sqrt.(A0))));
        end
        if state == "closing" && x[2] > Vlv # closing with backflow
            J[3,2]= (-x[3].*ts./h .+ Vs.*ts./(2.*Aann.*leff.*h.*zs).*
                ((Vlv./Vs .- x[2]) .- abs.(Vlv./Vs .- x[2])) .+ x[3].^2.*Ks.*E.*Vs.*ts.*zs.*
                Aann./(rho.*leff).*(Vlv./Vs .+ V0./Vs .-2.*x[2]));
        elseif state == "opening" || (state == "closing" && x[2] < Vlv) # closing w/ forward flow
            J[3,2] = (-x[3].*ts./h .+ Vs.*ts./(2.*Aann.*leff.*h.*zs).*
                ((-Vlv./Vs .+ x[2]) .- abs.(Vlv./Vs .- x[2])) .+ x[3].^2.*Ks.*E.*Vs.*ts.*zs.*
                Aann./(rho.*leff).*(Vlv./Vs .+ V0./Vs .- 2.*x[2]));
        end
        if state == "opening"
            J[3,3] = ((Vlv./Vs .- x[2]).*ts./h .- Aann.*ts./(rho.*leff.*
                Kvo.*Vs).*((-zeta.*x[3].^2 .+ zs.*x[3].^3)./((1./zs .- x[3]).^2) .+
                (3.*zs.*x[3].^2 .- 2.*zeta.*x[3])./(1./zs .- x[3])) .+
                2.*Ks.*E.*Aann.*ts.*Vs.*zs./(rho.*leff).*(x[2] .- V0./Vs).*
                (Vlv./Vs .- x[2]).*x[3] .- ts./Vs.*Q);
        else
            J[3,3] = ((Vlv./Vs .- x[2]).*ts./h .- Aann.*ts./(rho.*leff.*
                Kvc.*Vs).*(2.*x[3].*zs .- zeta) .+ 2.*Ks.*E.*Aann.*ts.*Vs.*zs./
                (rho.*leff).*(x[2] .- V0./Vs).*(Vlv./Vs .- x[2]).*x[3] .-
                ts./Vs.*Q);
        end
    else # reduced system (constant aortic valve state)
        J = zeros(2,2);
        J[1,1] = (-2.*ts.*vs.^5./Vs.*(rho./beta).^2.*(((x[1] .- W2root./vs)./8 .+
            c0./vs).^3.*((x[1] .+ W2root./vs)./2) .+ ((x[1] .- W2root./vs)./8 .+
            c0./vs).^4));
        J[1,2] = -1.*ts./h;

        if state == "closing" && x[2] > Vlv # closing w/ backflow
            J[2,2] = (-ts./h .+ Vs.*ts./(2.*Aann.*leff.*h.*zs).*((Vlv./Vs .- x[2]) .-
                abs.(Vlv./Vs .- x[2])) .+ Ks.*E.*Vs.*ts.*zs.*Aann./(rho.*leff).*
                (Vlv./Vs .+ V0./Vs .-2.*x[2])); # closing w/ forward flow
        elseif state == "opening" || (state == "closing" && x[2] < Vlv)
            J[2,2] = (-ts./h .+ Vs.*ts./(2.*Aann.*leff.*h.*zs).*((-Vlv./Vs .+ x[2]) .-
                abs.(Vlv./Vs .- x[2])) .+ Ks.*E.*Vs.*ts.*zs.*Aann./(rho.*leff).*
                (Vlv./Vs .+ V0./Vs .- 2.*x[2]));
        end
    end
    return J
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
__precompile__()
module CodecLz4

using Lz4_jll: liblz4
using TranscodingStreams
using TranscodingStreams: TranscodingStream, Memory, Error

export LZ4FrameCompressor, LZ4FrameCompressorStream,
    LZ4FrameDecompressor, LZ4FrameDecompressorStream,
    LZ4FastCompressor, LZ4FastCompressorStream,
    LZ4SafeDecompressor, LZ4SafeDecompressorStream,
    LZ4HCCompressor, LZ4HCCompressorStream,
    BlockSizeID, default_size, max64KB, max256KB, max1MB, max4MB,
    BlockMode, block_linked, block_independent,
    FrameType, normal_frame, skippable_frame,
    lz4_compress, lz4_hc_compress, lz4_decompress

struct LZ4Exception <: Exception
    src::AbstractString
    msg::AbstractString
end

Base.showerror(io::IO, ex::LZ4Exception) = print(io, "$(ex.src): $(ex.msg)")

include("headers/lz4frame.jl")
include("headers/lz4.jl")
include("headers/lz4hc.jl")
include("frame_compression.jl")
include("hc_compression.jl")
include("lz4_compression.jl")
include("simple_compression.jl")

end


text:
using BinaryVectorProbability
using DataFrames

Base.Random.srand(170017)

# N_neurons = 10
# Jtrue = rand(N_neurons, N_neurons); Jtrue = (Jtrue + Jtrue')/2
# P2 = IsingDistribution(Jtrue; comment="true distribution")

#testing if i can recover moments properly
neuron_numbers = 10:5:20
N_trials = 10
df = DataFrame(
    N_neurons=Int[],
    N_samples=Int[],
    trial=Int[],
    Jdiff_MPF=Float64[], Jdiff_LL=Float64[],
    thetadiff_MPF=Float64[], thetadiff_LL=Float64[],
    mudiff_MPF=Float64[], mudiff_LL=Float64[],
    sample_time=Float64[],
    fit_time_MPF=Float64[], fit_time_LL=Float64[])
Ptrue = Vector{IsingDistribution}(length(neuron_numbers))
mutrue = similar(Ptrue, Matrix{Float64})
Phat_MPF = Array{IsingDistribution,2}(length(neuron_numbers), N_trials)
Phat_LL = similar(Phat_MPF)
muhat_MPF = similar(Phat_MPF, Matrix{Float64})
muhat_LL = similar(muhat_MPF)
for (idx, N_neurons) in enumerate(neuron_numbers)
    Jtrue = 10 * randn(N_neurons, N_neurons); Jtrue = (Jtrue + Jtrue')/2
    Ptrue[idx] = IsingDistribution(Jtrue)
    mutrue[idx] = expectation_matrix(Ptrue[idx])
    for N_samples in [100N_neurons, 1000N_neurons, 5000N_neurons]
        for trial = 1:N_trials
            println("*** $N_neurons neurons, $N_samples samples, trial #$trial")
            print("  * Sampling...")
            tic()
            X = random(Ptrue[idx], N_samples)
            sample_time = toq()
            println("done: $sample_time s")

            print("  * Fitting using MPF...")
            tic()
            Phat_MPF[idx, trial] = second_order_model(X, verbose=0, fun="MPF")
            fit_time_MPF = toq()
            println("done: $fit_time_MPF s")

            print("  * Fitting using LL...")
            tic()
            Phat_LL[idx, trial] = second_order_model(X, verbose=0, fun="loglikelihood")
            fit_time_LL = toq()
            println("done:  $fit_time_LL s")

            muhat_MPF[idx, trial] = expectation_matrix(Phat_MPF[idx,trial])
            muhat_LL[idx, trial] = expectation_matrix(Phat_LL[idx,trial])

            Jdiff_MPF = norm(Phat_MPF[idx, trial].J - Ptrue[idx].J) / norm(Ptrue[idx].J)
            thetadiff_MPF = norm(Phat_MPF[idx, trial].theta - Ptrue[idx].theta) / norm(Ptrue[idx].theta)
            mudiff_MPF = norm(muhat_MPF[idx,trial] - mutrue[idx]) / norm(mutrue[idx])
            # maxmu_MPF = maximum(abs.((muhat_MPF[idx,trial] - mutrue[idx]) ./ mutrue[idx]))

            Jdiff_LL = norm(Phat_LL[idx, trial].J - Ptrue[idx].J) / norm(Ptrue[idx].J)
            thetadiff_LL = norm(Phat_LL[idx, trial].theta - Ptrue[idx].theta) / norm(Ptrue[idx].theta)
            mudiff_LL = norm(muhat_LL[idx,trial] - mutrue[idx]) / norm(mutrue[idx])
            # maxmu_LL = maximum(abs.((muhat_LL[idx,trial] - mutrue[idx]) ./ mutrue[idx]))

            # println(">>> Sampling:     $sample_time s")
            # println(">>> Fitting:      $fit_time s")
            println("  * |ΔJ/J| (MPF)  $Jdiff_MPF")
            println("  * |Δθ/θ| (MPF)  $thetadiff_MPF")
            println("  * |Δμ/μ| (MPF)  $mudiff_MPF")
            # println(">>> max Δμ/μ  $maxmu")
            println()
            println("  * |ΔJ/J| (LL)   $Jdiff_LL")
            println("  * |Δθ/θ| (LL)   $thetadiff_LL")
            println("  * |Δμ/μ| (LL)   $mudiff_LL")

            println()
            push!(df, [N_neurons,
                        N_samples,
                        trial,
                        Jdiff_MPF, Jdiff_LL,
                        thetadiff_MPF, thetadiff_LL,
                        mudiff_MPF, mudiff_LL,
                        sample_time,
                        fit_time_MPF, fit_time_LL])
        end
    end
end

# testing if I can fit data keeping H_1 > H_2 > H_n
using BinaryVectorProbability
using Spikes, Optim
Base.Random.srand(170017)

N_neurons = 10
funcs = [t -> 1 + k + (k/4) * sin(t + pi/k) + (k/4) * cos(t/k) for k = 1:N_neurons]
TT = [inhomogeneous_poisson_process(f.(0:0.01:2000), 0.01) for f in funcs]
Sp = SpikeTrains(TT)
X = transpose(raster(Sp, 0.020))

P_1 = first_order_model(X)
P_N = data_model(X)
results = [];
P_2_L = second_order_model(X; verbose=2, show_every=1, allow_f_increases=false, res=results, store_trace=true, J0=zeros(N_neurons,N_neurons))
P_2_f_incr = second_order_model(X; verbose=2, show_every=1, allow_f_increases=true, res=results, J0=zeros(N_neurons,N_neurons))
P_2_Hessian = second_order_model(X; algorithm=Optim.Newton(), verbose=2, show_every=1, allow_f_increases=false, res=results, store_trace=true, J0=zeros(N_neurons,N_neurons))
# J0_new = P_2_MPF.J + Diagonal(P_2_MPF.theta)
# J0 = P_2_MPF.metadata[:J0]
P_2_L = second_order_model(X; verbose=2, allow_f_increases=true, J0=J0)
P_2_different_m = second_order_model(X; verbose=2, J0=J0, algorithm=LBFGS(m=20))
P_2_different_m_increase = second_order_model(X; verbose=2, J0=J0, algorithm=LBFGS(m=20), allow_f_increases=true)
P_2_large_m = second_order_model(X; verbose=2, J0=J0, algorithm=LBFGS(m=40))
P_2_med_m = second_order_model(X; verbose=2, J0=J0, algorithm=LBFGS(m=30))
P_2_thirdtime = second_order_model(X; verbose=2)
P_2_fourthtime = second_order_model(X; verbose=2)
P_N = data_model(X)

H_1 = entropy2(P_1)
H_2_MPF = entropy2(P_2_MPF)
H_2_L = entropy2(P_2_L)
H_2_thirdtime = entropy2(P_2_thirdtime)
H_2_fourthtime = entropy(P_2_fourthtime)
H_N = entropy2(P_N)

mu_1 = expectation_matrix(P_1)
mu_2_MPF = expectation_matrix(P_2_MPF)
mu_2_L = expectation_matrix(P_2_L)
mu_N = expectation_matrix(P_N)

for P in [P_2_MPF, P_2_L, P_2_thirdtime, P_2_fourthtime]
    println("convergence: $(Optim.converged(P.metadata[:opt_res]))\tentropy in bounds: $(H_N < entropy2(P) < H_1)")
end

max_stars_repo_name: sekunder/BinaryVectorProbability

max_issues_repo_name: sekunder/BinaryVectorProbability

max_forks_repo_name: sekunder/BinaryVectorProbability


text:
using PowerModelsReliability
using PowerModels
using InfrastructureModels
using Ipopt
using Mosek
using Juniper
using Cbc
using CPLEX
using Gurobi
using JuMP
using SCS

scs = JuMP.with_optimizer(SCS.Optimizer, max_iters=100000)
ipopt = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)

cplex = JuMP.with_optimizer(CPLEX.Optimizer)
cbc = JuMP.with_optimizer(Cbc.Optimizer)
gurobi = JuMP.with_optimizer(Gurobi.Optimizer)
mosek = JuMP.with_optimizer(Mosek.Optimizer)

juniper = JuMP.with_optimizer(Juniper.Optimizer, nl_solver = ipopt, mip_solver= cbc, time_limit= 7200)

function build_mn_data(base_data)
    mp_data = PowerModels.parse_file(base_data)
    mp_data["load"] = mp_data["sc_load"]
    n_cont = length(mp_data["contingencies"])
    return InfrastructureModels.replicate(mp_data, n_cont)
end

function build_mn_data(base_data_1, base_data_2)
    mp_data_1 = PowerModels.parse_file(base_data_1)
    mp_data_2 = PowerModels.parse_file(base_data_2)

    @assert mp_data_1["per_unit"] == mp_data_2["per_unit"]
    @assert mp_data_1["baseMVA"] == mp_data_2["baseMVA"]

    mn_data = Dict{String,Any}(
        "name" => "$(mp_data_1["name"]) + $(mp_data_2["name"])",
        "multinetwork" => true,
        "per_unit" => mp_data_1["per_unit"],
        "baseMVA" => mp_data_1["baseMVA"],
        "nw" => Dict{String,Any}()
    )

    delete!(mp_data_1, "multinetwork")
    delete!(mp_data_1, "per_unit")
    delete!(mp_data_1, "baseMVA")
    mn_data["nw"]["1"] = mp_data_1

    delete!(mp_data_2, "multinetwork")
    delete!(mp_data_2, "per_unit")
    delete!(mp_data_2, "baseMVA")
    mn_data["nw"]["2"] = mp_data_2

    return mn_data
end

data = build_mn_data("./test/data/case5_scopf.m")
resultAC_rc = run_scunittfopf(data, ACPPowerModel, ipopt; multinetwork=true, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => true, "relax_absolute_value" => true))
data = build_mn_data("./test/data/case5_scopf.m")
resultAC = run_scunittfopf(data, ACPPowerModel, juniper; multinetwork=true, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => false, "relax_absolute_value" => true))
data = build_mn_data("./test/data/case5_scopf.m")
resultDC_rc = run_scunittfopf(data, DCPPowerModel, cplex; multinetwork=true, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => true, "relax_absolute_value" => false))
data = build_mn_data("./test/data/case5_scopf.m")
resultDC = run_scunittfopf(data, DCPPowerModel, cplex; multinetwork=true, setting = Dict("output" => Dict("branch_flows" => true),"relax_continuous" => false, "relax_absolute_value" => false))

# ,"relax_absolute_value" => true
#TODO
# Fixes to make everything compatible
# Extend data model with HVDC contingencies (Branch + Converter)
# Extend load shedding model to possible DC loads
# Start addig contingency constraints for the DC branches and converter models (first check on paper for issues such voltage etc, get inspiration from transformer model)
#

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 91e98d52-29f9-11eb-06af-0fe713025110
using Automata

# ╔═╡ 95d8a1d0-29f9-11eb-0da2-0d030aeeb8b9
using Plots

# ╔═╡ 96ff6e90-29f9-11eb-1c14-075e1e729948
using Images

# ╔═╡ 97f9ad10-29f9-11eb-0be4-93844ea7a8d3
using Flux

# ╔═╡ f58eca9e-29fe-11eb-34ae-d1bd20aea5e0
using Zygote

# ╔═╡ 99f81520-29f9-11eb-388f-27ffa5b411e5
md"
## 1D Automata Function Approximation
"

# ╔═╡ c3bf6cf0-29f9-11eb-28d2-dd6eb9fba2ca
θ = zeros(3)

# ╔═╡ df949cae-29ff-11eb-3371-6dfd3d77f8c7
function perform_step(current_state, map_fn, x)
	uns::Float32 = 0
	i = Int(x[1])
	if i == 1 || i == length(current_state)
		return 0.
	end
	return map_fn(Automata.state_at(current_state, uns, i-1) + Automata.state_at(current_state, uns, i) + Automata.state_at(current_state, uns, i+1))
end

# ╔═╡ f27427b0-29ff-11eb-0c41-21eb6e434cae
Zygote.@adjoint perform_step(x) = perform_step(x), Δ->(Δ,)

# ╔═╡ 27129d2e-2a00-11eb-225e-c9a511df0821
Zygote.@adjoint Base.Iterators.Enumerate(c) = Base.Iterators.Enumerate(c), Δ->(Δ,)

# ╔═╡ 836f62d0-29fa-11eb-17c5-4d74c9a9bb86
function continuous_custom(map_fn, steps::UInt64, initial_state::Vector{Float32})
    current_state = initial_state
    uns::Float32 = 0
    history = Automata.AutomataContinuousHistory()

    for i = 2:steps
		next_state = map((x) -> perform_step(current_state, map_fn, x), enumerate(current_state))
        current_state = next_state

    end

    return current_state
end

# ╔═╡ d1961cc0-29fe-11eb-264b-d7f5355c7c60
Zygote.@adjoint continuous_custom(a, b, c) = continuous_custom(a, b, c), Δ->(Δ, Δ, Δ)

# ╔═╡ 0f3f2d52-29ff-11eb-3722-b15a7bd3449b
Zygote.refresh()

# ╔═╡ d98f8830-29f9-11eb-1f67-a90a4c535ee8
ŷ(x) = continuous_custom(x -> sum(fill(x, length(θ)) .^ collect(0:(length(θ)-1)) .* θ), UInt64(5), Float32[0, 0, 0, 0, x, 0, 0, 0, 0])[5]

# ╔═╡ eb783010-29f9-11eb-18e9-c511e061617e
train_x = Float32.(collect(0:0.01:1))

# ╔═╡ 141d760e-29fa-11eb-2803-17f74479ac10
train_y = train_x .^ 2

# ╔═╡ 18480cee-29fa-11eb-02bb-e578f452e751
loss(x, y) = sum(abs.(y .- ŷ.(x)))

# ╔═╡ 35b55090-29fa-11eb-15c6-fb231004dae2
loss(train_x, train_y)

# ╔═╡ 67b9f9a0-2a00-11eb-1125-f5c97a4e9371
δ = 0.00001

# ╔═╡ dccd9800-2a00-11eb-048d-99648eb16f20
μ = 0.0001

# ╔═╡ 6f7630d0-2a02-11eb-0ecb-4f660a1cb96f
for i = 1:1000
	∇J = gradient(() -> loss(train_x, train_y), params(θ))
	θ .-= μ .* ∇J[θ]
	println(θ)
end

# ╔═╡ ada33e70-2a02-11eb-26a5-cf4310fefb23
@time gradient(() -> loss(train_x, train_y), params(θ))

# ╔═╡ aec82480-2a04-11eb-3925-fb5497d53bec
θ

# ╔═╡ eea21ba0-2a00-11eb-2f02-2330fae2c0e5
begin
	p = plot(train_x, train_y, ylim=(0, 1))
	plot!(p, train_x, ŷ.(train_x))
end

# ╔═╡ 3c2fedc0-2a01-11eb-296e-e128d057977d
ŷ.(train_x)

# ╔═╡ Cell order:
# ╟─99f81520-29f9-11eb-388f-27ffa5b411e5
# ╠═91e98d52-29f9-11eb-06af-0fe713025110
# ╠═95d8a1d0-29f9-11eb-0da2-0d030aeeb8b9
# ╠═96ff6e90-29f9-11eb-1c14-075e1e729948
# ╠═97f9ad10-29f9-11eb-0be4-93844ea7a8d3
# ╠═f58eca9e-29fe-11eb-34ae-d1bd20aea5e0
# ╠═c3bf6cf0-29f9-11eb-28d2-dd6eb9fba2ca
# ╠═df949cae-29ff-11eb-3371-6dfd3d77f8c7
# ╠═f27427b0-29ff-11eb-0c41-21eb6e434cae
# ╠═27129d2e-2a00-11eb-225e-c9a511df0821
# ╠═836f62d0-29fa-11eb-17c5-4d74c9a9bb86
# ╠═d1961cc0-29fe-11eb-264b-d7f5355c7c60
# ╠═0f3f2d52-29ff-11eb-3722-b15a7bd3449b
# ╠═d98f8830-29f9-11eb-1f67-a90a4c535ee8
# ╠═eb783010-29f9-11eb-18e9-c511e061617e
# ╠═141d760e-29fa-11eb-2803-17f74479ac10
# ╠═18480cee-29fa-11eb-02bb-e578f452e751
# ╠═35b55090-29fa-11eb-15c6-fb231004dae2
# ╠═67b9f9a0-2a00-11eb-1125-f5c97a4e9371
# ╠═dccd9800-2a00-11eb-048d-99648eb16f20
# ╠═6f7630d0-2a02-11eb-0ecb-4f660a1cb96f
# ╠═ada33e70-2a02-11eb-26a5-cf4310fefb23
# ╠═aec82480-2a04-11eb-3925-fb5497d53bec
# ╠═eea21ba0-2a00-11eb-2f02-2330fae2c0e5
# ╠═3c2fedc0-2a01-11eb-296e-e128d057977d


text:
"""
**gilljam_parameters**
"""
function gilljam_parameters(parameters, cost, specialistPrefMag, preferenceMethod)
  #preference parameters
  rewireP = preference_parameters(cost, specialistPrefMag, parameters[:A], preferenceMethod)
  #check preferenceMethod
  if preferenceMethod ∈ [:generalist, :specialist]
    rewireP[:preferenceMethod] = preferenceMethod
    rewireP[:specialistPref] = get_specialist_preferences(rewireP,parameters[:A])
  else
    error("Invalid value for preferenceMethod -- must be :generalist or :specialist")
  end
  parameters[:similarity] = rewireP[:similarity]
  parameters[:specialistPrefMag] = rewireP[:specialistPrefMag]
  #parameters[:extinctions] = rewireP[:extinctions]
  parameters[:preferenceMethod] = rewireP[:preferenceMethod]
  parameters[:cost] = rewireP[:cost]
  parameters[:costMat] = rewireP[:costMat]
  parameters[:specialistPref] = rewireP[:specialistPref]
end

"""
Diet similarity rewiring model
"""
function Gilljam(S::Int64, parameters::Dict{Symbol,Any}, biomass::Vector{Float64})
  #create matrix to work on
  preferenceMat = deepcopy(parameters[:A])

  #add extinction to the new matrix
  preferenceMat[:,parameters[:extinctions]] .= 0
  preferenceMat[parameters[:extinctions],:] .= 0

  #where do new links need to be formed?
  newLinks = trues(3,S)
  #Sp that have no prey
  newLinks[1,:] = (sum(preferenceMat, dims = 2) .== 0)
  #Sp that are not producers
  newLinks[2,:] = .!parameters[:is_producer]
  #Sp that are extinction
  for i in parameters[:extinctions]
      newLinks[3,i] = false
  end
  #rows where all conditions are met
  nl = mapslices(all,newLinks, dims = 1)
  newLinks = LinearIndices(nl)[findall(nl)]

  if !isempty(newLinks)
    #code adds new links based on jaccard similariy
    for i in newLinks #for each Sp that needs a new link
      similarities = parameters[:similarity][newLinks][1] #get the similarity ranks
      deleteat!(similarities,findall((in)(parameters[:extinctions]), similarities)) # remove extinctions
      spnoprey = LinearIndices(sum(preferenceMat, dims = 2))[findall(x -> x == 0, sum(preferenceMat, dims = 2))]
      deleteat!(similarities,findall((in)(spnoprey), similarities))#remove Sp without prey
      if !isempty(similarities) #TODO : find a cleaner way to do it
        newLink = sample(findall(x -> x == 1, preferenceMat[similarities[end],:])) #choose link to make
        preferenceMat[i,newLink] = 1 #add  link to the new predation matrix
        parameters[:costMat][i,newLink] = parameters[:cost]
      end
    end
  end
  return(preferenceMat,parameters)
end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
using .StaticArrays: SMatrix, SVector, MMatrix, MVector, SA

eval(load_split_static())
eval(load_genmat_2D_static())
eval(load_genmat_ballinf_static())
eval(load_genmat_hyperrectangle_static())


text:
#    Problem : GROUP A
#    *********
#    A simple quartic function.
#
#    Origonal SIF Source: problem 157 (p. 87) in
#    A.R. Buckley,
#   "Test functions for unconstrained minimization",
#    TR 1989CS-3, Mathematics, statistics and computing centre,
#    Dalhousie University, Halifax (CDN), 1989.
#
#    QUARTC.SIF classification OUR2-AN-V-0
#
#    Number of variables is variable
#
# Daniel Henderson, 08/2021

f = x -> begin
    return sum((x[i] - i)^4 for i in 1:lastindex(x))
end

g! = (g, x) -> begin
	for i in 1:lastindex(x)
		g[i] = 4(x[i] - i)^3
	end
    return g
end

fg! = (g, x) -> begin
	fx = 0.0
	for i in 1:lastindex(x)
		fx  += (x[i] - i)^4
		g[i] = 4(x[i] - i)^3
	end
    return fx, g
end

init = (n::Int=5000) -> begin
    return n, 2.0ones(n)
end

TestSet["QUARTC"] = UncProgram("QUARTC",  f, g!, fg!, init)


text:
using Pluck
using Base.Test

# ------------
#  test array
# ------------

init_array = [
  404,
  123,
  1,
  -1,
  0,
  42
]

# ------------
#  test pluck
# ------------

cur_array = copy(init_array)

for cur_index in 1:( 5 * length(cur_array) )
  @test in(pluck(cur_array), init_array)
end

# -------------
#  test pluck!
# -------------

cur_array = copy(init_array)

actual_values = []

for cur_index in 1:length(cur_array)
  cur_value = pluck!(cur_array)
  @test in(cur_value, init_array)
  push!(actual_values, cur_value)
end

@test_throws BoundsError pluck!(cur_array)

@test sort(init_array) == sort(actual_values)

# ------------
#  test dicts
# ------------

cur_dict = Dict(
  "foo" => +7,
  "bar" => -1
)

@test in(pluck(keys(cur_dict)), keys(cur_dict))

@test in(pluck(values(cur_dict)), values(cur_dict))

@test_throws MethodError pluck!(keys(cur_dict))

@test_throws MethodError pluck!(values(cur_dict))


text:
#///////////////////////////////////////
#// File Name: bic.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/02/25
#// Description: Buffered Input Cells and Related Functions for Benchmarking
#///////////////////////////////////////

import Convex: Variable, norm, minimize, dot, solve!
using GeometricalPredicates
using LazySets
using LinearAlgebra
using Printf
using Polyhedra
using Random
using RobotOS
using SCS
using VoronoiCells

# Buffered Voronoi Cells
# # helper function for sorting Voronoi vertices couterclock-wise
function vertex_direction(generator::IndexablePoint2D, vertex::Point2D)
    x = getx(vertex) - getx(generator)
    y = gety(vertex) - gety(generator)
    return atan(y, x)
end

# # sort Voronoi vertices counterclock-wise
function sort_vertices(generator::IndexablePoint2D, vertices::Vector{Point2D})
    return sort(vertices, by=x->vertex_direction(generator, x))
end

# # scale and shift points (to fit a certain rectangle, expecially [1,2]x[1,2])
function scale_and_shift(x::Real; scale_factor::Real, subtract::Real,
                         offset::Real)
    return x = scale_factor*(x - subtract) + offset
end

# # find extrema of coordinates of all points in the scene
function find_extrema(input_positions::Vector{Vector{Float64}})
    xmin = minimum([p[1] for p in input_positions]);
    xmax = maximum([p[1] for p in input_positions]);
    offset_x = 5. + (xmax - xmin);
    xmin -= offset_x; # additional offset
    xmax += offset_x; # additinal offset
    ymin = minimum([p[2] for p in input_positions]);
    ymax = maximum([p[2] for p in input_positions]);
    offset_y = 5. + (ymax - ymin);
    ymin -= offset_y; # additional offset
    ymax += offset_y; # additional offset
    return xmin, xmax, ymin, ymax
end

# # compute buffered voronoi cell from standard voronoi cell
function retract_cell(generator::IndexablePoint2D,
                      vertices::Vector{Vector{Float64}},
                      radius::Real)
    vpoly = VPolygon(vertices)
    ball = Ball2(zeros(2), radius);
    hpoly = minkowski_difference(vpoly, ball)
    point_vec_new = vertices_list(tovrep(hpoly))
    points = Point2D[]
    for point in point_vec_new
        push!(points, Point2D(point[1], point[2]))
    end
    if isempty(points) # bvc is an empty set. Returning generator.
        p = Point2D(getx(generator), gety(generator))
        return Polygon2D(p)
    else
        return Polygon2D(sort_vertices(generator, points)...)
    end
end

# # extract all agent positions from WorldState
function get_input_positions(w_init::WorldState)
    input_positions = Vector{Vector{Float64}}()
    # ego is always the first element
    push!(input_positions, get_position(w_init.e_state));
    for ap in values(w_init.ap_dict)
        push!(input_positions, ap)
    end
    return input_positions
end

# # main function to compute buffered voronoi cells (idx 1 is for ego robot)
function compute_buffered_voronoi_cells(w_init::WorldState,
                                        radius::Real)
    # compute input positions
    input_positions = get_input_positions(w_init);
    xmin, xmax, ymin, ymax = find_extrema(input_positions)

    # scale and shift
    scaled_positions = [[scale_and_shift(p[1], scale_factor=1. /(xmax - xmin),
                                               subtract=xmin, offset=1.),
                         scale_and_shift(p[2], scale_factor=1. /(ymax - ymin),
                                               subtract=ymin, offset=1.)]
                         for p in input_positions]
    generators = [IndexablePoint2D(p[1], p[2], ii) for (ii, p) in
                  enumerate(scaled_positions)];

    # compute voronoi cells (vertices)
    scaled_voronoi_cells = voronoicells(generators)

    # unshift and unscale
    vertices_vec = Vector{Vector{Vector{Float64}}}(undef, length(generators));
    for ii = 1:length(generators)
        vertices_vec[ii] = [[scale_and_shift(getx(p), scale_factor=(xmax - xmin),
                                                      subtract=1., offset=xmin),
                             scale_and_shift(gety(p), scale_factor=(ymax - ymin),
                                                      subtract=1., offset=ymin)]
                             for p in sort_vertices(generators[ii],
                                                    scaled_voronoi_cells[ii])];
    end

    # compute buffered voronoi cells (polygons)
    buffered_voronoi_cells = Dict{Int64, Polygon2D}()
    for ii = 1:length(generators)
        points = [Point2D(p[1], p[2]) for p in vertices_vec[ii]];
        #buffered_voronoi_cells[ii] = Polygon2D(points...)
        generator = IndexablePoint2D(input_positions[ii][1],
                                     input_positions[ii][2], ii)
        buffered_voronoi_cells[ii] = retract_cell(generator, vertices_vec[ii],
                                                  radius)
    end
    return buffered_voronoi_cells
end

# # robot's dynamics (exact double-integrator, not Euler approximation)
C() = [1. 0. 0. 0.;
       0. 1. 0. 0.]; # state [px, py, vx, vy] to position [px, py]
F(dt) = [1. 0. dt 0.;
         0. 1. 0. dt;
         0. 0. 1. 0.;
         0. 0. 0. 1.]; # dynamics matrix "A"

G(dt) = [0.5*dt^2  0.;
         0.        0.5*dt^2;
         dt        0.;
         0.        dt]; # dynamics matrix "B"

# # main function to compute buffered input cell
function compute_buffered_input_cell(bvc::Polygon2D, # buffered voronoi cell
                                     e_state::RobotState,
                                     dt::Float64)
    x = e_state.x;
    Cm = C();
    Fm = F(dt);
    Gm = G(dt);
    J = Cm*Gm;
    vertices_bvc = [[getx(p), gety(p)] for p in getpoints(bvc)];
    vertices_bic = similar(vertices_bvc);
    for ii = 1:length(vertices_bvc)
        vertices_bic[ii] = J\(vertices_bvc[ii] - Cm*Fm*x);
    end
    return Polygon2D([Point2D(p[1], p[2]) for p in vertices_bic]...)
end

# # project a given 2D point onto buffered input cell (represented as Polygon2D)
function project(point::Vector{Float64}, bic::Polygon2D)
    # QP solver to project target control input to bic
    vertices = [[getx(p), gety(p)] for p in getpoints(bic)];
    vpoly = VPolygon(vertices);
    hpoly = tohrep(vpoly);

    x = Variable(2);
    objective = norm(point - x, 2) # L2 norm between point and p.
    problem = minimize(objective)
    for c in hpoly.constraints
        problem.constraints += dot(c.a, x) <= c.b
    end
    solve!(problem, () -> SCS.Optimizer(verbose=false))
    return Float64[x.value[1], x.value[2]]
end

# LQ Tracker
function lq_tracking(w_init::WorldState,
                     target_trajectory::Trajectory2D,
                     sim_param::SimulationParameter)
    # Reference: Optimal Control by Frank L. Lewis, Table 4.4-1

    # retrieve waypoints from target_trajectory
    prediction_horizon = sim_param.dto*sim_param.prediction_steps;
    ts_times = get_measurement_schedule(w_init, prediction_horizon, sim_param);
    pushfirst!(ts_times, w_init.t); # [t_1, t_2, ..., t_N]
    target_state_array = [[get_position(target_trajectory, t)..., 0.0, 0.0]
                           for t in ts_times]; # [r_1, r_2, ..., r_N, r_N+1]
                                               # final time idx == N+1

    # compute gains
    Gm = G(sim_param.dto);
    Fm = F(sim_param.dto);

    S = Vector{Matrix{Float64}}(undef, sim_param.prediction_steps + 1) # [S_1, S_2, ..., S_N, S_N+1]
    K = Vector{Matrix{Float64}}(undef, sim_param.prediction_steps)     # [K_1, K_2, ..., K_N]
    v = Vector{Vector{Float64}}(undef, sim_param.prediction_steps + 1) # [v_1, v_2, ..., v_N, v_N+1]
    H = Vector{Matrix{Float64}}(undef, sim_param.prediction_steps)     # [H_1, H_2, ..., H_N]

    Q = zeros(4, 4); # state x is [px, py, vx, vy];
    Q[1:2, 1:2] = sim_param.cost_param.Cep;
    R = sim_param.cost_param.Cu;
    S[end] = Q;
    v[end] = Q*target_state_array[end];

    for ii = Iterators.reverse(1:sim_param.prediction_steps)
        K[ii] = (Gm'*S[ii + 1]*Gm + R)\Gm'*S[ii + 1]*Fm;
        S[ii] = Fm'*S[ii + 1]*(Fm - Gm*K[ii]) + Q;
        v[ii] = (Fm - Gm*K[ii])'*v[ii + 1] + Q*target_state_array[ii];
        H[ii] = (Gm'*S[ii + 1]*Gm + R)\Gm'
    end

    # compute optimal control
    u_array = Vector{Vector{Float64}}(undef, sim_param.prediction_steps);    # [u_1, u_2, ..., u_N]
    x_array = Vector{Vector{Float64}}(undef, sim_param.prediction_steps + 1);# [x_1, x_2, ..., x_N, x_N+1]
    x_array[1] = w_init.e_state.x;
    for ii = 1:sim_param.prediction_steps
        u_array[ii] = -K[ii]*x_array[ii] + H[ii]*v[ii + 1];
        x_array[ii + 1] = Fm*x_array[ii] + Gm*u_array[ii]
    end

    return u_array, x_array
end

# BIC Controller
# # bic control parameter & controller
struct BICControlParameter <: Parameter
    eamax::Float64 # Maximum absolute value of accleration
    tcalc::Float64 # Allocated control computation time
    dtr::Float64   # Replanning time interval
    min_dist::Float64 # Minimum allowable distance between ego and ados

    function BICControlParameter(eamax::Float64, tcalc::Float64, dtr::Float64,
                                 min_dist::Float64)
        @assert tcalc <= dtr
        return new(eamax, tcalc, dtr, min_dist)
    end
end

mutable struct BICController
    sim_param::SimulationParameter
    cnt_param::BICControlParameter

    tcalc_actual::Union{Nothing, Float64}
    u_value::Union{Nothing, Vector{Float64}}

    control_update_task::Union{Nothing, Task}

    tcalc_actual_tmp::Union{Nothing, Float64}
    u_value_tmp::Union{Nothing, Vector{Float64}}
end

function BICController(sim_param::SimulationParameter,
                       cnt_param::BICControlParameter)
    return BICController(sim_param, cnt_param, nothing, nothing, nothing,
                         nothing, nothing)
end

# # bic control update
@inline function bic_control_update(w_init::WorldState,
                                    target_trajectory::Trajectory2D,
                                    sim_param::SimulationParameter,
                                    cnt_param::BICControlParameter)
    # LQ Tracker
    tcalc_actual = @elapsed u_array, ~ =
        lq_tracking(w_init, target_trajectory, sim_param);

    # compute buffered input cell
    tcalc_actual +=
        @elapsed bvcs = compute_buffered_voronoi_cells(w_init,
                                                       cnt_param.min_dist/2.);

    tcalc_actual +=
        @elapsed bic_ego = compute_buffered_input_cell(bvcs[1], w_init.e_state,
                                                       sim_param.dto);

    # project u_array[1] to bic_ego
    tcalc_actual +=
        @elapsed u_projected = project(u_array[1], bic_ego);

    # saturate u_projected
    if norm(u_projected) > cnt_param.eamax
        u_projected = u_projected./norm(u_projected).*cnt_param.eamax;
    end

    if tcalc_actual >= cnt_param.tcalc
        # tcalc_actual has exceeded allowable computation time
        time = @sprintf "Time %.2f" round(to_sec(w_init.t), digits=5)
        @warn "$(time) [sec]: BIC computation took $(round(tcalc_actual, digits=3)) [sec], which exceeds the maximum computation time allowed."
    end

    return tcalc_actual, u_projected
end

# # main control function for bic
function control!(controller::BICController,
                  current_time::Time;
                  log::Union{Nothing, Vector{Tuple{Time, String}}}=nothing)

    if !isnothing(controller.control_update_task) &&
        istaskdone(controller.control_update_task)
        if !isnothing(log)
            msg = "New BIC control is available to the controller."
            push!(log, (current_time, msg))
        end
        controller.tcalc_actual = controller.tcalc_actual_tmp;
        controller.u_value = copy(controller.u_value_tmp);
        controller.control_update_task = nothing;
    end

    u = copy(controller.u_value);
    if !isnothing(log)
        msg = "Control: $(u) is applied to the system."
        push!(log, (current_time, msg))
    end
    return u
end

# # control update schedule function for bic
function schedule_control_update!(controller::BICController,
                                  w_init::WorldState,
                                  target_trajectory::Trajectory2D;
                                  log::Union{Nothing, Vector{Tuple{Time, String}}}=nothing)
    if !isnothing(log)
        msg = "New BIC control computation is scheduled."
        push!(log, (w_init.t, msg))
    end

    controller.control_update_task = @task begin
        controller.tcalc_actual_tmp, controller.u_value_tmp =
            bic_control_update(w_init, target_trajectory,
                               controller.sim_param,
                               controller.cnt_param);
    end
    schedule(controller.control_update_task);
end


text:
"""
    HarmonicAngle(i, j, k, th0, cth)

A bond angle between three atoms.
"""
struct HarmonicAngle{T} <: SpecificInteraction
    i::Int
    j::Int
    k::Int
    th0::T
    cth::T
end

# Sometimes domain error occurs for acos if the value is > 1.0 or < -1.0
acosbound(x::Real) = acos(clamp(x, -1, 1))

@inline @inbounds function force(a::HarmonicAngle,
                                  coords,
                                  s::Simulation)
    ba = vector(coords[a.j], coords[a.i], s.box_size)
    bc = vector(coords[a.j], coords[a.k], s.box_size)
    pa = normalize(ba × (ba × bc))
    pc = normalize(-bc × (ba × bc))
    angle_term = -a.cth * (acosbound(dot(ba, bc) / (norm(ba) * norm(bc))) - a.th0)
    fa = (angle_term / norm(ba)) * pa
    fc = (angle_term / norm(bc)) * pc
    fb = -fa - fc
    return [a.i, a.j, a.k], [fa, fb, fc]
end

@inline @inbounds function potential_energy(a::HarmonicAngle,
                                            s::Simulation)
    ba = vector(s.coords[a.j], s.coords[a.i], s.box_size)
    bc = vector(s.coords[a.j], s.coords[a.k], s.box_size)
    pa = normalize(ba × (ba × bc))
    pc = normalize(-bc × (ba × bc))

    a.cth / 2 * (acosbound(dot(ba, bc) / (norm(ba) * norm(bc))) - a.th0)^2
end


text:
using Flux

"""
    RMSPropTF(η, ρ)
Implements the RMSProp algortihm as implemented in tensorflow.
  - Learning Rate (η): Defaults to `0.001`.
  - Rho (ρ): Defaults to `0.9`.
  - Gamma (γ): Defaults to `0.0`.
  - Epsilon (ϵ): Defaults to `1e-6`
## Examples

## References
[RMSProp](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
[Tensorflow RMSProp](https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/RMSprop)
"""
mutable struct RMSPropTF
    eta::Float64
    rho::Float64
    gamma::Float64
    epsilon::Float64
    acc::IdDict
    mom::IdDict
end

RMSPropTF(η = 0.001, ρ = 0.9, γ = 0.0, ϵ = 1e-6) = RMSPropTF(η, ρ, γ, ϵ, IdDict(), IdDict())

function Flux.Optimise.apply!(o::RMSPropTF, x, Δ)
    η, ρ, γ, ϵ = o.eta, o.rho, o.gamma, o.epsilon
    acc = get!(o.acc, x, zero(x))::typeof(x)
    mom = get!(o.mom, x, zero(x))::typeof(x)
    @. acc = ρ * acc + (1 - ρ) * Δ^2
    @. mom = γ * mom + η * Δ/(sqrt(acc + ϵ))
    mom
end

"""
    RMSPropTFCentered(η, ρ)
Implements the Centered version of RMSProp algortihm as implemented in tensorflow.
  - Learning Rate (η): Defaults to `0.001`.
  - Rho (ρ): Defaults to `0.9`.
  - Gamma (γ): Defaults to `0.0`.
  - Epsilon (ϵ): Defaults to `1e-6`
## Examples

## References
[RMSProp](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)
[Tensorflow RMSProp](https://www.tensorflow.org/api_docs/python/tf/keras/optimizers/RMSprop)
"""
mutable struct RMSPropTFCentered
    eta::Float64
    rho::Float64
    gamma::Float64
    epsilon::Float64
    grad::IdDict
    acc::IdDict
    mom::IdDict
end

RMSPropTFCentered(η = 0.001, ρ = 0.9, γ = 0.0, ϵ = 1e-6) =
    RMSPropTFCentered(η, ρ, γ, ϵ, IdDict(), IdDict(), IdDict())

function Flux.Optimise.apply!(o::RMSPropTFCentered, x, Δ)
    η, ρ, γ, ϵ = o.eta, o.rho, o.gamma, o.epsilon
    grad = get!(o.grad, x, zero(x))::typeof(x)
    acc = get!(o.acc, x, zero(x))::typeof(x)
    mom = get!(o.mom, x, zero(x))::typeof(x)
    @. grad = ρ * grad + (1 - ρ) * Δ
    @. acc = ρ * acc + (1 - ρ) * Δ^2
    @. mom = γ * mom + η * Δ/(sqrt(acc - grad^2 + ϵ))
    mom
end


text:
begin
    u = collect(-5:0.5:5)
    ω = ones(21) ./ 21
    prim = [1.0, 0.0, 1.0]
    w = KitML.KitBase.prim_conserve(prim, 5 / 3)
    M = KitML.KitBase.maxwellian(u, prim)
    τ = 0.01
    f = deepcopy(M)
    nn = Chain(Dense(21, 21, tanh), Dense(21, 21))
    p = initial_params(nn)
    nn1 = FastChain(FastDense(21, 21, tanh), FastDense(21, 21))
    p1 = initial_params(nn1)
    fw = zeros(3)
    fh = zeros(21)
    nn2 = Chain(Dense(42, 42, tanh), Dense(42, 42, tanh))
    p2 = initial_params(nn2)
end

KitML.step_ube!(
    fw,
    fh,
    fh,
    w,
    prim,
    f,
    f,
    fw,
    fh,
    fh,
    u,
    ω,
    (2.0, 5 / 3, 0.01, 0.5, 1.0, 0.1, 0.001, zeros(3), zeros(3), nn2, p2);
    mode = :bgk,
)
KitML.step_ube!(
    fw,
    fh,
    fh,
    w,
    prim,
    f,
    f,
    fw,
    fh,
    fh,
    u,
    ω,
    (2.0, 5 / 3, 0.01, 0.5, 1.0, 0.1, 0.001, zeros(3), zeros(3), nn2, p2);
    mode = :shakhov,
)
KitML.step_ube!(
    fw,
    fh,
    fh,
    w,
    prim,
    f,
    f,
    fw,
    fh,
    fh,
    u,
    ω,
    (2.0, 5 / 3, 0.01, 0.5, 1.0, 0.1, 0.001, zeros(3), zeros(3), nn2, p2);
    mode = :nn,
)


text:
module AnalyticsTools

include("plotting/Plotting.jl")

end


text:
using Test
using BinaryProvider, FileIO, Random, Pkg
using MakieGallery
using Makie, AbstractPlotting
using Statistics
database = (get(ENV, "MAKIEGALLERY_MINIMAL", "false") == "true") ? MakieGallery.load_tests() : MakieGallery.load_database()

# THese examples download additional data - don't want to deal with that!
# to_skip = ["WorldClim visualization", "Image on Geometry (Moon)", "Image on Geometry (Earth)"]
# # we directly modify the database, which seems easiest for now
# filter!(entry-> !(entry.title in to_skip), database)

printstyled("Creating ", color = :green, bold = true)
println("recording folders")

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")

println("Diff path  : $tested_diff_path")
println("Record path: $test_record_path")

rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)

rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)

# rerecord = [
#     "Robot Arm"
# ]
# filter!(x-> x.title in rerecord, database)

printstyled("Recording ", color = :green, bold = true)
println("examples")

examples = MakieGallery.record_examples(test_record_path)
if length(examples) != length(database)
    @warn "Not all examples recorded"
end
# MakieGallery.generate_preview(test_record_path, joinpath(homedir(), "Desktop", "index.html"))
# MakieGallery.generate_thumbnails(test_record_path)
# MakieGallery.gallery_from_recordings(test_record_path, joinpath(test_record_path, "index.html"))

printstyled("Running ", color = :green, bold = true)
println("visual regression tests")

MakieGallery.run_comparison(test_record_path, tested_diff_path)


text:
using Documenter, SparseArrays, MaxPlus

# Hack Documenter.jl seems not finding files outside docs/ folder: so copy them
cp(normpath(@__FILE__, "../../README.md"), normpath(@__FILE__, "../src/home.md"); force=true)
cp(normpath(@__FILE__, "../../tutorial/README.md"), normpath(@__FILE__, "../src/tutorials.md"); force=true)

makedocs(
    modules = [MaxPlus],
    format = Documenter.HTML(),
    sitename = "MaxPlus.jl",
#   repo="https://github.com/Lecrapouille/MaxPlus.jl/blob/{commit}{path}#L{line}",
    repo="https://github.com/Lecrapouille/MaxPlus.jl",
    authors="Quentin Quadrat [AT gmail. COM]",
    pages = Any[
        "Home" => "home.md",
        "API" => "API.md",
        "Tutorials" => "tutorials.md",
        "ScicosLab to Julia" => "functions.md",
        "Bibliography" => "bibliography.md",
    ],
)

deploydocs(
    repo = "github.com/Lecrapouille/MaxPlus.jl.git",
    branch = "gh-pages",
    devbranch = "dev",
    devurl = "dev",
)

# Remove copied files
rm(normpath(@__FILE__, "../src/home.md"))
rm(normpath(@__FILE__, "../src/tutorials.md"))


text:
using Base.Test
using DoubleArrayTrie

lines = open(readlines, joinpath(pwd(),"words.txt"))
sort!(lines)

words = map(chomp, lines)
keys = map(w -> [Int(c) for c in w], words)
trie = DATrie(keys, words)

for i = 1:length(words)
  v = get(trie, keys[i], "")
  @test v == words[i]
end


text:
# Modules represent the top-level structure in an LLVM program.

export dispose,
       name, name!,
       triple, triple!,
       datalayout, datalayout!,
       context, inline_asm!,
       set_used!, set_compiler_used!

# forward definition of Module in src/core/value/constant.jl

function Base.unsafe_convert(::Type{API.LLVMModuleRef}, mod::Module)
    # modules can get destroyed, so be sure to check for validity
    mod.ref == C_NULL && throw(UndefRefError())
    mod.ref
end

Base.:(==)(x::Module, y::Module) = (x.ref === y.ref)

# forward declarations
@checked struct DataLayout
    ref::API.LLVMTargetDataRef
end

Module(name::String; ctx::Context) =
    Module(API.LLVMModuleCreateWithNameInContext(name, ctx))

Module(mod::Module) = Module(API.LLVMCloneModule(mod))
Base.copy(mod::Module) = Module(mod)

dispose(mod::Module) = API.LLVMDisposeModule(mod)

function Module(f::Core.Function, args...; kwargs...)
    mod = Module(args...; kwargs...)
    try
        f(mod)
    finally
        dispose(mod)
    end
end

function Base.show(io::IO, mod::Module)
    output = string(mod)
    print(io, output)
end

function name(mod::Module)
    out_len = Ref{Csize_t}()
    ptr = convert(Ptr{UInt8}, API.LLVMGetModuleIdentifier(mod, out_len))
    return unsafe_string(ptr, out_len[])
end
name!(mod::Module, str::String) =
    API.LLVMSetModuleIdentifier(mod, str, Csize_t(length(str)))

triple(mod::Module) = unsafe_string(API.LLVMGetTarget(mod))
triple!(mod::Module, triple) = API.LLVMSetTarget(mod, triple)

datalayout(mod::Module) = DataLayout(API.LLVMGetModuleDataLayout(mod))
datalayout!(mod::Module, layout::String) = API.LLVMSetDataLayout(mod, layout)
datalayout!(mod::Module, layout::DataLayout) =
    API.LLVMSetModuleDataLayout(mod, layout)

inline_asm!(mod::Module, asm::String) =
    API.LLVMSetModuleInlineAsm(mod, asm)

context(mod::Module) = Context(API.LLVMGetModuleContext(mod))

set_used!(mod::Module, values::GlobalVariable...) =
    API.LLVMExtraAppendToUsed(mod, collect(values), length(values))

set_compiler_used!(mod::Module, values::GlobalVariable...) =
    API.LLVMExtraAppendToCompilerUsed(mod, collect(values), length(values))

## named metadata iteration

export metadata, NamedMDNode, operands

# a named metadata note, tying together a name and a MDNode

struct NamedMDNode
    mod::LLVM.Module # not exposed by the API
    ref::API.LLVMNamedMDNodeRef
end

Base.unsafe_convert(::Type{API.LLVMNamedMDNodeRef}, node::NamedMDNode) = node.ref

function name(node::NamedMDNode)
    len = Ref{Csize_t}()
    data = API.LLVMGetNamedMetadataName(node, len)
    unsafe_string(convert(Ptr{Int8}, data), len[])
end

function Base.show(io::IO, mime::MIME"text/plain", node::NamedMDNode)
    print(io, "!$(name(node)) = !{")
    for (i, op) in enumerate(operands(node))
        i > 1 && print(io, ", ")
        show(io, mime, op)
    end
    print(io, "}")
    return io
end

function operands(node::NamedMDNode)
    nops = API.LLVMExtraGetNamedMetadataNumOperands2(node)
    ops = Vector{API.LLVMMetadataRef}(undef, nops)
    if nops > 0
        API.LLVMExtraGetNamedMetadataOperands2(node, ops)
    end
    return [Metadata(op) for op in ops]
end

Base.push!(node::NamedMDNode, val::MDNode) =
    API.LLVMExtraAddNamedMetadataOperand2(node, val)

# module metadata iteration

struct ModuleMetadataIterator <: AbstractDict{String,NamedMDNode}
    mod::Module
end

"""
    metadata(mod)

Fetch the module-level named metadata. This can be inspected using a Dict-like interface.
Mutation is different: There is no `setindex!` method, as named metadata is append-only.
Instead, fetch the named metadata node using `getindex`, and `push!` to it.
"""
metadata(mod::Module) = ModuleMetadataIterator(mod)

function Base.show(io::IO, mime::MIME"text/plain", iter::ModuleMetadataIterator)
    print(io, "ModuleMetadataIterator for module $(name(iter.mod))")
    if !isempty(iter)
        print(io, ":")
        for (key,val) in iter
            print(io, "\n  ")
            show(io, mime, val)
        end
    end
    return io
end

function Base.iterate(iter::ModuleMetadataIterator, state=API.LLVMGetFirstNamedMetadata(iter.mod))
    if state == C_NULL
        nothing
    else
        node = NamedMDNode(iter.mod, state)
        (name(node) => node, API.LLVMGetNextNamedMetadata(state))
    end
end

Base.last(iter::ModuleMetadataIterator) =
    NamedMDNode(iter.mod, API.LLVMGetLastNamedMetadata(iter.mod))

Base.isempty(iter::ModuleMetadataIterator) =
    API.LLVMGetLastNamedMetadata(iter.mod) == C_NULL

Base.IteratorSize(::ModuleMetadataIterator) = Base.SizeUnknown()

function Base.haskey(iter::ModuleMetadataIterator, name::String)
    return API.LLVMGetNamedMetadata(iter.mod, name, length(name)) != C_NULL
end

function Base.getindex(iter::ModuleMetadataIterator, name::String)
    ref = API.LLVMGetOrInsertNamedMetadata(iter.mod, name, length(name))
    @assert ref != C_NULL
    node = NamedMDNode(iter.mod, ref)
    return node
end

## global variable iteration

export globals

struct ModuleGlobalSet
    mod::Module
end

globals(mod::Module) = ModuleGlobalSet(mod)

Base.eltype(::ModuleGlobalSet) = GlobalVariable

function Base.iterate(iter::ModuleGlobalSet, state=API.LLVMGetFirstGlobal(iter.mod))
    state == C_NULL ? nothing : (GlobalVariable(state), API.LLVMGetNextGlobal(state))
end

Base.last(iter::ModuleGlobalSet) =
    GlobalVariable(API.LLVMGetLastGlobal(iter.mod))

Base.isempty(iter::ModuleGlobalSet) =
    API.LLVMGetLastGlobal(iter.mod) == C_NULL

Base.IteratorSize(::ModuleGlobalSet) = Base.SizeUnknown()

# partial associative interface

function Base.haskey(iter::ModuleGlobalSet, name::String)
    return API.LLVMGetNamedGlobal(iter.mod, name) != C_NULL
end

function Base.getindex(iter::ModuleGlobalSet, name::String)
    objref = API.LLVMGetNamedGlobal(iter.mod, name)
    objref == C_NULL && throw(KeyError(name))
    return GlobalVariable(objref)
end

## function iteration

export functions

struct ModuleFunctionSet
    mod::Module
end

functions(mod::Module) = ModuleFunctionSet(mod)

Base.eltype(::ModuleFunctionSet) = Function

function Base.iterate(iter::ModuleFunctionSet, state=API.LLVMGetFirstFunction(iter.mod))
    state == C_NULL ? nothing : (Function(state), API.LLVMGetNextFunction(state))
end

Base.last(iter::ModuleFunctionSet) =
    Function(API.LLVMGetLastFunction(iter.mod))

Base.isempty(iter::ModuleFunctionSet) =
    API.LLVMGetLastFunction(iter.mod) == C_NULL

Base.IteratorSize(::ModuleFunctionSet) = Base.SizeUnknown()

# partial associative interface

function Base.haskey(iter::ModuleFunctionSet, name::String)
    return API.LLVMGetNamedFunction(iter.mod, name) != C_NULL
end

function Base.getindex(iter::ModuleFunctionSet, name::String)
    objref = API.LLVMGetNamedFunction(iter.mod, name)
    objref == C_NULL && throw(KeyError(name))
    return Function(objref)
end

## module flag iteration
# TODO: doesn't actually iterate, since we can't list the available keys

export flags

struct ModuleFlagDict <: AbstractDict{String,Metadata}
    mod::Module
end

flags(mod::Module) = ModuleFlagDict(mod)

Base.haskey(iter::ModuleFlagDict, name::String) =
    API.LLVMGetModuleFlag(iter.mod, name, length(name)) != C_NULL

function Base.getindex(iter::ModuleFlagDict, name::String)
    objref = API.LLVMGetModuleFlag(iter.mod, name, length(name))
    objref == C_NULL && throw(KeyError(name))
    return Metadata(objref)
end

function Base.setindex!(iter::ModuleFlagDict, val::Metadata,
                        (name, behavior)::Tuple{String, API.LLVMModuleFlagBehavior})
    API.LLVMAddModuleFlag(iter.mod, behavior, name, length(name), val)
end