text:
function _add_system_balance_slacks!(
    psi_container::PSIContainer,
    slack_name::String,
    expression::Symbol,
    single_first_axes::Bool = false,
)
    time_steps = model_time_steps(psi_container)
    expression_array = get_expression(psi_container, expression)
    var_name_up = make_variable_name(slack_name, SLACK_UP)
    var_name_dn = make_variable_name(slack_name, SLACK_DN)
    single_first_axes && (first_index = [axes(expression_array)[1][1]])
    !single_first_axes && (first_index = axes(expression_array)[1])
    variable_up = add_var_container!(psi_container, var_name_up, first_index, time_steps)
    variable_dn = add_var_container!(psi_container, var_name_dn, first_index, time_steps)
    for ix in first_index, jx in time_steps
        variable_up[ix, jx] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "$(var_name_up)_{$(ix), $(jx)}",
            lower_bound = 0.0
        )
        variable_dn[ix, jx] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "$(var_name_dn)_{$(ix), $(jx)}",
            lower_bound = 0.0
        )
        add_to_expression!(expression_array, ix, jx, variable_up[ix, jx], 1.0)
        add_to_expression!(expression_array, ix, jx, variable_dn[ix, jx], -1.0)
        JuMP.add_to_expression!(
            psi_container.cost_function,
            (variable_dn[ix, jx] + variable_up[ix, jx]) * BALANCE_SLACK_COST,
        )
    end
    return
end

function add_slacks!(psi_container::PSIContainer, ::Type{CopperPlatePowerModel})
    _add_system_balance_slacks!(psi_container, ACTIVE_POWER, :nodal_balance_active, true)
    return
end

function add_slacks!(
    psi_container::PSIContainer,
    ::Type{T},
) where {T <: PM.AbstractActivePowerModel}
    _add_system_balance_slacks!(psi_container, ACTIVE_POWER, :nodal_balance_active)
    return
end

function add_slacks!(
    psi_container::PSIContainer,
    ::Type{T},
) where {T <: PM.AbstractPowerModel}
    _add_system_balance_slacks!(psi_container, ACTIVE_POWER, :nodal_balance_active)
    _add_system_balance_slacks!(psi_container, REACTIVE_POWER, :nodal_balance_reactive)
    return
end


text:
using Test
using ProjectIO
using DataFrames

data_dir = joinpath(dirname(dirname(pathof(ProjectIO))), "test/data")
spectrum_file = data_dir * "/liquid_crystal_in_etalon.csv"
angles_dir = data_dir * "/angle_resolved"

println(pwd())

@testset "read JASCO FTIR csv file" begin
    df = ProjectIO.read_spectrum(spectrum_file)
    @test df[:, 1][1] == 999.9101
    @test round(df[:, 1][end], sigdigits=8) == 6000.4248
    @test round(df[:, 2][1], sigdigits=3) ≈ 0.00673
    @test size(df) == (10373, 2)
    @test names(df)[1] == "Wavenumber"
    @test names(df)[2] == "Transmittance"
end

@testset "read angle-resolved data" begin
    datalist = ProjectIO.read_angle_data_from_dir(angles_dir)
    @test length(datalist) == 6
    @test [typeof(i[1]) == Int for i in datalist] == ones(length(datalist))
    @test typeof(datalist[1][2]) == DataFrame
    for d in datalist
        @test length(d) == 2
    end
end

ext: jl

max_stars_repo_name: garrekstemo/projectIO

max_issues_repo_name: garrekstemo/projectIO

max_forks_repo_name: garrekstemo/projectIO


text:
type PAWarpParams
    shape::Matrix{Float64}             # shape to warp to
    dst_size::@compat(Tuple{Int, Int}) # destination size
    trigs::Matrix{Int}                 # triangulation of a shape
    warp_map::Matrix{Int}              # map from pixel in destination to
                                       #   corresponding triangle
    alpha_coords::Matrix{Float64}      # 1st barycentric coordinate w.r.t.
                                       #   corrsponding triangle in warp_map
    beta_coords::Matrix{Float64}       # 2nd barycentric coordinate w.r.t.
                                       #   corrsponding triangle in warp_map
end

function pa_warp_params(shape::Matrix{Float64}, trigs::Matrix{Int},
                        dst_size::@compat(Tuple{Int, Int}))
    h, w = dst_size
    warp_map = zeros(Int, h, w)
    alpha_coords = zeros(h, w)
    beta_coords = zeros(h, w)
    for j=1:w
        for i=1:h
            for k=1:size(trigs, 1)
                t = trigs[k, :]
                i1 = shape[t[1], 1]
                j1 = shape[t[1], 2]
                i2 = shape[t[2], 1]
                j2 = shape[t[2], 2]
                i3 = shape[t[3], 1]
                j3 = shape[t[3], 2]

                den = (i2 - i1) * (j3 - j1) - (j2 - j1) * (i3 - i1)
                alpha = ((i - i1) * (j3 - j1) - (j - j1) * (i3 - i1)) / den
                beta = ((j - j1) * (i2 - i1) - (i - i1) * (j2 - j1)) / den

                if alpha >= 0 && beta >= 0 && (alpha + beta) <= 1
                    # found the triangle, save data to the bitmaps and break
                    warp_map[i, j] = k
                    alpha_coords[i, j] = alpha
                    beta_coords[i,j] = beta
                    break
                end
            end
        end
    end
    return PAWarpParams(copy(shape), dst_size, trigs,
                        warp_map, alpha_coords, beta_coords)
end

function pa_warp{T,N}(params::PAWarpParams, src_img::Array{T,N},
                    src_shape::Matrix{Float64};
                    interp=:bilinear)
    nc = size(src_img, 3)
    h, w = params.dst_size
    warp_map = params.warp_map
    alpha_coords, beta_coords = params.alpha_coords, params.beta_coords
    dst_img = zeros(T, h, w, nc)
    for j=1:w
        for i=1:h
            t = warp_map[i, j]
            # if t <= 0, pixel is out of destination shape
            if t > 0
                # index of first vertex of the triangle
                v1 = params.trigs[t, 1]
                i1 = src_shape[v1, 1]
                j1 = src_shape[v1, 2]

                v2 = params.trigs[t, 2]
                i2 = src_shape[v2, 1]
                j2 = src_shape[v2, 2]

                v3 = params.trigs[t, 3]
                i3 = src_shape[v3, 1]
                j3 = src_shape[v3, 2]

                wi = (i1 +
                      alpha_coords[i, j] * (i2 - i1) +
                      beta_coords[i, j] * (i3 - i1))
                wj = (j1 +
                      alpha_coords[i, j] * (j2 - j1) +
                      beta_coords[i, j] * (j3 - j1))

                if wi < 1 || wi > size(src_img, 1) ||
                    wj < 1 || wj > size(src_img, 2)
                    # throw(@compat(BoundsError("Warp pixel is out of bounds: " *
                    #                   "wi=$wi, wj=$wj")))
                    throw(BoundsError())
                end

                lli = convert(Int, floor(wi))
                llj = convert(Int, floor(wj))
                uri = lli + 1
                urj = llj + 1

                f0 = (uri - wi) * (urj - wj)
                f1 = (wi - lli) * (urj - wj)
                f2 = (uri - wi) * (wj - llj)
                f3 = (wi - lli) * (wj - llj)

                for c=1:nc
                    if interp == :bilinear
                        interpolated = (src_img[lli, llj, c] * f0 +
                                        src_img[uri, llj, c] * f1 +
                                        src_img[lli, urj, c] * f2 +
                                        src_img[uri, urj, c] * f3)
                        dst_img[i, j, c] = interpolated
                    elseif interp == :nearest
                        dst_img[i, j, c] = src_img[convert(Int, round(wi)),
                                                   convert(Int, round(wj)), c]
                    else
                        error("Unknown interpolation type: $interp")
                    end
                end
            end
        end
    end
    # for 2D images, fix dimensions
    if N == 2
        dst_img = squeeze(dst_img, (3,))
    end
    return dst_img
end


text:
# # Replication of Sims (2010)

# This example replicates [Sims (2010)](http://sims.princeton.edu/yftp/RIMP/handbookChapterRI2.pdf) from the Handbook of Monetary Economics using the [DRIPs](https://github.com/afrouzi/DRIPs.jl) package.

# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/afrouzi/DRIPs.jl/binder?filepath=examples) to run and modify the following code (no software is needed on the local machine).

# See [Afrouzi and Yang (2020)](http://www.afrouzi.com/dynamic_inattention.pdf) for background on the theory.

# ## Contents
# * [Setup](@ref sims2010_setup)
# * [Initialization](@ref sims2010_param)
# * [Solution](@ref sims2010_solution)
#     * [Benchmark Parameterization](@ref sims2010_benchmark)
#     * [Lower Cost of Attention: $\omega = 0.1$](@ref sims2010_lowomega)
#     * [Other Discount Factors: $\beta \in \{0,1\}$](@ref sims2010_betas)
# * [Impulse Response Functions](@ref sims2010_figures)
#     * [Benchmark Parameterization](@ref sims2010_fig_benchmark)
#     * [Lower Cost of Attention: $\omega = 0.1$](@ref sims2010_fig_lowomega)
#     * [Other Discount Factors: $\beta \in \{0,1\}$](@ref sims2010_fig_betas)
# * [Extensions](@ref sims2010_extensions)
#     * [Transition Dynamics of Attention](@ref sims2010_trip)
#     * [Impulse Response Functions with Information Treatment](@ref sims2010_trip_irfs)

# ## [Setup](@id sims2010_setup)

# The problem in [Sims (2011)](http://sims.princeton.edu/yftp/RIMP/handbookChapterRI2.pdf), as it appears on page 21, with slight change of notation,
# ```math
# \begin{aligned}
#             & \min_{\{\Sigma_{t|t}\succeq 0\}_{t\geq 0}} \mathbb{E}_0\left[\sum_{t=0}^\infty
#   \beta^t \left(tr(\Sigma_{t|t}\mathbf{H}\mathbf{H}')+\omega\log\left(\frac{|\Sigma_{t|t-1}|}{|\Sigma_{t|t}|}\right)\right)\right] \\
#   s.t.\quad &
#   \Sigma_{t+1|t}=\mathbf{A}\Sigma_{t|t}\mathbf{A}'+\mathbf{Q}\mathbf{Q}'\\
#             & \Sigma_{t|t-1}-\Sigma_{t|t} \text{ positive semi-definite}
# \end{aligned}
# ```

# where
# ```math
# \begin{aligned}
#   \mathbf{H} = \left[\begin{array}{c} 1 \\ 1\end{array}\right],
#     \quad
#     \mathbf{A} = \left[\begin{array}{cc}
#                           0.95 & 0\\
#                           0 & 0.4\\
#                      \end{array}\right],
#   \quad
#   \mathbf{Q} = \left[\begin{array}{cc}
#                           \sqrt{0.0975} & 0\\
#                           0           & \sqrt{0.86}\\
#                      \end{array}\right]
# \end{aligned}
# ```
# We have renamed the parameters so that the problem directly maps to a D.R.I.P. Otherwise, the problem is the same.

# ## [Initialization](@id sims2010_param)
# Include the package:
using DRIPs;

# Set parameters:

β = 0.9;
ω = 1.0;
A = [0.95 0.0; 0.0 0.4];
Q = [√0.0975 0.0; 0.0 √0.86];
H = [1.0; 1.0];

# ## [Solution and Performance](@id sims2010_solution)
# ### [Benchmark Parameterization](@id sims2010_benchmark)
# Solve and display the optimal posterior covariance matrix:

sol_bp = Drip(ω,β,A,Q,H);
sol_bp.ss.Σ_p

# Performance for random values of $\omega\in [0,2]$:

using BenchmarkTools;
@benchmark Drip(ω,β,A,Q,H) setup = (ω = 2*rand())

# Performance for random values of $\beta\in[0,1]$:

@benchmark Drip(ω,β,A,Q,H) setup = (β = rand())

# ### [Lower Cost of Attention: $\omega = 0.1$](@id sims2010_lowomega)
# Solve and display the optimal posterior covariance matrix:

sol_lω = Drip(0.1,β,A,Q,H);
sol_lω.ss.Σ_p

# ### [Different Discount Factors: $\beta \in \{0,1\}$](@id sims2010_betas)
# Solve the model for $\beta=0$ and $\beta=1$ to compare with the benchmark value of $\beta=0.9$:

# ``\beta = 0``
sol_lβ = Drip(ω,0,A,Q,H);
sol_lβ.ss.Σ_p

# ``\beta = 1``:
sol_hβ = Drip(ω,1,A,Q,H);
sol_hβ.ss.Σ_p

# ## [Impulse Response Functions](@id sims2010_figures)
# ### [Benchmark Parameterization](@id sims2010_fig_benchmark)
# Get the IRFs:

T = 25;
irfs_bp = irfs(sol_bp,T = T);

# Plot IRFs:

using Plots, LaTeXStrings; pyplot();
p1 = plot(1:T, [irfs_bp.x[1,1,:], irfs_bp.a[1,1,:]],
    title             = L"IRFs to Slow-Moving Shock ($\rho = 0.95$)",
    label             = ["Shock" "Price"],
    color             = [:darkgray :black],
    marker            = [:circle :square],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p2 = plot(1:T, [irfs_bp.x[2,2,:], irfs_bp.a[1,2,:]],
    title             = L"IRFs to Fast-Moving Shock ($\rho = 0.4$)",
    label             = ["Shock" "Price"],
    color             = [:darkgray :black],
    marker            = [:circle :square],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p = plot(p1,p2,
    layout     = (2,1),
    xlabel     = "Time",
    lw         = 2,
    xticks     = (1:2:T),
    xlim       = (0,T+1),
    fontfamily = "serif",
    legend     = :topright,
    legendfont = font(12),
    tickfont   = font(12),
    size       = (900,550),
    framestyle = :box)

# ### [Lower Cost of Attention: $\omega=0.1$](@id sims2010_fig_lowomega)
# Get the IRFs:

T = 25; #length of IRFs
irfs_lω = irfs(sol_lω,T = T);

# Plot IRFs:

p1 = plot(1:T, [irfs_lω.x[1,1,:], irfs_lω.a[1,1,:]],
    title             = L"IRFs to Slow-Moving Shock ($\rho = 0.95$)",
    label             = ["Shock" "Price"],
    color             = [:darkgray :black],
    marker            = [:circle :square],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p2 = plot(1:T, [irfs_lω.x[2,2,:], irfs_lω.a[1,2,:]],
    title             = L"IRFs to Fast-Moving Shock ($\rho = 0.4$)",
    label             = ["Shock" "Price"],
    color             = [:darkgray :black],
    marker            = [:circle :square],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p = plot(p1,p2,
    layout     = (2,1),
    xlabel     = "Time",
    lw         = 2,
    xticks     = (1:2:T),
    xlim       = (0,T+1),
    fontfamily = "serif",
    legend     = :topright,
    legendfont = font(12),
    tickfont   = font(12),
    size       = (900,550),
    framestyle = :box)

# ### [Other Discount Factors: $\beta\in\{0,1\}$](@id sims2010_fig_betas)
# Get the IRFs:

T = 25; #length of IRFs
irfs_lβ = irfs(sol_lβ,T = T);
irfs_hβ = irfs(sol_hβ,T = T);

# Plot IRFs:

p1 = plot(1:T, [irfs_bp.x[1,1,:],irfs_hβ.a[1,1,:], irfs_lβ.a[1,1,:]],
    title             = L"IRFs to Slow-Moving Shock ($\rho = 0.95$)",
    label             = ["Shock" L"Price ($\beta=1$)" L"Price ($\beta=0$)"],
    color             = [:darkgray :black :gray50],
    marker            = [:circle :square :utriangle],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p2 = plot(1:T, [irfs_bp.x[2,2,:],irfs_hβ.a[1,2,:], irfs_lβ.a[1,2,:]],
    title             = L"IRFs to Fast-Moving Shock ($\rho = 0.4$)",
    label             = ["Shock" L"Priceblack ($\beta=1$)" L"Price ($\beta=0$)"],
    color             = [:darkgray :black :gray50],
    marker            = [:circle :square :utriangle],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p = plot(p1,p2,
    layout     = (2,1),
    xlabel     = "Time",
    lw         = 2,
    xticks     = (1:2:T),
    xlim       = (0,T+1),
    fontfamily = "serif",
    legend     = :topright,
    legendfont = font(12),
    tickfont   = font(12),
    size       = (900,550),
    framestyle = :box)

# ## [Extensions](@id sims2010_extensions)

# ### [Transition Dynamics of Attention](@id sims2010_trip)

# In this section, we solve for the transition dynamics of the optimal posterior covariance matrix starting from an initial prior that is different from the steady state prior.

# For instance let us consider a case where the firm is at the steady state of the rational inattention problem at time 0, with prior covariance matrix $\bar{\Sigma}_{-1}$, and it receives a one time treatment with a perfectly informative signal about its optimal price:

# $$s_0 = \mathbf{H}'\vec{x}_0$$

# #### Solve for the transition dynamics
# The function `Trip` solves for the transition dynamics automatically given the initial signal. Start by initializing the initial signal:

s0 = DRIPs.Signal(H,0.0);

# Solve for the transition dynamics given $s_0$:

Tss     = 15; # guess for time until convergence
bp_trip = Trip(sol_bp, s0; T = Tss);

# Performance for solving the transition dynamics for a random signal:

@benchmark Trip(sol_bp, S; T = 30) setup = (S = DRIPs.Signal(rand(2),0.0))

# #### Plot Transition Path of Eigenvalues

# Plot the marginal values of information. In this problem the state is two dimensional. At any time, for every orthogonalized dimension, the agent weighs the **marginal value** of acquiring information in that dimension against the **marginal cost** of attention which is the parameter $\omega$.**The number of signals that the agent acquires at any time is the number of marginal values that are larger than $\omega$.**

p = plot(0:Tss-1,[bp_trip.Ds[1,1:Tss],bp_trip.Ds[2,1:Tss],bp_trip.p.ω*ones(Tss,1)],
    label             = ["Low marginal value dim." "High marginal value dim." "Marginal cost of attention"],
    size              = (900,275),
    title             = "Marginal Value of Information",
    xlabel            = "Time",
    color             = [:darkgray :black :black],
    line              = [:solid :solid :dash],
    marker            = [:circle :square :none],
    markercolor       = false,
    markerstrokecolor = :match,
    markersize        = 6,
    xlim              = (-1,Tss),
    xticks            = 0:2:Tss-1,
    legend            = :outertopright,
    fontfamily        = "serif",
    framestyle        = :box)

# ### [Impulse Response Functions with Information Treatment](@id sims2010_trip_irfs)
# Get the IRFs in the transition path after treatment:

T = 30;

tirfs_bp = irfs(sol_bp,s0,T = T); # irfs with treatment
irfs_bp  = irfs(sol_bp,T = T);    # irfs in the Ss (without treatment)

# Plot IRFs:

p1 = plot(1:T, [irfs_bp.x[1,1,:], tirfs_bp.a[1,1,:], irfs_bp.a[1,1,:]],
    title             = L"IRFs to Slow-Moving Shock ($\rho = 0.95$)",
    label             = ["Shock" "Price (w/ treatment)" "Price (w/o treatment)"],
    color             = [:darkgray :black :gray80],
    marker            = [:circle :square :utriangle],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p2 = plot(1:T, [tirfs_bp.x[2,2,:], tirfs_bp.a[1,2,:], irfs_bp.a[1,2,:]],
    title             = L"IRFs to Fast-Moving Shock ($\rho = 0.4$)",
    label             = ["Shock" "Price (w/ treatment)" "Price (w/o treatment)"],
    color             = [:darkgray :black :gray80],
    marker            = [:circle :square :utriangle],
    markerstrokecolor = :match,
    markercolor       = false,
    markersize        = 6)
p = plot(p1,p2,
    layout     = (2,1),
    xlabel     = "Time",
    lw         = 2,
    xticks     = (1:2:T),
    xlim       = (0,T+1),
    fontfamily = "serif",
    legend     = :topright,
    legendfont = font(12),
    tickfont   = font(12),
    size       = (900,550),
    framestyle = :box)


text:
__precompile__()

module RBT

export
    setdb,
    getdb,
    @q_str,
    query,
    @query

import Base:
    show, showcompact, showall, call, convert, get, isnull, isempty, length, eltype,
    start, next, done, endof, getindex, max, typemin, foreach, eachindex,
    ==, +, -, *, /, >>, ^, .==, .!=, .>, .>=, .<, .<=, ~, &, |

include("immdict.jl")
include("monetary.jl")
include("kind.jl")
include("domain.jl")
include("syntax.jl")
include("pipe.jl")
include("database.jl")
include("scope.jl")
include("binding.jl")
include("compile.jl")

has_dataframes =
    try
        using DataFrames
        true
    catch ArgumentError
        false
    end

if has_dataframes
    include("df.jl")
end

global DB = ToyDatabase(Schema(), Instance(Dict(), Dict()))

function setdb(db)
    global DB = db
end

function getdb()
    return DB
end

function query(state, expr; params...)
    params_signature = ([(name, typeof(param)) for (name, param) in params]...)
    return execute(prepare(state, expr; params_signature...); params...)
end

function query(expr; params...)
    return query(DB, expr; params...)
end

macro query(args...)
    L = endof(args)
    while L > 0 && isa(args[L], Expr) && args[L].head == :kw
        L = L-1
    end
    @assert 1 <= L <= 2
    if L == 2
        db, expr = args[1], args[2]
    else
        db, expr = DB, args[1]
    end
    expr = syntax(expr)
    params = args[L+1:end]
    return quote
        query($(esc(db)), $expr; $(params...))
    end
end

macro q_str(str)
    query(str)
end

function prepare(expr; params_signature...)
    return prepare(DB, expr; params_signature...)
end

macro prepare(args...)
    L = endof(args)
    while L > 0 && isa(args[L], Expr) && args[L].head == :kw
        L = L-1
    end
    @assert 1 <= L <= 2
    if L == 2
        db, expr = args[1], args[2]
    else
        db, expr = DB, args[1]
    end
    expr = syntax(expr)
    params = args[L+1:end]
    return quote
        prepare($(esc(db)), $expr; $(params...))
    end
end

end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using ObjectDetector
using Test, PrettyTables
using FileIO, ImageCore, ImageTransformations

pretrained_list = [
                    YOLO.v2_tiny_416_COCO,
                    # YOLO.v2_608_COCO,
                    YOLO.v3_tiny_416_COCO,
                    YOLO.v3_320_COCO,
                    YOLO.v3_416_COCO,
                    YOLO.v3_608_COCO,
                    # YOLO.v3_608_spp_COCO
                    ]

IMG = load(joinpath(@__DIR__,"images","dog-cycle-car.png"))

header = ["Model" "loaded?" "load time (s)" "ran?" "run time (s)" "objects detected"]
table = Array{Any}(undef, length(pretrained_list), 6)
for (i, pretrained) in pairs(pretrained_list)
    modelname = string(pretrained)
    table[i,:] = [modelname false "-" "-" "-" "-"]
    @testset "Pretrained Model: $modelname" begin
        global table

        t_load = @elapsed begin
            yolomod = pretrained(silent=true)
        end
        table[i, 2] = true
        table[i, 3] = round(t_load, digits=3)
        @info "$modelname: Loaded in $(round(t_load, digits=2)) seconds."

        batch = emptybatch(yolomod)
        batch[:,:,:,1] .= gpu(resizePadImage(IMG, yolomod))

        res = yolomod(batch) #run once
        val, t_run, bytes, gctime, memallocs = @timed yolomod(batch);
        table[i, 4] = true
        table[i, 5] = round(t_run, digits=4)
        table[i, 6] = size(res, 2)

        imgBoxes = drawBoxes(IMG, res)
        save(joinpath(@__DIR__,"results","$(modelname)_dog-cycle-car.jpg"), imgBoxes)

        @test size(res,2) > 0

        @info "$modelname: Ran in $(round(t_run, digits=2)) seconds. (bytes $bytes, gctime $gctime)"
    end
    GC.gc()
end
pretty_table(table, header)
@info "Times approximate. For more accurate benchmarking run ObjectDetector.benchmark()"


text:
__precompile__()
module SeqHax

using ArgParse
using Bio.Seq

include("ProgressLoggers.jl")
include("utils.jl")

include("comp.jl")
include("length.jl")
include("interleave.jl")

function parse_cli()
    s = ArgParseSettings()
    ## Global
    @add_arg_table s begin
        "comp"
            help = "Calculate nucleotide composition of sequences"
            action = :command
        "join"
            help = "Join separate R1/R2 files into an interleaved file"
            action = :command
        "length"
            help = "Count read lengths"
            action = :command
        "split"
            help = "Split an interleaved file into separate R1/R2 files"
            action = :command
    end
    Comp.add_args(s)
    Length.add_args(s)
    Interleave.add_join_args(s)
    Interleave.add_split_args(s)
    return parse_args(s)
end

function main()
    cli = parse_cli()
    cmd = cli["%COMMAND%"]
    mainfuncs = Dict{AbstractString, Any}(
        "comp" => Comp.main,
        "join" => Interleave.join_main,
        "length" => Length.main,
        "split" => Interleave.split_main,
        "preappend" => PreApp.main,
    )
    return mainfuncs[cmd](cli[cmd])
end

end # module SeqHax


text:
function trees_from_arg(arg::ARG)
	trees = Vector{Any}(undef, 2)
	for clr in 1:2
		trees[clr] = tree_from_arg(arg, clr)
	end
	return (trees[1], trees[2])
end

function tree_from_arg(arg::ARG, clr::Int)
	r = tree_from_arg!(nothing, arg.roots[clr], clr)
	return TreeTools.node2tree(r)
end

function tree_from_arg!(t_a, a_n::ARGNode, clr)
	@assert hascolor(a_n, clr)
	# Create new tree node
	n = TreeTools.TreeNode(; label = label(a_n), tau = branch_length(a_n, clr))
	if !isroot(a_n, clr)
		@assert !isnothing(t_a)
		TreeTools.graftnode!(t_a, n)
	end
	# Recurse
	for a_c in children(a_n, clr)
		tree_from_arg!(n, a_c, clr)
	end

	return n
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
using AdventOfCode2018
using AdventOfCode2018.Day6

inp = get_aoc_input(6)

coords = parse_input(inp)
println("Solution part 1: ", part_1(coords))
println("Solution part 2: ", part_2(coords))


text:
using Turing: @VarName, invlogit, logit, randcat, kl, align
using Distributions: Normal
using Test

i = 1
@test @VarName(s)[1:end-1] == (:s,())
@test @VarName(x[1,2][1+5][45][3][i])[1:end-1] == (:x,([1,2],[6],[45],[3],[1]))
@test invlogit(1.1) == 1.0 / (exp.(-1.1) + 1.0)
@test logit(0.3) ≈ -0.8472978603872036 atol=1e-9
@test isnan.(logit(1.0)) == false
@test isinf.(logit(1.0)) == true
@test isnan.(logit(0.)) == false
@test isinf.(logit(0.)) == true
randcat([0.1, 0.9])
@test kl(Normal(0, 1), Normal(0, 1)) == 0
@test align([1, 2, 3], [1]) == ([1,2,3],[1,0,0])
@test align([1], [1, 2, 3]) == ([1,0,0],[1,2,3])


text:
"""
    AbstractLinRange

Supertype for mutable or static ranges with `len` linearly spaced elements
between its `start` and `stop`.
"""
abstract type AbstractLinRange{T} <: AbstractRange{T} end

Base.firstindex(::AbstractLinRange) = 1

"""
    LinSRange

A static range with len linearly spaced elements between its start and stop.
The size of the spacing is controlled by len, which must be an Int.
"""
struct LinSRange{T,B,E,L,D} <: AbstractLinRange{T}
    function LinSRange{T}(start, stop, len) where T
        len >= 0 || throw(ArgumentError("srange($start, stop=$stop, length=$len): negative length"))
        if len == 1
            start == stop || throw(ArgumentError("srange($start, stop=$stop, length=$len): endpoints differ"))
            return new{T,start, stop, 1, 1}()
        end
        return new{T, T(start), T(stop), len, max(len-1,1)}()
    end
end

function LinSRange(start, stop, len::Integer)
    return LinSRange{typeof((stop-start)/len)}(start, stop, len)
end

function Base.getproperty(r::LinSRange, s::Symbol)
    if s === :start
        return first(r)
    elseif s === :stop
        return last(r)
    elseif s === :len
        return length(r)
    elseif s === :lendiv
        return lendiv(r)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

LinSRange{T}(r::AbstractRange) where {T} = LinSRange{T}(first(r), last(r), length(r))

"""
    LinMRange

A mutable range with len linearly spaced elements between its start and stop.
The size of the spacing is controlled by len, which must be an Int.
"""
mutable struct LinMRange{T} <: AbstractLinRange{T}
    start::T
    stop::T
    len::Int
    lendiv::Int

    function LinMRange{T}(start, stop, len) where T
        len >= 0 || throw(ArgumentError("mrange($start, stop=$stop, length=$len): negative length"))
        if len == 1
            start == stop || throw(ArgumentError("mrange($start, stop=$stop, length=$len): endpoints differ"))
            return new(start, stop, 1, 1)
        end
        return new(start, stop, len, max(len-1,1))
    end
end

function LinMRange(start, stop, len::Integer)
    return LinMRange{typeof((stop-start)/len)}(start, stop, len)
end

LinMRange{T}(r::AbstractRange) where {T}= LinMRange{T}(first(r), last(r), length(r))

function Base.setproperty!(r::LinMRange, s::Symbol, val)
    if s === :start
        return set_first!(r, val)
    elseif s === :stop
        return set_last!(r, val)
    elseif s === :len
        return set_length!(r, val)
    elseif s === :lendiv
        return set_lendiv!(r, val)
    else
        error("type $(typeof(r)) has no property $s")
    end
end

for (F,f) in ((:M,:m), (:S,:s))
    LR = Symbol(:Lin, F, :Range)
    frange = Symbol(f, :range)

    @eval begin
        Base.:(-)(r::$(LR)) = $(LR)(-first(r), -last(r), length(r))

        $(LR){T}(r::$(LR){T}) where {T} = r
        #$(LR){T}(r::AbstractRange) where {T} = $(LR){T}(first(r), last(r), length(r))
        $(LR)(r::AbstractRange{T}) where {T} = $(LR){T}(r)

        function Base.:(-)(r1::$(LR){T}, r2::$(LR){T}) where T
            len = length(r1)
            (len == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match")))
            return $(LR){T}(
                convert(T, -(first(r1), first(r2))),
                convert(T, -(last(r1), last(r2))),
                len
               )
        end

        function Base.:(+)(r1::$(LR){T}, r2::$(LR){T}) where T
            len = length(r1)
            (len == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match")))
            return $(LR){T}(
                convert(T, +(first(r1), first(r2))),
                convert(T, +(last(r1), last(r2))),
                len
               )
        end
    end
end


text:
nin(v::AbstractVertex) = nin(trait(v), v)
nin(t::DecoratingTrait, v::AbstractVertex) = nin(base(t), v)
nin(t::FixedSizeTrait, v::AbstractVertex) = nin(t, base(v), v)

nin(t::MutationTrait, v::AbstractVertex, vo::AbstractVertex=v) = nin(t, base(v), vo)
nin(t::MutationTrait, v::CompVertex, vo::AbstractVertex) = nin(v.computation, t, vo)
nin(::MutationTrait, ::InputSizeVertex) = []

# Give users a chance to fallback to smoe generic property of the vertex
nin(f, ::MutationTrait, ::AbstractVertex) = nin(f)
nin(f, ::FixedSizeTrait, ::AbstractVertex) = nin(f)
nin(f, ::SizeTransparent, v::AbstractVertex) = nout.(inputs(v))

nout(v::AbstractVertex) = nout(trait(v), v)
nout(t::DecoratingTrait, v::AbstractVertex) = nout(base(t), v)
nout(t::FixedSizeTrait, v::AbstractVertex) = nout(t, base(v), v)

nout(t::MutationTrait, v::AbstractVertex, vo::AbstractVertex=v) = nout(t, base(v), vo)
nout(t::MutationTrait, v::CompVertex, vo::AbstractVertex) = nout(v.computation, t, vo)
nout(::MutationTrait, v::InputSizeVertex) = v.size

# Give users a chance to fallback to smoe generic property of the vertex or trait
nout(f, ::MutationTrait, ::AbstractVertex) = nout(f)

# SizeTransparent might not care about size
nout(f, ::SizeInvariant, v::AbstractVertex) = isempty(nin(v)) ? 0 : nin(v)[1]
nout(f, ::SizeStack, v::AbstractVertex) = isempty(nin(v)) ? 0 : sum(nin(v))

# Add possibly for shape traits for computations
nout(f) = nout(shapetrait(f), f)
nin(f) = nin(shapetrait(f), f)

shapetrait(f::T) where T = T

struct NoShapeTraitError <: Exception
    msg::String
end
Base.showerror(io::IO, e::NoShapeTraitError) = print(io, e.msg)

noutnin_errmsg(f, name) = "$name not defined for $(f)! Either implement NaiveNASlib.$name(::$(typeof(f))) or implement NaiveNASlib.$name(::ST, ::$(typeof(f))) where ST is the return type of NaiveNASlib.shapetrait(::$(typeof(f)))"

nout(::Any, f) = throw(NoShapeTraitError(noutnin_errmsg(f, nout)))
nin(::Any, f) = throw(NoShapeTraitError(noutnin_errmsg(f, nin)))

"""
    ScalarSize

Treat vertices as having a scalar size when formulating the size change problem.
"""
struct ScalarSize end

# Main entry point for computing and applying size changes with ScalarSize

function Δsize!(case::ScalarSize, s::AbstractΔSizeStrategy, vertices::AbstractVector{<:AbstractVertex})
    execute, nins, nouts = newsizes(s, vertices)
    if execute
        applyΔsize!(case, s, nins, nouts, vertices)
    end
    return execute
end

# Main entry point for applying size changes with ScalarSize after we have compted them
# Peel off all DecoratingJuMPΔSizeStrategies until we find either an AbstractAfterΔSizeStrategy (which is needed in after_Δnin etc)
# or we hit the bottom (in the form of an AbstractΔSizeStrategy).

function applyΔsize!(case::ScalarSize, s::DecoratingJuMPΔSizeStrategy, nins::AbstractDict, nouts::AbstractVector, vertices::AbstractVector{<:AbstractVertex})
    applyΔsize!(case, base(s), nins, nouts, vertices)
end
function applyΔsize!(case::ScalarSize, s::AbstractAfterΔSizeStrategy, nins::AbstractDict, nouts::AbstractVector, vertices::AbstractVector{<:AbstractVertex})
    _applyΔsize!(case, s, nins, nouts, vertices)
end
function applyΔsize!(case::ScalarSize, s::AbstractΔSizeStrategy, nins::AbstractDict, nouts::AbstractVector, vertices::AbstractVector{<:AbstractVertex})
    _applyΔsize!(case, s, nins, nouts, vertices)
end

function _applyΔsize!(case::ScalarSize, s, nins::AbstractDict, nouts::AbstractVector, vertices::AbstractVector{<:AbstractVertex})
    Δnouts = nouts .- nout.(vertices)
    Δnins = [coalesce.(get(() -> nin(vi), nins, vi), nin(vi)) .- nin(vi) for vi in vertices]

    for (i, vi) in enumerate(vertices)
        insizes = get(() -> nin(vi), nins, vi)
        insizes = coalesce.(insizes, nin(vi))
        applyΔsize!(case, s, vi, insizes, nouts[i])
    end

    for (i, vi) in enumerate(vertices)
        after_Δnin(s, vi, Δnins[i], any(!=(0), Δnins[i]))
        after_Δnout(s, vi, Δnouts[i], Δnouts[i] != 0)
    end
end

function applyΔsize!(::ScalarSize, v::InputSizeVertex, insizes::AbstractVector{<:Integer}, outs::Integer; kwargs...)
    if !all(isempty, skipmissing(ins))
        throw(ArgumentError("Try to change input size of InputVertex $(name(v)) to $insizes"))
    end
    if outsize != nout(v)
        throw(ArgumentError("Try to change output size of InputVertex $(name(v)) to $outsize"))
    end
end

# Main entry point for applying size changes to a single vertex with ScalarSize

function applyΔsize!(case::ScalarSize, s::WithKwargs, v::AbstractVertex, insizes::AbstractVector{<:Integer}, outsize::Integer; kws...)
    applyΔsize!(case, v, insizes, outsize; kws..., s.kwargs...)
end
function applyΔsize!(case::ScalarSize, s::DecoratingJuMPΔSizeStrategy, v::AbstractVertex, insizes::AbstractVector{<:Integer}, outsize::Integer; kws...)
    applyΔsize!(case, base(s), v, insizes, outsize; kws...)
end
function applyΔsize!(case::ScalarSize, s::AbstractΔSizeStrategy, v::AbstractVertex, insizes::AbstractVector{<:Integer}, outsize::Integer; kws...)
    applyΔsize!(case, v, insizes, outsize; kws...)
end
function applyΔsize!(case::ScalarSize, v::AbstractVertex, insizes::AbstractVector{<:Integer}, outsize::Integer; kws...)
    applyΔsize!(case, base(v), insizes, outsize; kws...)
end
function applyΔsize!(case::ScalarSize, v::CompVertex, insizes::AbstractVector{<:Integer}, outsize::Integer; kws...)
    # I don't love how we go back to Δsize! from here, but I don't want to export applyΔsize! and I can't think of another name right now...
    Δsize!(case, v.computation, insizes, outsize; kws...)
end
Δsize!(::ScalarSize, f, insizes::AbstractVector{<:Integer}, outsize::Integer; kws...) = Δsize!(f, insizes, outsize; kws...)
function Δsize!(f, ins::AbstractVector{<:Integer}, outs::Integer; kws...) end

after_Δnin(::AbstractΔSizeStrategy, v, Δs, changed) = after_Δnin(trait(v), v, Δs, changed)
after_Δnin(t::DecoratingTrait, v, Δs, changed) = after_Δnin(base(t), v, Δs, changed)
after_Δnin(s::AbstractAfterΔSizeStrategy, args...) = _after_Δnin(s, args...)
function after_Δnin(t::AfterΔSizeTrait, v, Δs, changed)
    after_Δnin(t.strategy, v, Δs, changed)
    after_Δnin(base(t), v, Δs, changed)
end
_after_Δnin(s::DecoratingJuMPΔSizeStrategy, args...) = _after_Δnin(base(s), args...)
function _after_Δnin(::AbstractΔSizeStrategy, args...) end
function _after_Δnin(s::AfterΔSizeCallback, v, Δs, changed)
    s.cbfun(v, Δs, :nin, changed)
    _after_Δnin(base(s), v, Δs, changed)
end
function after_Δnin(t, v, Δs, changed) end

after_Δnout(::AbstractΔSizeStrategy, v, Δ, changed) = after_Δnout(trait(v), v, Δ, changed)
after_Δnout(t::DecoratingTrait, v, Δ, changed) = after_Δnout(base(t), v, Δ, changed)
after_Δnout(s::AbstractAfterΔSizeStrategy, args...) = _after_Δnout(s, args...)
function after_Δnout(t::AfterΔSizeTrait, v, Δs, changed)
    after_Δnout(t.strategy, v, Δs, changed)
    after_Δnout(base(t), v, Δs, changed)
end

_after_Δnout(s::DecoratingJuMPΔSizeStrategy, args...) = _after_Δnout(base(s), args...)
function _after_Δnout(::AbstractΔSizeStrategy, args...) end
function _after_Δnout(s::AfterΔSizeCallback, v, Δ, changed)
    s.cbfun(v, Δ, :nout, changed)
    _after_Δnout(base(s), v, Δ, changed)
end
function after_Δnout(t, v, Δ, changed) end

newsizes(s::ThrowΔSizeFailError, vs::AbstractVector{<:AbstractVertex}) = throw(ΔSizeFailError(s.msgfun(vs)))
newsizes(::ΔSizeFailNoOp, vs::AbstractVector{<:AbstractVertex}) = false, Dict(v => nin(v) for v in vs), nout.(vs)
function newsizes(s::LogΔSizeExec, vertices::AbstractVector{<:AbstractVertex})
    @logmsg s.level s.msgfun(vertices[1])
    return newsizes(base(s), vertices)
end

"""
    newsizes(s::AbstractΔSizeStrategy, vertices::AbstractArray{<:AbstractVertex})

Return a vector of new outputs sizes for and a `Dict` of new input sizes for all provided `vertices` using the strategy `s`.

Result vector is index aligned with `vertices`.
Result `Dict` has a vector of input sizes for each element of `vertices` which has an input (i.e everything except input vertices).
"""
function newsizes(s::AbstractJuMPΔSizeStrategy, vertices::AbstractVector{<:AbstractVertex})
    case = ScalarSize()
    model = sizemodel(s, vertices)

    noutvars = @variable(model, noutvars[i=1:length(vertices)], Int)

    noutdict = Dict(zip(vertices, noutvars))
    for v in vertices
        vertexconstraints!(case, v, s, (;model, noutdict))
    end

    sizeobjective!(case, s, vertices, (;model, noutdict))

    JuMP.optimize!(model)

    if accept(case, s, model)
        return true, ninsandnouts(s, vertices, noutvars)...
    end
    fbstrat = fallback(s)
    return newsizes(fbstrat, add_participants!(fbstrat, copy(vertices)))
end

"""
    sizemodel(s::AbstractJuMPΔSizeStrategy, vertices)

Return a `JuMP.Model` for executing strategy `s` on `vertices`.
"""
sizemodel(::AbstractJuMPΔSizeStrategy, vertices) = JuMP.Model(JuMP.optimizer_with_attributes(Cbc.Optimizer, "loglevel"=>0))
function sizemodel(s::TimeLimitΔSizeStrategy, vs)
    model = sizemodel(base(s), vs)
    JuMP.set_time_limit_sec(model, s.limit)
    return model
end

"""
    accept(case, ::AbstractJuMPΔSizeStrategy, model::JuMP.Model)

Return true of the solution for `model` is accepted using strategy `s`.
"""
accept(case, s::DecoratingJuMPΔSizeStrategy, model::JuMP.Model) = accept(case, base(s), model)
accept(case, ::AbstractJuMPΔSizeStrategy, model::JuMP.Model) = JuMP.termination_status(model) != MOI.INFEASIBLE && JuMP.primal_status(model) == MOI.FEASIBLE_POINT # Beware: primal_status seems unreliable for Cbc. See MathOptInterface issue #822
function accept(case, s::TimeOutAction, model::JuMP.Model)
    if  JuMP.termination_status(model) == MOI.TIME_LIMIT
        return s.action(model)
    end
    return accept(case, base(s), model)
end

# Just a shortcut for broadcasting on dicts
getall(d::Dict, ks, deffun=() -> missing) = get.(deffun, [d], ks)

# First we dispatch on trait in order to filter out immutable vertices
# First dispatch on traits to sort out things like immutable vertices
vertexconstraints!(case, v::AbstractVertex, s::AbstractJuMPΔSizeStrategy, data) = vertexconstraints!(case, trait(v), v, s, data)
vertexconstraints!(case, t::DecoratingTrait, v, s::AbstractJuMPΔSizeStrategy, data) = vertexconstraints!(case, base(t), v, s, data)
vertexconstraints!(case, ::MutationSizeTrait, v, s::AbstractJuMPΔSizeStrategy, data) = vertexconstraints!(case, s, v, data) # Now dispatch on strategy

function vertexconstraints!(case::ScalarSize, ::Immutable, v, s, data)
    @constraint(data.model, data.noutdict[v] == nout(v))
    @constraint(data.model, getall(data.noutdict, inputs(v)) .== nin(v))
    vertexconstraints!(case, s, v, data)
end

# This must be applied to immutable vertices as well
function vertexconstraints!(currcase::ScalarSize, v::AbstractVertex, s::AlignNinToNout, data, case=currcase)
    vertexconstraints!(case, v, s.vstrat, data)
    # Code below secretly assumes vo is in data.noutdict (ninarr will be left with undef entries otherwise).
    # ninvar shall be used to set new nin sizes after problem is solved
    for vo in filter(vo -> vo in keys(data.noutdict), outputs(v))
        ninvar = @variable(data.model, integer=true)
        @constraint(data.model, data.noutdict[v] == ninvar)

        ninarr = get!(() -> Vector{JuMP.VariableRef}(undef, length(inputs(vo))), s.nindict, vo)
        ninarr[inputs(vo) .== v] .= ninvar
    end
end

function vertexconstraints!(currcase::ScalarSize, v::AbstractVertex, s::AlignNinToNoutVertices, data, case=currcase)
    # Any s.ininds which are mapped to vertices which are not yet added in s.vstrat.nindict[s.vout] will cause undefined reference below
    # There is a check to wait for them to be added, but if they are not in data.noutdict they will never be added, so we need to check for that too. Not 100% this can even happen, so I'm just waiting for this in itself to trigger a bug. Sigh...
    neededinds = filter(i -> i !== nothing, indexin(keys(data.noutdict), inputs(s.vout)))

    condition() = s.vout in keys(s.vstrat.nindict) && all(i -> isassigned(s.vstrat.nindict[s.vout], i), neededinds)

    # Just to make sure we only do this once without having to store hasadded as a field in AlignNinToNoutVertices:
    #  - Check the condition before vertexconstraints! -> Only when condition changes from false to true do we add as we assume this is the first time the condition becomes fulfilled
    hasadded = condition()
    vertexconstraints!(case, v, s.vstrat, data)

    if !hasadded && condition()
        @constraint(data.model, data.noutdict[s.vin] .== s.vstrat.nindict[s.vout][s.ininds])
    end
end

"""
    vertexconstraints!(s::AbstractJuMPΔSizeStrategy, v, data)

Add constraints for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
function vertexconstraints!(case::ScalarSize, s::AbstractJuMPΔSizeStrategy, v, data)
    ninconstraint!(case, s, v, data)
    compconstraint!(case, s, v, (data..., vertex=v))
    sizeconstraint!(case, s, v, data)
end

"""
    sizeconstraint!(s::AbstractJuMPΔSizeStrategy, v, data)

Add size constraints for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
sizeconstraint!(case::ScalarSize, s::DecoratingJuMPΔSizeStrategy, v, data) = sizeconstraint!(case, base(s), v, data)
function sizeconstraint!(::ScalarSize, s::AbstractJuMPΔSizeStrategy, v, data)
    # If we find any outputs to v which are not part of the problem to solve, we need to prevent v from changing size
    if any(vo -> vo ∉ keys(data.noutdict), outputs(v))
        # Ok, this is a bit speculative, but should we find that there is a size mismatch (e.g. when create/remove vertex/edge)
        # we might make the problem infeasible somehow by doing this. Instead we just skip it and hope some other constraint deals
        # with aligning the size
        allnins = unique!(mapreduce(vo -> nin(vo)[v .== inputs(vo)], vcat, outputs(v)))
        if length(allnins) == 1
            @constraint(data.model, data.noutdict[v] == allnins[1])
            return
        end
    end
    @constraint(data.model, data.noutdict[v] >= 1)
end

function sizeconstraint!(case::ScalarSize, s::ΔNout{Exact}, v, data)
    Δ = get(s.Δs, v, nothing)
    if Δ !== nothing
        @constraint(data.model, data.noutdict[v] == nout(v) + Δ)
    else
        sizeconstraint!(case, DefaultJuMPΔSizeStrategy(), v, data)
    end
end

"""
    ninconstraint!(case, s, v, data)

Add input size constraints for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
ninconstraint!(case, s::DecoratingJuMPΔSizeStrategy, v, data) = ninconstraint!(case, base(s), v, data)
ninconstraint!(case, s, v, data) = ninconstraint!(case, s, trait(v), v, data)
ninconstraint!(case, s, t::DecoratingTrait, v, data) = ninconstraint!(case, s, base(t), v, data)
function ninconstraint!(case, s, ::MutationTrait, v, data) end
ninconstraint!(case, s, ::SizeStack, v, data) = @constraint(data.model, sum(getall(data.noutdict, inputs(v))) == data.noutdict[v])
ninconstraint!(case, s, ::SizeInvariant, v, data) = @constraint(data.model, getall(data.noutdict, unique(inputs(v))) .== data.noutdict[v])

"""
    compconstraint!(case, s, v, data)

Add constraints on the computation (e.g. neural network layer) for [`AbstractVertex`](@ref) `v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
compconstraint!(case, s, v::AbstractVertex, data) = compconstraint!(case, s, base(v), data)
compconstraint!(case, s, v::CompVertex, data) = compconstraint!(case, s, v.computation, data)
function compconstraint!(case, s, v::InputVertex, data) end
function compconstraint!(case, s, f, data) end

"""
    sizeobjective!(case, s::AbstractJuMPΔSizeStrategy, model, noutvars, sizetargets)

Add the objective for `noutvars` using strategy `s`.
"""
sizeobjective!(case::ScalarSize, s::DecoratingJuMPΔSizeStrategy, vertices, data) = sizeobjective!(case, base(s), vertices, data)
sizeobjective!(case::ScalarSize, s::AbstractJuMPΔSizeStrategy, vertices, data) = @objective(data.model, Min, objective!(case, s, vertices, data))

function objective!(::ScalarSize, s, vertices, data)
    model = data.model
    noutvars = map(v -> data.noutdict[v], vertices)
    sizetargets = nout.(vertices)
    # L1 norm prevents change in vertices which does not need to change.
    # Max norm tries to spread out the change so no single vertex takes most of the change.
    return norm!(SumNorm(0.1 => L1NormLinear(), 0.8 => MaxNormLinear()), model, @expression(model, objective[i=1:length(noutvars)], noutvars[i] - sizetargets[i]), sizetargets)
end

objective!(case::ScalarSize, s::ΔNout{Relaxed}, vertices, data) = noutrelax!(case, s.Δs, vertices, data)

function noutrelax!(case, Δs, vertices, data)
    def_obj = objective!(case, DefaultJuMPΔSizeStrategy(), setdiff(vertices, keys(Δs)), data)

    model = data.model
    # Force it to change as Δ might be too small (unless Δ is exactly 0)
    # Trick from http://lpsolve.sourceforge.net/5.1/absolute.htm
    Δnout_const = @expression(model, [data.noutdict[v] - nout(v) for v in keys(Δs) if Δs[v] != 0])

    B = @variable(model, [1:length(Δnout_const)], Bin)
    M = 1e5
    ϵ = 1e-2 # abs(Δnout_const) must be larger than this
    @constraint(model, Δnout_const .+ M .* B .>= ϵ)
    @constraint(model, Δnout_const .+ M .* B .<= M .- ϵ)

    sizediff = @expression(model, [data.noutdict[v] - nout(v) - Δ for (v, Δ) in Δs])
    Δnout_obj = norm!(SumNorm(0.1 => L1NormLinear(), 0.8 => MaxNormLinear()), model, sizediff)

    return @expression(model, def_obj + 1e6*sum(Δnout_obj))
end

function ninsandnouts(::AbstractJuMPΔSizeStrategy, vs, noutvars)
    nouts = round.(Int, JuMP.value.(noutvars))
    mapnout(i::Integer) = nouts[i]
    mapnout(i::Nothing) = missing

    nins = Dict(vs .=> map(vi -> mapnout.(indexin(inputs(vi), vs)), vs))
    return nins, nouts
end

function ninsandnouts(s::AlignNinToNout, vs, noutvars)
    nouts = round.(Int, JuMP.value.(noutvars))
    nins = Dict(key => round.(Int, JuMP.value.(value)) for (key, value) in s.nindict)
    return nins,nouts
end

ninsandnouts(s::AlignNinToNoutVertices, vs, noutvars) = ninsandnouts(s.vstrat, vs, noutvars)

# TODO: Remove since only used for debugging. If only it wasn't so bloody cumbersome to just list the constraints in a JuMP model....
nconstraints(model) = mapreduce(tt -> JuMP.num_constraints.(model,tt...), +,  filter(tt -> tt != (JuMP.VariableRef, MOI.Integer), JuMP.list_of_constraint_types(model)), init=0)

# TODO: Remove since only used for debugging. If only it wasn't so bloody cumbersome to just list the constraints in a JuMP model....
function list_constraints(model)
    for tt in filter(tt -> tt != (JuMP.VariableRef, MOI.Integer), JuMP.list_of_constraint_types(model))
        display(JuMP.all_constraints(model, tt...))
    end
end


text:
# some functions used by all physics modules, particularly durin ginitialization

module SolverCommon

using PDESolver
using Utils
using ODLCommonTools
using LinearSolvers
using NonlinearSolvers
using Input
using MPI
using SummationByParts
using PdePumiInterface

# from common.jl
export createMeshAndOperator, loadRestartState, call_nlsolver, getDataTypes

include("common.jl")

end # end module


text:
struct SpriteImage
    pattern::Cairo.CairoPattern
    width::Real
    height::Real
end

function SpriteImage(pattern::Cairo.CairoPattern; filter = Cairo.FILTER_NEAREST)
    surface = Cairo.pattern_get_surface(pattern)
    w, h = surface.width, surface.height
    Cairo.pattern_set_filter(pattern, filter)
    return SpriteImage(pattern, w, h)
end

function loadsprite(path::String, filter = Cairo.FILTER_NEAREST)
    img = readpng(path)
    pattern = Cairo.CairoPattern(img)
    return SpriteImage(pattern, filter = filter)
end

function drawsprite(sprite::SpriteImage,
    xpos::Real = 0, ypos::Real = 0,
    width::Real = sprite.width,
    height::Real = sprite.height; centered = false)

    gsave()
    cr = Luxor.get_current_cr()
    w, h = sprite.width, sprite.height

    if centered
        xpos, ypos = xpos - (width/2), ypos - (height/2)
    end

    translate(xpos, ypos)
    scale(width/w, height/h)
    Cairo.set_source(cr, sprite.pattern)
    Cairo.paint(cr)
    grestore()
end

function drawsprite(sprite::SpriteImage, pos::Point, size::UnitRange; centered = false)
    drawsprite(sprite, pos.x, pos.y, size.start, size.stop, centered = centered)
end


text:
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

function get_axis_rulers()
    label_style = cv_black → cv_fontface("sans-serif") → cv_fontsize(20)

    ticks1, ticks1h = "%.0f" ⇒ -4:4, "" ⇒ -4.5:4.5

    app_l = CV_TickLabelAppearance(; label_style, tick_length=10)
    app_s = CV_TickLabelAppearance(; label_style, tick_length=6)

    return (app_l ↦ ticks1, app_s ↦ ticks1h,)
end

function create_colorbar(setup, can_codomain_l, winding_painter) # {{{
    layout = setup.layout
    width = cv_width(can_codomain_l.rectangle)
    mid = cv_half(can_codomain_l.rectangle.right + can_codomain_l.rectangle.left)
    bottom = can_codomain_l.rectangle.bottom
    height = 30

    colorbar_data = cv_create_winding_colorbar(width, height,
        winding_painter, -2, 5)
    cont_colorbar = colorbar_data.colorbar_container
    cont_colorbar_l = cv_add_canvas!(layout,
        cont_colorbar,
        cv_anchor(cont_colorbar, :slider_south, cont_colorbar, :south),
        (mid, bottom - Int32(20)))

    return cv_setup_winding_colorbar(setup, colorbar_data, cont_colorbar_l)
end # }}}

function create_scene(trafo, can_codomain, curve_painter; padding=30) # {{{
    layout = CV_StateLayout(CV_2DLayout(), CV_CyclicValue(2))

    can_codomain_l = cv_add_canvas!(layout, can_codomain, (0, 0), (0, 0))

    cc_can_codomain = cv_create_context(can_codomain)

    rulers = get_axis_rulers()
    codomain_re_axis = cv_ticks_labels(layout, can_codomain_l, cv_south, rulers)
    codomain_im_axis = cv_ticks_labels(layout, can_codomain_l, cv_west, rulers)
    codomain_border = cv_border(layout, can_codomain_l, 2)

    fill_painter = CV_FillPainter()
    styled_fill_painter = cv_color(1,1,1) ↦ fill_painter

    wind_painter = CV_WindingPainter(trafo, curve_painter.segments)

    style = cv_color(0,0,1) → cv_linewidth(4) → cv_antialias_best
    styled_curve_painter = style ↦ curve_painter

    dir_painter = CV_DirectionPainter(trafo,
        curve_painter.segments, true; every_len=3)
    dir_style = cv_color(0,0,1) → cv_op_source
    styled_dir_painter = dir_style ↦ dir_painter

    setup = CV_MinimalSetupChain(layout)

    actionpixel_update = (px, py, future_layout) -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        cv_paint(cc_can_codomain, styled_fill_painter)
        cv_paint(cc_can_codomain, wind_painter)
        cv_paint(cc_can_codomain, styled_curve_painter)
        cv_paint(cc_can_codomain, styled_dir_painter)
        can_codomain_l(cc_can_layout)
        return nothing
    end

    draw_once_func = future_layout -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        codomain_re_axis(cc_can_layout)
        codomain_im_axis(cc_can_layout)
        codomain_border(cc_can_layout)
        return nothing
    end
    setup = cv_combine(setup; draw_once_func, actionpixel_update)

    setup = create_colorbar(setup, can_codomain_l, wind_painter)
    padding > 0 && cv_add_padding!(setup.layout, padding)
    setup = cv_setup_2dminimal_scene(setup)
    setup.layout.actionpixel_update(Int32(0), Int32(0))
    setup.layout.redraw_func()
    return setup.layout
end # }}}

function create_curve()
    phi = collect(LinRange(0.0, 2*pi, 700))[1:end-1]
    curve = 2.0*exp.(1im.*phi)
    return [curve,] :: CV_LineSegments
end

trafo = z -> sin(z^2)/10 - 5/z
codomain = CV_Math2DCanvas(-5.0 + 5.0im, 5.0 - 5.0im, 80)
curve_painter = CV_LinePainter(trafo, create_curve(), true)

scene = create_scene(trafo, codomain, curve_painter)

cv_save_image(cv_get_can_layout(scene), "./Winding01.png")

if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
function ed_model(ed_system, optimizer, step::Int64=1)

    m = JuMP.Model(optimizer)

    #Time Information
    time_periods = PowerSystems.get_forecasts_horizon(ed_system)
    time_periods_set = 1:time_periods
    data_first_step = PowerSystems.get_forecast_initial_times(ed_system)[step]
    minutes_per_period = Dates.Minute(PowerSystems.get_forecasts_resolution(ed_system))/Dates.Minute(1)

    # Thermal Generation
    thermal_generators = PowerSystems.get_components(PowerSystems.ThermalStandard, ed_system)
    gen_pwl_points = Dict(PowerSystems.get_name(g) => 1:length(g.op_cost.variable) for g in thermal_generators)
    thermal_gen_names = [PowerSystems.get_name(g) for g in thermal_generators]

    JuMP.@variable(m, cg[thermal_gen_names,time_periods_set])
    JuMP.@variable(m, pg[thermal_gen_names,time_periods_set] >= 0)
    JuMP.@variable(m, rg[thermal_gen_names,time_periods_set] >= 0)
    JuMP.@variable(m, ug[thermal_gen_names,time_periods_set] >= 0)
    JuMP.@variable(m, slack_reserve[time_periods_set] >= 0)
    JuMP.@variable(m, slack_load[time_periods_set] >= 0)
    JuMP.@variable(m, slack_generation[time_periods_set] >= 0)
    JuMP.@variable(m, 0 <= lambda_lg[g in thermal_gen_names, gen_pwl_points[g], time_periods_set] <= 1);

    # Renewable Generation
    renewable_forecasts = PowerSystems.get_component_forecasts(PowerSystems.RenewableDispatch, ed_system, data_first_step)
    renewable_gen_names = PowerSystems.get_forecast_component_name.(renewable_forecasts)
    JuMP.@variable(m, pw[renewable_gen_names,time_periods_set] >= 0)

    # Loads
    fix_load_forecasts =  PowerSystems.get_component_forecasts(PowerSystems.PowerLoad, ed_system, data_first_step)

    interruptible_load_forecasts =  PowerSystems.get_component_forecasts(PowerSystems.InterruptibleLoad, ed_system, data_first_step)
    interruptible_load_names = PowerSystems.get_forecast_component_name.(interruptible_load_forecasts)
    JuMP.@variable(m, pl[interruptible_load_names, time_periods_set] >= 0)

    JuMP.@objective(m, Min,
    sum(
        sum(cg[PowerSystems.get_name(g),t] + (PowerSystems.get_op_cost(g) |> PowerSystems.get_fixed )*ug[PowerSystems.get_name(g),t] for g in thermal_generators) + 1e6*(slack_reserve[t] + slack_load[t]+slack_generation[t])
       - sum(PowerSystems.get_component(il).op_cost.variable.cost[2]*pl[PowerSystems.get_forecast_component_name(il), t] for il in interruptible_load_forecasts)
       - sum(pw[PowerSystems.get_forecast_component_name(ren), t] for ren in renewable_forecasts)
        for t in time_periods_set)
        )

    # Constraints for first time period that require initial conditions
    for g in thermal_generators
        name = PowerSystems.get_name(g)
        power_output_t0 = PowerSystems.get_activepower(g)
        unit_on_t0 = 1.0*(power_output_t0 > 0)
        activepowerlimits = PowerSystems.get_tech(g) |> PowerSystems.get_activepowerlimits
        time_minimum = PowerSystems.get_tech(g) |> PowerSystems.get_timelimits
        ramplimits = PowerSystems.get_tech(g) |> PowerSystems.get_ramplimits

        # Ramp Constraints
        JuMP.@constraint(m, pg[name,1] + rg[name,1] - unit_on_t0*(power_output_t0 - activepowerlimits.min) <= ramplimits.up*minutes_per_period)

        JuMP.@constraint(m, unit_on_t0*(power_output_t0 - activepowerlimits.min) - pg[name,1] <= ramplimits.down*minutes_per_period)

    end

    #Remaning Constraints.

    for t in time_periods_set

        JuMP.@constraint(m,
            sum( pg[PowerSystems.get_name(g),t] + g.tech.activepowerlimits.min*ug[PowerSystems.get_name(g),t] for g in thermal_generators) + slack_generation[t] +
            sum( pw[PowerSystems.get_forecast_component_name(g),t] for g in renewable_forecasts)
            == sum(PowerSystems.get_component(load).maxactivepower*PowerSystems.get_forecast_value(load, t) for load in fix_load_forecasts) + sum(pl[PowerSystems.get_forecast_component_name(l),t] for l in interruptible_load_forecasts) + slack_load[t]
        )

        for il in interruptible_load_forecasts
            load_value = PowerSystems.get_component(il).maxactivepower*PowerSystems.get_forecast_value(il, t)
            JuMP.set_upper_bound(pl[PowerSystems.get_forecast_component_name(il), t], load_value)
        end

        for reserve in PowerSystems.get_component_forecasts(PowerSystems.StaticReserve, ed_system, data_first_step)
            JuMP.@constraint(m, sum(rg[name,t] for name in thermal_gen_names) + slack_reserve[t] >= PowerSystems.get_component(reserve).requirement*PowerSystems.get_forecast_value(reserve, t)) # (3)
        end

         for g in thermal_generators
            name = PowerSystems.get_name(g)
            power_output_t0 = PowerSystems.get_activepower(g)
            unit_on_t0 = 1.0*(power_output_t0 > 0)
            activepowerlimits = PowerSystems.get_tech(g) |> PowerSystems.get_activepowerlimits
            time_minimum = PowerSystems.get_tech(g) |> PowerSystems.get_timelimits
            ramplimits = PowerSystems.get_tech(g) |> PowerSystems.get_ramplimits
            piecewise_production = PowerSystems.get_op_cost(g) |> PowerSystems.get_variable

            #JuMP.@constraint(m, pg[name,t] + rg[name,t] - pg[name,t-1] <= ramplimits.up*minutes_per_period) # (19)
            if t > 1
                JuMP.@constraint(m, pg[name,t] - pg[name,t-1] <= ramplimits.up*minutes_per_period) # (19)
                JuMP.@constraint(m, pg[name,t-1] - pg[name,t] <= ramplimits.down*minutes_per_period) # (20)
            end

            JuMP.@constraint(m, pg[name,t] + rg[name,t] <= (activepowerlimits.max - activepowerlimits.min)*ug[name,t])

            JuMP.@constraint(m, pg[name,t] == sum((piecewise_production[l][1] - piecewise_production[1][1])*lambda_lg[name,l,t] for l in gen_pwl_points[name])) # (21)
            JuMP.@constraint(m, cg[name,t] == sum((piecewise_production[l][2] - piecewise_production[1][2])*lambda_lg[name,l,t] for l in gen_pwl_points[name])) # (22)
            JuMP.@constraint(m, ug[name,t] == sum(lambda_lg[name,l,t] for l in gen_pwl_points[name])) # (23)
        end

        for rgen in renewable_forecasts
            name = PowerSystems.get_forecast_component_name(rgen)
            ub = rgen.component.tech.rating*PowerSystems.get_forecast_value(rgen,t)
            JuMP.set_upper_bound(pw[name,t], ub)
        end

    end

    return m

end


text:
using Pkg

# let julia handle python internally
#ENV["PYTHON"] = ""

# What has already been installed
packages = keys(Pkg.installed())

# add
package_names = [
    "FileIO",
    "JLD2",
    "GridInterpolations",
    "PyCall",
    "PyPlot",
    "HDF5",
    "LinearAlgebra"
]
for name in package_names
    if !in(name, packages)
        Pkg.add(name)
    end
end

# SISL packages

if !in("Vec", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/Vec.jl.git"))
end
if !in("Records", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/Records.jl.git"))
end
if !in("NGSIM", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/NGSIM.jl.git"))
end
if !in("BayesNets", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/BayesNets.jl.git"))
end
if !in("AutomotiveDrivingModels", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutomotiveDrivingModels.jl.git"))
end
if !in("AutoViz", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutoViz.jl.git"))
end
if !in("AutoRisk", packages)
    Pkg.add(PackageSpec(url="https://github.com/sisl/AutoRisk.jl.git", rev="v0.7fixes"))
end


text:
abstract type AbstractLShapedSolver{T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver} <: AbstractStructuredModel end

nscenarios(lshaped::AbstractLShapedSolver) = lshaped.nscenarios

# Initialization #
# ======================================================================== #
function init!(lshaped::AbstractLShapedSolver{T,A,M,S},subsolver::AbstractMathProgSolver) where {T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
    # Initialize progress meter
    lshaped.progress.thresh = lshaped.parameters.τ
    # Prepare the master optimization problem
    prepare_master!(lshaped)
    # Finish initialization based on solver traits
    init_solver!(lshaped)
    init_subproblems!(lshaped,subsolver)
end
# ======================================================================== #

# Functions #
# ======================================================================== #
function set_params!(lshaped::AbstractLShapedSolver; kwargs...)
    for (k,v) in kwargs
        setfield!(lshaped.parameters,k,v)
    end
end

function update_solution!(lshaped::AbstractLShapedSolver)
    ncols = lshaped.structuredmodel.numCols
    x = getsolution(lshaped.mastersolver)
    lshaped.mastervector[:] = x
    lshaped.x[1:ncols] = x[1:ncols]
    lshaped.θs[:] = x[end-lshaped.nscenarios+1:end]
    nothing
end

function update_structuredmodel!(lshaped::AbstractLShapedSolver)
    lshaped.structuredmodel.colVal = copy(lshaped.x)
    lshaped.structuredmodel.objVal = lshaped.c⋅lshaped.x + sum(lshaped.subobjectives)

    for i in 1:lshaped.nscenarios
        m = subproblem(lshaped.structuredmodel,i)
        m.colVal = copy(getsolution(lshaped.subproblems[i].solver))
        m.objVal = getobjval(lshaped.subproblems[i].solver)
    end
    nothing
end

function calculate_estimate(lshaped::AbstractLShapedSolver)
    return lshaped.c⋅lshaped.x + sum(lshaped.θs)
end

function current_objective_value(lshaped::AbstractLShapedSolver,Qs::AbstractVector)
    return lshaped.c⋅lshaped.x + sum(Qs)
end
current_objective_value(lshaped) = current_objective_value(lshaped,lshaped.subobjectives)

function get_solution(lshaped::AbstractLShapedSolver)
    return lshaped.x
end

function get_objective_value(lshaped::AbstractLShapedSolver)
    if !isempty(lshaped.Q_history)
        return lshaped.Q_history[end]
    else
        return calculate_objective_value(lshaped)
    end
end

function prepare_master!(lshaped::AbstractLShapedSolver)
    # θs
    for i = 1:lshaped.nscenarios
        addvar!(lshaped.mastersolver.lqmodel,-Inf,Inf,1.0)
        push!(lshaped.mastervector,-1e10)
    end
end

function resolve_subproblems!(lshaped::AbstractLShapedSolver{T,A,M,S}) where {T <: Real, A <: AbstractVector, M <: LQSolver, S <: LQSolver}
    # Update subproblems
    update_subproblems!(lshaped.subproblems,lshaped.x)

    # Solve sub problems
    for subproblem ∈ lshaped.subproblems
        cut::SparseHyperPlane{T} = subproblem()
        if !bounded(cut)
            warn("Subproblem ",cut.id," is unbounded, aborting procedure.")
            return -Inf
        end
        addcut!(lshaped,cut)
        lshaped.subobjectives[cut.id] = cut(lshaped.x)
    end

    # Return current objective value
    return current_objective_value(lshaped)
end

function iterate_nominal!(lshaped::AbstractLShapedSolver)
    # Resolve all subproblems at the current optimal solution
    Q = resolve_subproblems!(lshaped)
    if Q == -Inf
        return :Unbounded
    end
    lshaped.solverdata.Q = Q
    # Update the optimization vector
    take_step!(lshaped)
    # Resolve master
    try
        lshaped.mastersolver(lshaped.mastervector)
    catch
        # Master problem could not be solved for some reason. Try to continue
        return :Valid
    end
    if status(lshaped.mastersolver) == :Infeasible
        warn("Master is infeasible. Aborting procedure.")
        return :Infeasible
    end
    # Update master solution
    update_solution!(lshaped)
    lshaped.solverdata.θ = calculate_estimate(lshaped)
    # Log progress
    log!(lshaped)
    # Check optimality
    if check_optimality(lshaped)
        # Optimal
        lshaped.solverdata.Q = calculate_objective_value(lshaped,lshaped.x)
        push!(lshaped.Q_history,lshaped.solverdata.Q)
        return :Optimal
    end
    # Project (if applicable)
    project!(lshaped)
    # Just return a valid status for this iteration
    return :Valid
end

function log!(lshaped::AbstractLShapedSolver)
    @unpack Q,θ = lshaped.solverdata
    push!(lshaped.Q_history,Q)
    push!(lshaped.θ_history,θ)
    lshaped.solverdata.iterations += 1

    log_regularization!(lshaped)

    gap = abs(θ-Q)/(1+abs(Q))
    if lshaped.parameters.log
        ProgressMeter.update!(lshaped.progress,gap,
                              showvalues = [
                                  ("Objective",Q),
                                  ("Gap",gap),
                                  ("Number of cuts",length(lshaped.cuts))
                              ])
    end
end

function log!(lshaped::AbstractLShapedSolver,t::Integer)
    @unpack Q,θ = lshaped.solverdata
    lshaped.Q_history[t] = Q
    lshaped.θ_history[t] = θ

    log_regularization!(lshaped,t)
end

function check_optimality(lshaped::AbstractLShapedSolver)
    @unpack τ = lshaped.parameters
    @unpack Q,θ = lshaped.solverdata
    return θ > -Inf && abs(θ-Q) <= τ*(1+abs(Q))
end
# ======================================================================== #

# Cut functions #
# ======================================================================== #
active(lshaped::AbstractLShapedSolver,hyperplane::HyperPlane) = active(hyperplane,lshaped.x,lshaped.parameters.τ)
active(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut}) = optimal(cut,lshaped.x,lshaped.θs[cut.id],lshaped.parameters.τ)
satisfied(lshaped::AbstractLShapedSolver,hyperplane::HyperPlane) = satisfied(hyperplane,lshaped.x,lshaped.parameters.τ)
satisfied(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut}) = satisfied(cut,lshaped.x,lshaped.θs[cut.id],lshaped.parameters.τ)
violated(lshaped::AbstractLShapedSolver,hyperplane::HyperPlane) = !satisfied(lshaped,hyperplane)
gap(lshaped::AbstractLShapedSolver,hyperplane::HyperPlane) = gap(hyperplane,lshaped.x)
gap(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut}) = gap(cut,lshaped.x,lshaped.θs[cut.id])

function addcut!(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut},Q::Real)
    θ = lshaped.θs[cut.id]
    @unpack τ = lshaped.parameters

    if θ > -Inf && abs(θ-Q) <= τ*(1+abs(Q))
        # Optimal with respect to this subproblem
        return false
    end

    process_cut!(lshaped,cut)
    addconstr!(lshaped.mastersolver.lqmodel,lowlevel(cut)...)
    push!(lshaped.cuts,cut)
    return true
end
addcut!(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut},x::AbstractVector) = addcut!(lshaped,cut,cut(x))
addcut!(lshaped::AbstractLShapedSolver,cut::HyperPlane{OptimalityCut}) = addcut!(lshaped,cut,lshaped.x)

function addcut!(lshaped::AbstractLShapedSolver,cut::HyperPlane{FeasibilityCut})
    process_cut!(lshaped,cut)
    addconstr!(lshaped.mastersolver.lqmodel,lowlevel(cut)...)
    push!(lshaped.cuts,cut)
    return true
end

function show(io::IO, lshaped::AbstractLShapedSolver)
    println(io,typeof(lshaped).name.name)
    println(io,"State:")
    show(io,lshaped.solverdata)
    println(io,"Parameters:")
    show(io,lshaped.parameters)
end

function show(io::IO, ::MIME"text/plain", lshaped::AbstractLShapedSolver)
    show(io,lshaped)
end
# ======================================================================== #

# Plot recipe #
# ======================================================================== #
@recipe f(lshaped::AbstractLShapedSolver) = lshaped,-1
@recipe function f(lshaped::AbstractLShapedSolver, time::Real; showθ = false)
    length(lshaped.Q_history) > 0 || error("No solution data. Has solver been run?")
    showθ && (length(lshaped.θ_history) > 0 || error("No solution data. Has solver been run?"))
    Qmin = showθ ? minimum(lshaped.θ_history) : minimum(lshaped.Q_history)
    Qmax = maximum(lshaped.Q_history)
    increment = std(lshaped.Q_history)

    linewidth --> 4
    linecolor --> :black
    tickfontsize := 14
    tickfontfamily := "sans-serif"
    guidefontsize := 16
    guidefontfamily := "sans-serif"
    titlefontsize := 22
    titlefontfamily := "sans-serif"
    xlabel := time == -1 ? "Iteration" : "Time [s]"
    ylabel := "Q"
    ylims --> (Qmin-increment,Qmax+increment)
    if time == -1
        xlims --> (1,length(lshaped.Q_history)+1)
        xticks --> 1:5:length(lshaped.Q_history)
    else
        xlims --> (0,time)
        xticks --> linspace(0,time,ceil(Int,length(lshaped.Q_history)/5))
    end
    yticks --> Qmin:increment:Qmax
    xformatter := (d) -> @sprintf("%.1f",d)
    yformatter := (d) -> begin
        if abs(d) <= sqrt(eps())
            "0.0"
        elseif (log10(abs(d)) < -2.0 || log10(abs(d)) > 3.0)
            @sprintf("%.4e",d)
        elseif log10(abs(d)) > 2.0
            @sprintf("%.1f",d)
        else
            @sprintf("%.2f",d)
        end
    end

    @series begin
        label --> "Q"
        seriescolor --> :black
        if time == -1
            1:1:length(lshaped.Q_history),lshaped.Q_history
        else
            linspace(0,time,length(lshaped.Q_history)),lshaped.Q_history
        end
    end

    if showθ
        @series begin
            label --> "θ"
            linestyle --> :dash
            seriescolor --> :red
            linecolor --> :red
            linewidth --> 2
            if time == -1
                1:1:length(lshaped.θ_history),lshaped.θ_history
            else
                linspace(0,time,length(lshaped.θ_history)),lshaped.θ_history
            end
        end
    end
end
# ======================================================================== #


text:
function scaledbessel(xij, qi, qj, ti, tj, mu, eta, kappa)

dim, N = size(xij);
rcut = eta[1]
alpha = eta[2]
nbf = kappa[1]

dij = sqrt.(xij[1,:].*xij[1,:] .+ xij[2,:].*xij[2,:] .+ xij[3,:].*xij[3,:]);
dr = xij./reshape(dij,(1,N));
x =  (1.0 .- exp.(-alpha*dij/rcut))/(1-exp(-alpha));
dx = -((alpha/rcut)*exp.(-(alpha*dij/rcut)))/(exp(-alpha) - 1);

rbf = zeros(N, nbf);
drbf = zeros(dim, N, nbf);
for i = 1:nbf
    a = i*pi;
    rbf[:,i] = (sqrt(2/rcut)/i)*(sin.(a*x)./dij);
    b = (sqrt(2/rcut)/i)*(sin.(a*x)./(dij.^2) .- (a*cos.(a*x).*dx)./dij);
    drbf[:,:,i] = -dr.*reshape(b,(1,N))
    # for j = 1:3
    #     drbf[j,:,i] = (xij[i,:].*b)./dij;
    # end
end

return rbf, drbf

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
module eqjl
export OUT_put, IN_get, output_get, input_q, output_q

using DataStructures

    EQJL_ABORT = "EQJL_ABORT"

    global input_q = Deque{Any}()
    global output_q = Deque{Any}()

    function OUT_put(string_params)
        @async push!(output_q, string_params)
    end

    function IN_get()
        #here must be blocking
        while isempty(input_q)
            yield()
        end
        pop!(input_q)
    end

    function output_get()
        while isempty(output_q)
            yield()
        end
        pop!(output_q)
    end

end

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


text:
module contactForceLaw_newtons_cradle

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmatGraphics = Modia3D.Material(color="LightBlue" , transparency=0.1)    # material of Graphics
vmatSolids = Modia3D.Material(color="Red" , transparency=0.0)         # material of solids

cmat = "Steel"

Lx = 0.3
Ly = 4.1
Lz = 0.1

diameter = 1.0
lengthWire = 4.0

@assembly NewtonComparison begin
  world = Modia3D.Object3D(visualizeFrame=false)
  cradle = NewtonsCradle(world)
#  cradle2 = NewtonsCradleSolidWire(world)
end

@assembly NewtonsCradle(world) begin
  box    = Modia3D.Object3D(world, Modia3D.Box(Lx, Ly, Lz; material=vmatGraphics)   ; r=[0.0, 0.0, 0.0] )

  frame1 = Modia3D.Object3D(box; r=[0.0, -Ly/2, 0.0])
  frame2 = Modia3D.Object3D(box; r=[0.0, -Ly/4, 0.0])
  frame3 = Modia3D.Object3D(box; r=[0.0, 0.0, 0.0])
  frame4 = Modia3D.Object3D(box; r=[0.0, Ly/4, 0.0])
  frame5 = Modia3D.Object3D(box; r=[0.0, Ly/2, 0.0])

  pendulum1 = Pendulum()
  pendulum2 = Pendulum()
  pendulum3 = Pendulum()
  pendulum4 = Pendulum()
  pendulum5 = Pendulum()

  rev1   = Modia3D.Revolute(frame1, pendulum1.frame1; axis = 1 , phi_start =  -pi/3)
  rev2   = Modia3D.Revolute(frame2, pendulum2.frame1; axis = 1 )#, phi_start =  -pi/3)
  rev3   = Modia3D.Revolute(frame3, pendulum3.frame1; axis = 1 ) #, phi_start =  -pi/3)
  rev4   = Modia3D.Revolute(frame4, pendulum4.frame1; axis = 1 )#, phi_start =  pi/3 ) #
  rev5   = Modia3D.Revolute(frame5, pendulum5.frame1; axis = 1 )# , phi_start =  pi/3 ) #
end

@assembly Pendulum() begin
  wire = Modia3D.Object3D( Modia3D.Cylinder(diameter/5, lengthWire; material=vmatGraphics) )

  frame1 = Modia3D.Object3D(wire; r=[0.0, 0.0, lengthWire/2])
  frame2 = Modia3D.Object3D(wire; r=[0.0, 0.0, -lengthWire/2])

  cyl    = Modia3D.Cylinder(0.15,Lx; material=vmatGraphics)
  cylinder   = Modia3D.Object3D(frame1, cyl; R=ModiaMath.rot2(-pi/2), visualizeFrame=false)

  sphere = Modia3D.Object3D(frame2, Modia3D.Solid(Modia3D.SolidSphere(diameter) , "Steel", vmatSolids; contactMaterial = cmat) )
end

@assembly NewtonsCradleSolidWire(world) begin
  # world = Modia3D.Object3D(visualizeFrame=false)
  box    = Modia3D.Object3D(world, Modia3D.Box(Lx, Ly, Lz; material=vmatGraphics); r=[0.0, -8.0, 0.0] )

  frame1 = Modia3D.Object3D(box; r=[0.0, -Ly/2, 0.0])

  frame2 = Modia3D.Object3D(box; r=[0.0, -Ly/4, 0.0])
  frame3 = Modia3D.Object3D(box; r=[0.0, 0.0, 0.0])
  frame4 = Modia3D.Object3D(box; r=[0.0, Ly/4, 0.0])
  frame5 = Modia3D.Object3D(box; r=[0.0, Ly/2, 0.0])

  pendulum1 = Pendulum1()
  pendulum2 = Pendulum1()
  pendulum3 = Pendulum1()
  pendulum4 = Pendulum1()
  pendulum5 = Pendulum1()

  rev1   = Modia3D.Revolute(frame1, pendulum1.frame1; axis = 1 , phi_start =  -pi/3)
  rev2   = Modia3D.Revolute(frame2, pendulum2.frame1; axis = 1 , phi_start =  -pi/3)
  rev3   = Modia3D.Revolute(frame3, pendulum3.frame1; axis = 1 ) #, phi_start =  -pi/3)
  rev4   = Modia3D.Revolute(frame4, pendulum4.frame1; axis = 1 ) #, phi_start =  pi/3 ) #
  rev5   = Modia3D.Revolute(frame5, pendulum5.frame1; axis = 1  , phi_start =  pi/3 ) #
end

@assembly Pendulum1() begin
  wire = Modia3D.Object3D(Modia3D.Solid(Modia3D.SolidCylinder(diameter/5, lengthWire), "Steel", vmatSolids) )

  frame1 = Modia3D.Object3D(wire; r=[0.0, 0.0, lengthWire/2])
  frame2 = Modia3D.Object3D(wire; r=[0.0, 0.0, -lengthWire/2])

  cyl    = Modia3D.Cylinder(0.15,Lx; material=vmatGraphics)
  cylinder   = Modia3D.Object3D(frame1, cyl; R=ModiaMath.rot2(-pi/2), visualizeFrame=false)

  sphere = Modia3D.Object3D(frame2, Modia3D.Solid(Modia3D.SolidSphere(diameter) , "Steel", vmatSolids; contactMaterial = cmat) )
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1])
newton = NewtonComparison(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=false, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=false, visualizeSupportPoints=false))

model = Modia3D.SimulationModel( newton )
result = ModiaMath.simulate!(model; stopTime=10.0, tolerance=1e-8,interval=0.001, log=false)

println("... success of contactForceLaw_newtons_cradle.jl!")
end


text:
load("../../mutation", dirname(@__FILE__))

type TwoOptMutation <: Mutation
  size::Integer
  max_attempts::Integer

  TwoOptMutation(r::Representation) = new(r.size)
end

num_inputs(o::TwoOptMutation) = 1
num_outputs(o::TwoOptMutation) = 1

function operate!{T}(o::TwoOptMutation,
  inputs::Vector{IndividualStage{Vector{T}}})

  attempts = 0
  solution = input[1]
  fitness = # this is a bit of a problem.

  while (attempts < o.max_attempts)

    for i = 1:size
      for j = (i + 1):size
        swap(i, j)

      end
    end

    adjusted_fitness = ##

    ##
    if improves(adjusted_fitness, best_fitness)
      attempts = 0
      solution = adjusted
      best_fitness = adjusted_fitness
    else
      attempts += 1
    end
  end

end

register(joinpath(dirname(@__FILE__), "manifest.yml"))


text:
using CausalityTools
using CausalityToolsBase
using TimeseriesSurrogates
using PerronFrobenius
using PyCall, Conda
using CrossMappings
using HypothesisTests
using Documenter, DocumenterMarkdown
#Conda.add("scipy")
using Plots
using DynamicalSystems
using SimpleDiffEq
using DynamicalSystems
using Distributions
using StaticArrays
using Statistics
using StatsBase
using TransferEntropy
using UncertainData

ENV["GKSwstype"] = "100"

PAGES = [
    "index.md",
    "Syntax overview" => "syntax_overview.md",
    "CHANGELOG.md",
    "Causality tests" => [
        # From time series
        "causalitytests/causality_from_time_series.md",

        # From uncertain data
        "causalitytests/causality_from_uncertain_data_naive.md",
        "causalitytests/causality_from_uncertain_data_binneddatacausalitytest.md",
        "causalitytests/causality_from_uncertain_data_strictlyincreasing_interpolated.md",
        "causalitytests/causality_from_uncertain_data_naive_constrained.md",
        "causalitytests/causality_from_uncertain_data_InterpolateBinTest.md",
        "causalitytests/causality_from_uncertain_data_RandomSequencesTest.md",
        "causalitytests/causality_from_uncertain_data.md",

        # From dynamical systems
        "causalitytests/causality_from_dynamical_systems.md",

        # Tests
        "causalitytests/causality_tests.md",

        "causalitytests/ConvergentCrossMappingTest.md",
        "causalitytests/CrossMappingTest.md",
        "causalitytests/JointDistanceDistributionTest.md",
        "causalitytests/SMeasureTest.md",

        "causalitytests/TransferOperatorGridTest.md",
        "causalitytests/VisitationFrequencyTest.md",
        "causalitytests/ApproximateSimplexIntersectionTest.md",
        "causalitytests/ExactSimplexIntersectionTest.md",

        "causalitytests/PredictiveAsymmetryTest.md"
    ],
    "CausalityToolsBase" => [
        "Discretization" => "causalitytoolsbase/discretization.md",
        "Delay reconstructions" => "causalitytoolsbase/delay_reconstructions.md"
    ],
    "Transfer operator estimation" => [
        "Transfer operator" => "perronfrobenius/transferoperator.md",
        "Invariant measure" => "perronfrobenius/invariantmeasure.md"
    ],
    "Transfer entropy" => [
        "Estimators" => "transferentropy/transferentropy_estimators.md",
        "Convenience functions" => "transferentropy/convenience_functions_te.md",
        "TEVars" => "transferentropy/TEVars.md",
        "Effect of discretization scheme" => "transferentropy/examples_TE_different_partitionings.md",
    ],
    "Distance based measures" => [
        "CCM" => [
            "Overview" => "crossmappings/ccm/overview.md",
            "Cross mapping" => "crossmappings/ccm/crossmapping.md",
            "Converent cross mapping" => "crossmappings/ccm/convergentcrossmapping.md"
        ],
        "Joint distance distribution" => "algorithms/joint_distance_distribution.md"
    ],
    "Worked example" => [
        "worked_examples/worked_example_transferentropy.md"
    ],

    "Example systems" => [
        "example_systems/noise.md",
        "example_systems/example_systems_overview.md",
        "example_systems/example_systems_discrete.md",
        "example_systems/example_systems_continuous.md"
    ],

    "Surrogate data" => [
        "surrogates/iaaft_docs.md",
        "surrogates/aaft_docs.md",
        "surrogates/randomphases_docs.md",
        "surrogates/randomamplitudes_docs.md",
        "surrogates/randomshuffle_docs.md",
        "surrogates/surrogates_overview.md"
    ],
    "UncertainData" => [
        "uncertaindata/BinnedDataCausalityTest.md"
    ],
    "Tutorials" => [
        "causality" => [
            "tutorials/list_of_tutorials.md",
            "tutorials/causality/binned_uncertain_data/tutorial_BinnedDataCausalityTest_PredictiveAsymmetryTest_BinnedResampling.md",
            "tutorials/causality/binned_uncertain_data/tutorial_BinnedDataCausalityTest_PredictiveAsymmetryTest_BinnedMeanResampling.md"
        ],
        "tutorials/tutorial_SMeasureTest.md"
    ]
]

makedocs(
    sitename = "CausalityTools.jl documentation",
    #modules = [CausalityTools, TransferEntropy, PerronFrobenius, CrossMappings, CausalityToolsBase, UncertainData],
    modules = [CausalityTools, TransferEntropy, PerronFrobenius, CrossMappings, CausalityToolsBase, UncertainData],
    format = DocumenterMarkdown.Markdown(),
    linkcheck = false,
    pages = PAGES,
    highlightsig = true
)

if !Sys.iswindows()
    deploydocs(
        deps   = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4",
        "python-markdown-math", "pygments", "pymdown-extensions"),
        repo   = "github.com/kahaaga/CausalityTools.jl.git",
        target = "site",
        make = () -> run(`mkdocs build`)
    )
end


text:
# Autogenerated wrapper script for XZ_jll for aarch64-linux-musl
export liblzma, lzmadec, lzmainfo, xz, xzdec

JLLWrappers.@generate_wrapper_header("XZ")
JLLWrappers.@declare_library_product(liblzma, "liblzma.so.5")
JLLWrappers.@declare_executable_product(lzmadec)
JLLWrappers.@declare_executable_product(lzmainfo)
JLLWrappers.@declare_executable_product(xz)
JLLWrappers.@declare_executable_product(xzdec)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        liblzma,
        "lib/liblzma.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        lzmadec,
        "bin/lzmadec",
    )

    JLLWrappers.@init_executable_product(
        lzmainfo,
        "bin/lzmainfo",
    )

    JLLWrappers.@init_executable_product(
        xz,
        "bin/xz",
    )

    JLLWrappers.@init_executable_product(
        xzdec,
        "bin/xzdec",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
Helper functions to generate fake simulation data for Distributed AdaOja's algorithm
    with MPI

Author: Wayne Wang
Last modified: 08/02/2021
"""

function gen_spiked_cov_data(Σ::AbstractArray; N::Int = 1)
    """
    Generate fake data from a spiked covariance model s.t.
        X ∈ d × 1 ∼ N(0, Σ) where Σ = A diag(w)^2 A^T + σ^2 I.

    Here, A ∈ d × k is a set of k d-dimensional orthonormal vectors.
    Set w ∈ k × 1 s.t. w_i ∼ Unif(0,1) and w_1 ≥ w_1 ≥ ... ≥ w_k and w_i
    is then scaled by w_1. The σ here is a noise parameter.

    Output: X ∈ d × N
    """
    d = size(Σ, 1)
    Z = randn((d, N))

    C = cholesky(Σ)
    Z .= C.L * Z #follows N(0,Σ)

    return (N == 1) ? Z[:] : Z
end

function spiked_cov(k::Int, d::Int, σ::Real = 0.5)
    """
    The spiked covariance model: Σ = A diag(w)^2 A^T + σ^2 I.

    Here, A ∈ d × k is a set of k d-dimensional orthonormal vectors.
    Set w ∈ k × 1 s.t. w_i ∼ Unif(0,1) and w_1 ≥ w_1 ≥ ... ≥ w_k and w_i
    is then scaled by w_1. The σ here is a noise parameter.

    Output: A ∈ d × k, Σ ∈ d × d
    """
    A = randn((d, k))
    A .= Array(qr!(A).Q)

    w = [rand(Uniform(0,1)) for _ = 1:k]
    sort!(w, rev = true)
    w ./= w[1]

    Σ = A * (diagm(w).^2) * copy(A') .+ σ^2 * I(d)
    Σ .= (Σ .+ copy(Σ')) ./ 2

    return A, Σ
end

function gen_spiked_cov_data_par(num_par::Int, out_dir::AbstractString;
    N::Int = 1000, d::Int = 128*128, k::Int = 5, σ::Real = 0.5)
    """
    Generate data with spiked covariance and save into partitions
    """
    A, Σ = spiked_cov(k, d, σ)
    X = gen_spiked_cov_data(Σ, N = N)
    d_par = Int(d / num_par)
    X_par = zeros(d_par, N)
    for i = 1:num_par
        X_par .= X[d_par * (i - 1) + 1:d_par * i, :]
        data_dir = string(out_dir, "scaling_N", N, "_d", d, "_k", k, "/scaling_", num_par)
        if !isdir(data_dir)
            mkpath(data_dir)
        end
        save(string(data_dir, "/data_par_", i, ".jld2"), Dict("data_par" => X_par))
    end
    return nothing
end

function gen_data_par(num_par::Int, out_dir::AbstractString;
    N::Int = 1000, d::Int = 512*256, k::Int = 128)
    """
    Generate fake partitioned data
    """
    d_par = Int(d / num_par)
    X_par = zeros(d_par, N)
    for i = 1:num_par
        X_par .= randn((d_par, N))
        data_dir = string(out_dir, "scaling_N", N, "_d", d, "_k", k, "/scaling_", num_par)
        if !isdir(data_dir)
            mkpath(data_dir)
        end
        save(string(data_dir, "/data_par_", i, ".jld2"), Dict("data_par" => X_par))
    end
    return nothing
end

#
# Start generating data
#
using TributaryPCA

num_proc = parse(Int64, ARGS[1]) #take the cml input for number of partitions
out_dir = ARGS[2] #take the cml input for data directory
gen_data_par(num_proc, out_dir)

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
int main () {
  printInt(fact(7)) ;
  printInt(factr(7)) ;
  return 0 ;
} ;

// iteracyjnie
int fact (int n) {
  int i,r ;
  i = 1 ;
  r = 1 ;
  while (i < n+1) {
    r = r * i ;
    i++ ;
  } ;
  return r ;
} ;

// rekurencyjnie
int factr (int n) {
  if (n < 2)
    return 1 ;
  else
    return (n * factr(n-1)) ;
} ;

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
@testset "NumericalVI" begin
    b = 10
    @test NumericalVI(:quad) isa QuadratureVI
    @test NumericalVI(:mc) isa MCIntegrationVI
    @test_throws ErrorException NumericalVI(:blah)
    @test NumericalSVI(b, :quad) isa QuadratureVI
    @test NumericalSVI(b, :mc) isa MCIntegrationVI
    @test_throws ErrorException NumericalSVI(b, :blah)
    @test repr(NumericalVI(:quad)) == "Numerical Inference by Quadrature"
    @test repr(NumericalVI(:mc)) == "Numerical Inference by Monte Carlo Integration"
end


text:
@testset "Minimum spanning trees in a complete graph" begin
  @testset "Constructor with $i nodes" for i in [2, 5, 10]
    # Edge weights: high gain likelihood for the edges 1->2->3...n, low for the other edges.
    n = i
    ε = 1 / n
    graph = complete_graph(n)
    reward = Dict{Tuple{Int, Int}, Distribution}(
      (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
      for e in edges(graph))
    instance = SpanningTree(graph, reward, SpanningTreeNoSolver())

    @test instance.n_arms == n * (n - 1) / 2
    @test instance.reward == reward
  end

  @testset "State with $i nodes" for i in [2, 5, 10]
    n = i
    ε = 1 / n
    graph = complete_graph(n)
    reward = Dict{Tuple{Int, Int}, Distribution}(
      (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
      for e in edges(graph))
    instance = SpanningTree(graph, reward, SpanningTreeNoSolver())

    state = initial_state(instance)

    @test state.round == 0
    @test state.regret == 0.0
    @test state.reward == 0.0
    @test length(state.arm_counts) == ne(graph)
    @test length(state.arm_reward) == ne(graph)
    @test length(state.arm_average_reward) == ne(graph)

    for e in keys(reward)
      @test state.arm_counts[e...] == 0
      @test state.arm_reward[e...] == 0.0
      @test state.arm_average_reward[e...] == 0.0
    end
  end

  @testset "Trace with $i nodes" for i in [2, 5, 10]
    n = i
    ε = 1 / n
    graph = complete_graph(n)
    reward = Dict{Tuple{Int, Int}, Distribution}(
      (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
      for e in edges(graph))
    instance = SpanningTree(graph, reward, SpanningTreeNoSolver())

    trace = initial_trace(instance)

    @test length(trace.states) == 0
    @test length(trace.arms) == 0
    @test length(trace.reward) == 0
    @test length(trace.policy_details) == 0
    @test length(trace.time_choose_action) == 0

    @test eltype(trace.states) == State{Tuple{Int, Int}}
    @test eltype(trace.arms) == Vector{Tuple{Int, Int}}
    @test eltype(trace.reward) == Vector{Float64}
    @test eltype(trace.time_choose_action) == Int
  end

  @testset "Pull with $i nodes" for i in [2, 5, 10]
    n = i
    graph = complete_graph(n)
    reward = Dict{Tuple{Int, Int}, Distribution}(
      (src(e), dst(e)) => Bernoulli((src(e) == 1 && dst(e) == n) ? 1 : 0)
      for e in edges(graph))
    instance = SpanningTree(graph, reward, SpanningTreeLightGraphsPrimSolver())

    Random.seed!(1)
    @test pull(instance, [(1, i)]) == ([1.0], 0.0)
  end

  @testset "Check feasibility with 3 nodes" begin
    n = 3
    ε = 1 / n
    graph = complete_graph(n)
    reward = Dict{Tuple{Int, Int}, Distribution}(
      (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
      for e in edges(graph))
    instance = SpanningTree(graph, reward, SpanningTreeNoSolver())

    @test ! is_feasible(instance, [(1, 3)])
    @test is_feasible(instance, [(1, 2), (2, 3)])
    @test is_feasible(instance, [(1, 3), (1, 3)])
    @test is_feasible(instance, [(2, 1), (2, 3)])
    @test ! is_feasible(instance, Tuple{Int, Int}[])
    @test ! is_feasible(instance, [(1, 2), (2, 3), (3, 1)])
  end

  @testset "Check partial acceptability with 3 nodes" begin
    n = 3
    ε = 1 / n
    graph = complete_graph(n)
    reward = Dict{Tuple{Int, Int}, Distribution}(
      (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
      for e in edges(graph))
    instance = SpanningTree(graph, reward, SpanningTreeNoSolver())

    @test is_partially_acceptable(instance, [(1, 3)])
    @test is_partially_acceptable(instance, [(1, 2), (2, 3)])
    @test is_partially_acceptable(instance, [(1, 3), (1, 2)])
    @test is_partially_acceptable(instance, [(2, 1), (2, 3)])
    @test is_partially_acceptable(instance, Tuple{Int, Int}[])
    @test ! is_partially_acceptable(instance, [(1, 3), (1, 3)])
    @test ! is_partially_acceptable(instance, [(1, 2), (2, 3), (3, 1)])
    @test ! is_partially_acceptable(instance, [(1, 2), (2, 1)])
  end

  @testset "LightGraphs.jl Prim solver" begin
    @testset "Constructor" for i in [2, 3]#, 5, 10]
      n = i
      ε = 1 / n
      graph = complete_graph(n)
      reward = Dict{Tuple{Int, Int}, Distribution}(
        (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
        for e in edges(graph))
      instance = SpanningTree(graph, reward, SpanningTreeLightGraphsPrimSolver())

      @test instance.solver != nothing
    end

    @testset "Solve with $i nodes" for i in [2, 5, 10]
      n = i
      ε = 1 / (n + 1)
      graph = complete_graph(n)
      reward = Dict{Tuple{Int, Int}, Distribution}(
        (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
        for e in edges(graph))
      instance = SpanningTree(graph, reward, SpanningTreeLightGraphsPrimSolver())

      Random.seed!(i)
      drawn = Dict(k => rand() for (k, _) in reward)
      solution = solve_linear(instance, drawn)
      @test is_feasible(instance, solution)
    end
  end

  if ! is_travis
    @testset "LP solver" begin
      @testset "Constructor" for i in [2, 5, 10]
        n = i
        ε = 1 / n
        graph = complete_graph(n)
        reward = Dict{Tuple{Int, Int}, Distribution}(
          (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
          for e in edges(graph))
        instance = SpanningTree(graph, reward, SpanningTreeLPSolver(Gurobi.Optimizer))

        @test instance.solver != nothing
      end

      @testset "Solve with $i nodes" for i in [2, 5, 10]
        n = i
        ε = 1 / (n + 1)
        graph = complete_graph(n)
        reward = Dict{Tuple{Int, Int}, Distribution}(
          (src(e), dst(e)) => Bernoulli((src(e) + 1 == dst(e)) ? (1 - ε) : (ε))
          for e in edges(graph))
        instance = SpanningTree(graph, reward, SpanningTreeLPSolver(Gurobi.Optimizer))

        Random.seed!(i)
        drawn = Dict(k => rand() for (k, _) in reward)
        solution = solve_linear(instance, drawn)
        @test is_feasible(instance, solution)
      end
    end
  end
end


text:
import GreyDecision.GreyNumbers: GreyNumber
import GreyDecision.Utility: makeminmax
import GreyDecision.Electre: electre

@testset "Electre with white numbers" begin
    tol = 0.00001
        w = [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077]
        Amat = [
            100 92 10 2 80 70 95 80
            80 70 8 4 100 80 80 90
            90 85 5 0 75 95 70 70
            70 88 20 18 60 90 95 85
        ]
        greyamat = Array{GreyNumber{Float64}, 2}(undef, 4, 8)
        for i in 1:4
            for j in 1:8
                greyamat[i,j] = GreyNumber{Float64}(Float64(Amat[i,j]), Float64(Amat[i, j]))
            end
        end
        dmat = greyamat
        fns = makeminmax([maximum for i = 1:8])
        result = electre(dmat, w, fns)

        @test isa(result, Dict)
        @test isa(result[:bestindex], Tuple)
        @test result[:bestindex][1] == 4

        @test isapprox(result[:C][1], GreyNumber(0.36936937,0.36936937), atol = tol)
        @test isapprox(result[:C][2], GreyNumber(0.01501502, 0.01501502),  atol = tol)
        @test isapprox(result[:C][3], GreyNumber(-2.47347347, -2.47347347),  atol = tol)
        @test isapprox(result[:C][4], GreyNumber(2.08908909, 2.08908909), atol = tol)

        @test isapprox(result[:D][1], GreyNumber(0.1914244, 0.1914244), atol = tol)
        @test isapprox(result[:D][2], GreyNumber(-0.1903929,-0.1903929),  atol = tol)
        @test isapprox(result[:D][3], GreyNumber(2.8843076, 2.8843076),  atol = tol)
        @test isapprox(result[:D][4], GreyNumber(-2.8853391, -2.8853391), atol = tol)
end


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    x = y = sqrt.(1:10:101)
    z = abs.(randn(10,10))
    fig = Figure(resolution = (450, 400))
    ax = Axis(fig, aspect = 1, xlabel = "x", ylabel = "y")
    hmap = heatmap!(x,y,z, colormap = :magma)
    cbar = Colorbar(fig, hmap, label = "z values", width = 15, ticksize=15,
                tickalign = 1, height = Relative(0.96))
    fig[1, 1] = ax
    fig[1, 2] = cbar
    colgap!(fig.layout, 7) # HIDE
    fig
    save(joinpath(@__DIR__, "output", "heatmapIrregular.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module MNPDynamics

using LinearAlgebra
using SparseArrays
using OrdinaryDiffEq
using SparseDiffTools
#using Sundials

include("Neel.jl")
include("Brown.jl")

end # module


text:
####################################################
# Metropolis–Hasting Step
####################################################

function MHStep!(CPM::CellPotts)

    #unpack current step structure update
    stepInfo = CPM.stepInfo #does this allocate/slow things down? Parameters.jl ?

    #Loop through until a good source target is found
    searching = true
    while searching
        #Pick a random location on the graph
        stepInfo.sourceNode = rand(1:nv(CPM.graph.network))
        #What cell does it belong to?
        stepInfo.sourceCell = CPM.graph.σ[stepInfo.sourceNode]

        #Get all of the unique cell IDs neighboring this Node
        stepInfo.sourceNodeNeighbors = neighbors(CPM.graph.network, stepInfo.sourceNode)
        stepInfo.possibleCellTargets = unique(CPM.graph.σ[stepInfo.sourceNodeNeighbors])
        #Choose a target
        stepInfo.targetCell = rand(stepInfo.possibleCellTargets)

        #Some checks before attempting a flip
            #In the middle of a cell
            inMiddle = all(x -> x == stepInfo.sourceCell, stepInfo.possibleCellTargets)
            #will fragment the cell
            isArticulation = CPM.graph.isArticulation[stepInfo.sourceNode]
            #target is the same as source cell
            isSource = stepInfo.targetCell == stepInfo.sourceCell

        #if all checks pass, attempt flip
        if !(inMiddle | isArticulation | isSource)
            searching = false
        else
            CPM.stepCounter += 1
        end
    end

    #Calculate the change in energy when source node is modified
    ΔH =  sum([f(CPM, stepInfo) for f in CPM.M.penalties])

    #Calculate an acceptance ratio
    acceptRatio = min(1.0,exp(-ΔH/CPM.M.temperature))

    if rand() < acceptRatio #If the acceptance ratio is large enough

        #Need to update all cell and graph properties
        #---Cell properties---

        #Update cell volumes
        CPM.cell.volumes[stepInfo.sourceCell] -= 1
        CPM.cell.volumes[stepInfo.targetCell] += 1

        #---Graph properties---

        #Cell IDs
        CPM.graph.σ[stepInfo.sourceNode] = stepInfo.targetCell
        CPM.graph.τ[stepInfo.sourceNode] = CPM.cell.types[stepInfo.targetCell]

        #Articulation points
        UpdateConnections!(CPM.graph)

        #---Overall properties---
        #Update energy
        CPM.energy += ΔH
        #Update visual
        CPM.visual[stepInfo.sourceNode] = stepInfo.targetCell

    end

    return nothing
end

function MHStep_naive!(CPM::CellPotts)

    #Create a structure to hold the source and target information
    stepInfo = MHStepInfo() #should add as input, might save a little time

    stepInfo.sourceNode = rand(1:nv(CPM.graph.network))
    stepInfo.sourceCell = CPM.graph.σ[stepInfo.sourceNode]
    stepInfo.sourceNodeNeighbors = neighbors(CPM.graph.network, stepInfo.sourceNode)
    stepInfo.targetCell = rand(CPM.graph.σ[stepInfo.sourceNodeNeighbors])

    #Calculate the change in energy when source node is modified
    ΔH =  sum([f(CPM, stepInfo) for f in CPM.M.penalties])

    #Calculate an acceptance ratio
    acceptRatio = min(1.0,exp(-ΔH/CPM.M.temperature))

    if rand() < acceptRatio #If the acceptance ratio is large enough

        #Need to update all cell and graph properties
        #---Cell properties---

        #Update cell volumes
        CPM.cell.volumes[stepInfo.sourceCell] -= 1
        CPM.cell.volumes[stepInfo.targetCell] += 1

        #---Graph properties---

        #Cell IDs
        CPM.graph.σ[stepInfo.sourceNode] = stepInfo.targetCell
        CPM.graph.τ[stepInfo.sourceNode] = CPM.cell.types[stepInfo.targetCell]

        #---Overall properties---
        #Update energy
        CPM.energy += ΔH
        #Update visual
        CPM.visual[stepInfo.sourceNode] = stepInfo.targetCell

    end
    return nothing
end