text:
mutable struct RandomPositionsXXZWithXField <: Model
    posdatadescriptor
    interaction
    fields
    scale_interactions
    basis
end

# small constructor - useful for loading
RandomPositionsXXZWithXField(posdatadescriptor, interaction) = RandomPositionsXXZWithXField(posdatadescriptor, interaction, nothing, nothing, nothing)

function Base.getproperty(xxz::RandomPositionsXXZWithXField, p::Symbol)
    if hasfield(typeof(xxz), p)
        return getfield(xxz, p)
    else
        return getproperty(getfield(xxz, :posdatadescriptor), p)
    end
end

model_name(::RandomPositionsXXZWithXField) = "Jxxz_with_field"
model_fileprefix(xxz::RandomPositionsXXZWithXField) = @sprintf("%s_%id_alpha_%.1f_N_%02i", xxz.geometry, xxz.dimension, xxz.interaction.α, xxz.system_size)

struct PreparedJXXZWithField{JT,SB}
    J::JT
    fields::Vector{Float64}
    basis::SB
end

# shot, field, rho
parameter_dims(model::PreparedJXXZWithField) = (size(model.J)[3], length(model.fields), size(model.J)[4])

function initialize_model!(model::RandomPositionsXXZWithXField, array_initializer)
    posdata = load(model.posdatadescriptor)
    shots = posdata.shots
    N = posdata.system_size
    ρs = posdata.ρs
    J = fill!(array_initializer(Float64,N,N,shots, length(ρs)), 0)
    _compute_interactions!(J, posdata, model.interaction, model.scale_interactions)

    return PreparedJXXZWithField(J, model.fields, model.basis)
end

function _ensemble_J_mean(interaction, geometry, positions)
    mapreduce(pos->sum(interaction_matrix(interaction, geometry, pos)), +, eachslice(positions; dims=3))/size(positions,2)/size(positions,3)
end

function _compute_interactions!(arr, posdata, interaction, scaling)
    logmsg("Calculating interactions (scaling = $(scaling))")
    N = size(arr, 1)
    nshots = size(arr, 3)
    for (i, ρ) in enumerate(posdata.ρs)
        geom = geometry_from_density(posdata.geometry, ρ, N, posdata.dimension)
        for shot in 1:nshots
            arr[:,:, shot, i] .= interaction_matrix(interaction, geom, posdata[:,:,shot,i])
            if scaling == :shot
                arr[:,:, shot, i] ./= sum(view(arr, :,:, shot, i))/N
            end
        end
        if scaling == :ensemble
            factor = sum(@view arr[:,:,:,i])/nshots/N
            logmsg("rho = $ρ -> J_mean = $factor")
            arr[:,:,:,i] ./= factor
        end
    end
end

function _workload_chunks(total, nchunks)
    splits = [round(Int, s) for s in range(0; stop=total, length=nchunks+1)]
    [(splits[i]+1):splits[i+1] for i in 1:nchunks]
end

function split_workload(model::PreparedJXXZWithField, parts)
    nshots = size(model.J, 3)
    nρs = size(model.J, 4)
    return _workload_chunks(nshots*nρs, parts)
end

function do_parameters(diag_callback, model::PreparedJXXZWithField, parameter_chunk)
    N = size(model.J, 1)
    spin_ops = real.(symmetrize_operator.(op_list(σx/2, N), Ref(model.basis)))
    field_operator = sum(spin_ops)
    nshots = size(model.J, 3)
    for index in parameter_chunk
        logmsg("Starting index $(index) of $(parameter_chunk)")
        ρindex, shot = _flat_to_indices(index, nshots)

        J = @view model.J[:,:, shot, ρindex]
        H_int = symmetrize_operator(xxzmodel(J, -0.73), model.basis)

        for (k, h) in enumerate(model.fields)
            logmsg("Doing #rho=$ρindex #shot=$shot #field=$k")
            parameterI = CartesianIndex(shot, k, ρindex)
            diag_callback(parameterI, H_int + h*field_operator)
        end
    end
end


text:
# NOTE: this model file is run from `/deps/generate_models.jl` or in the build phase.

# Model as in https://www.sas.upenn.edu/~jesusfv/ARE_Estimation.pdf

function FVGQ20()
    # Parameter values
    # δ = 0.025
    # ε = 10
    # ϕ = 0
    # γ2 = 0.001

    # Estimated parameters, taken from FV(2010), Table 3, p. 38, median estimate parameters
    # β = 0.998
    # h = 0.97
    # ϑ = 1.17
    # κ = 9.51
    # α = 0.21

    # θp = 0.82
    # χ = 0.63
    # γR = 0.77
    # γy = 0.19

    # γΠ = 1.29
    # Πbar = 1.01
    # ρd  = 0.12
    # ρφ = 0.93
    # σ_A = -3.97
    # σ_d = -1.51

    # σ_φ = -2.36
    # σ_μ = -5.43
    # σ_m = -5.85
    # Λμ = 3.4e-3
    # ΛA = 2.8e-3

    # 7 variable relabelings
    # 7 observations
    # 38 equations

    ∞ = Inf
    # 28 parameters. 23 of them to estimate, 5 parameters as constants
    @variables δ, ε, ϕ, γ2, β, h, ϑ, κ, α, θp, χ, γR, γy, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA, Ω_ii
    # States and Controls, 14 + 24 = 38
    @variables t::Integer, c_m(..), x_m(..), Π_m(..), w_m(..), R_m(..), yd_m(..), vp_m(..), k(..), d(..), φ(..), μ_I(..), μ_A(..), mshock(..), g(..), c(..), x(..), Π(..), w(..), R(..), yd(..), vp(..), Πstar(..), u(..), q(..), λ(..), μ_z(..), r(..), l(..), mc(..), g1(..), g2(..), ob_μz(..), ob_Π(..), ob_R(..), ob_dw(..), ob_dy(..), ob_l(..), ob_μ(..)

    # Precompute some parameters, and steady state values
    ΛYd = (ΛA + α * Λμ) / (1 - α)
    Λx = exp(ΛYd)
    μ_z_val = exp(ΛYd)
    μ_I_val = exp(Λμ)
    γ1 = μ_z_val * μ_I_val / β - (1 - δ)
    Rbar = Πbar * μ_z_val / β
    Πstar_val = ((1 - θp * Πbar^(-(1 - ε) * (1 - χ))) / (1 - θp))^(1 / (1 - ε))
    vp_val = (1 - θp) / (1 - θp * Πbar^((1 - χ) * ε)) * Πstar_val^(-ε)
    mc_val = (1 - β * θp * Πbar^((1 - χ) * ε)) / (1 - β * θp * Πbar^(-(1 - ε) * (1 - χ))) * (ε - 1) / ε * Πstar_val
    w_val = (1 - α) * (mc_val * (α / γ1)^α)^(1 / (1 - α))
    k_val = α / (1 - α) * w_val / γ1 * μ_z_val * μ_I_val
    yd_val = (exp(ΛA) / μ_z_val * k_val^α - ϕ) / vp_val
    x_val = (μ_z_val * μ_I_val - (1 - δ)) / (μ_z_val * μ_I_val) * k_val
    c_val = (1 - g_bar) * yd_val - x_val
    λ_val = (1 - h * β / μ_z_val) / (1 - h / μ_z_val) / c_val
    ψ = w_val * λ_val
    g1_val = λ_val * yd_val / (1 - β * θp * Πbar^(-(1 - ε) * (1 - χ))) * (ε - 1) / ε * Πstar_val
    g2_val = λ_val * yd_val / (1 - β * θp * Πbar^(-(1 - ε) * (1 - χ))) * Πstar_val

    x_sym = [c_m, x_m, Π_m, w_m, R_m, yd_m, vp_m, k, d, φ, μ_I, μ_A, mshock, g]
    y_sym = [c, x, Π, w, R, yd, vp, Πstar, u, q, λ, μ_z, r, l, mc, g1, g2, ob_μz, ob_Π, ob_R, ob_dw, ob_dy, ob_l, ob_μ]
    # Merge p and p_f
    p = [β, h, ϑ, κ, α, θp, χ, γR, γy, γΠ, Πbar, ρd, ρφ, ρg, g_bar, σ_A, σ_d, σ_φ, σ_μ, σ_m, σ_g, Λμ, ΛA, δ, ε, ϕ, γ2, Ω_ii]

    H = [
         # Household Problem
         d(t) * (c(t) - h * c_m(t) * μ_z(t)^(-1))^(-1) - h * β * d(t+1) * (c(t+1) * μ_z(t+1) - h * c(t))^(-1) - λ(t),
         ψ * φ(t) * l(t)^ϑ - w(t) * λ(t),
         λ(t) - β * λ(t+1) * μ_z(t+1)^(-1) * R(t) / Π(t+1),
         γ1 + γ2 * (u(t) - 1) - r(t),
         q(t) - β * λ(t+1) / λ(t) / μ_z(t+1) / μ_I(t+1) * ((1 - δ) * q(t+1) + r(t+1) * u(t+1) - (γ1 * (u(t+1) - 1) + γ2 / 2 * (u(t+1) - 1)^2)),
         q(t) * (1 - (κ / 2 * (x(t) / x_m(t) * μ_z(t) - Λx)^2) - (κ * (x(t) / x_m(t) * μ_z(t) - Λx) * x(t) / x_m(t) * μ_z(t))) + β * q(t+1) * λ(t+1) / λ(t) / μ_z(t+1) * κ * (x(t+1) / x(t) * μ_z(t+1) - Λx) * (x(t+1) / x(t) * μ_z(t+1))^2 - 1,

         # Firm Problem
         λ(t) * mc(t) * yd(t) + β * θp * (Π(t)^χ / Π(t+1))^(-ε) * g1(t+1) - g1(t),
         λ(t) * Πstar(t) * yd(t) + β * θp * (Π(t)^χ / Π(t+1))^(1 - ε) * Πstar(t) / Πstar(t+1) * g2(t+1) - g2(t),
         ε * g1(t) - (ε - 1) * g2(t),
         u(t) * k(t) / l(t) - α / (1 - α) * w(t) / r(t) * μ_z(t) * μ_I(t),
         (1 / (1 - α))^(1 - α) * (1 / α)^α * w(t)^(1 - α) * r(t)^α - mc(t),

         # Price Evolution
         θp * (Π_m(t)^χ / Π(t))^(1 - ε) + (1 - θp) * Πstar(t)^(1 - ε) - 1,

         # Taylor Rule
         R(t) / Rbar - (R_m(t) / Rbar)^γR * ((Π(t) / Πbar)^γΠ * ((yd(t) / yd_m(t) * μ_z(t)) / exp(ΛYd))^γy)^(1 - γR) * mshock(t),

         # Market Clear
         c(t) + x(t) + g(t) + μ_z(t)^(-1) * μ_I(t)^(-1) * (γ1 * (u(t) - 1) + γ2 / 2 * (u(t) - 1)^2) * k(t) - yd(t),
         (μ_A(t) / μ_z(t) * (u(t) * k(t))^α * l(t)^(1 - α) - ϕ) / vp(t) - yd(t),
         θp * (Π_m(t)^χ / Π(t))^(-ε) * vp_m(t) + (1 - θp) * Πstar(t)^(-ε) - vp(t),
         k(t+1) * μ_z(t) * μ_I(t) - (1 - δ) * k(t) - μ_z(t) * μ_I(t) * (1 - κ / 2 * (x(t) / x_m(t) * μ_z(t) - Λx)^2) * x(t),

         # Variable Mapping
         c(t) - c_m(t+1), x(t) - x_m(t+1), Π(t) - Π_m(t+1), w(t) - w_m(t+1), R(t) - R_m(t+1), yd(t) - yd_m(t+1), vp(t) - vp_m(t+1),

         # Shock Evolution
         μ_z(t) - μ_A(t)^(1 / (1 - α)) * μ_I(t)^(α / (1 - α)),
         log(d(t+1)) - ρd * log(d(t)),
         log(φ(t+1)) - ρφ * log(φ(t)),
         log(μ_I(t+1)) - Λμ,
         log(μ_A(t+1)) - ΛA,
         log(mshock(t+1)) - 0,
         log(g(t+1)) - ρg * log(g(t)) - (1 - ρg) * log(g_bar * yd_val),

         # Observation equations Π, R, dw, dy, l, μ_I(-1)
         ob_μz(t) - log(μ_z(t)), ob_Π(t) - log(Π(t)), ob_R(t) - log(R(t)), ob_dw(t) - (log(w(t)) - log(w_m(t))), ob_dy(t) - (log(yd(t)) - log(yd_m(t))), ob_l(t) - log(l(t)), ob_μ(t) - log(μ_I(t))
    ]

    steady_states = [c_m(∞) ~ c_val, x_m(∞) ~ x_val, Π_m(∞) ~ Πbar, w_m(∞) ~ w_val, R_m(∞) ~ Rbar, yd_m(∞) ~ yd_val, vp_m(∞) ~ vp_val, k(∞) ~ k_val, d(∞) ~ 1.0, φ(∞) ~ 1.0, μ_I(∞) ~ exp(Λμ), μ_A(∞) ~ exp(ΛA), mshock(∞) ~ 1.0, g(∞) ~ g_bar * yd_val, c(∞) ~ c_val, x(∞) ~ x_val, Π(∞) ~ Πbar, w(∞) ~ w_val, R(∞) ~ Rbar, yd(∞) ~ yd_val, vp(∞) ~ vp_val, Πstar(∞) ~ Πstar_val, u(∞) ~ 1.0, q(∞) ~ 1.0, λ(∞) ~ λ_val, μ_z(∞) ~ exp(ΛYd), r(∞) ~ γ1, l(∞) ~ 1.0, mc(∞) ~ mc_val, g1(∞) ~ g1_val, g2(∞) ~ g2_val, ob_μz(∞) ~ ΛYd, ob_Π(∞) ~ log(Πbar),ob_R(∞) ~ log(Rbar), ob_dw(∞) ~ 0.0, ob_dy(∞) ~ 0.0, ob_l(∞) ~ 0.0, ob_μ(∞) ~ Λμ]

    n_ϵ = 6
    n_x = length(x_sym)
    n_y = length(y_sym)
    n_z = 6
    Γ = zeros(Num, n_ϵ, n_ϵ) # Make sure it is not a float64 matrix
    Γ[1, 1] = σ_d
    Γ[2, 2] = σ_φ
    Γ[3, 3] = σ_μ
    Γ[4, 4] = σ_A
    Γ[5, 5] = σ_m
    Γ[6, 6] = σ_g
    η_mat = zeros(n_x, n_ϵ)
    # The 9-14 states are: d φ μ_I μ_A mshock g
    η_mat[9, 1] = 1
    η_mat[10, 2] = 1
    η_mat[11, 3] = 1
    η_mat[12, 4] = 1
    η_mat[13, 5] = 1
    η_mat[14, 6] = 1
    # Define Q
    Q = zeros(n_z, n_x + n_y)
    Q[:, (n_y - 5):n_y] = [1.0 0.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 0.0 0.0 1.0]
    # Define Ω
    Ω = fill(Ω_ii, n_z)

	return H, (; t, x = x_sym, y = y_sym, p, steady_states, Γ, η = η_mat, Q, Ω), "FVGQ20"

end

ext: jl


text:
struct A1 <: Layout
    a::Bool
    b::Cint
    c::Float64
    d::Cstring
end

@testset "field access" begin
    a = fill(UInt64(0), 3)
    p = pointer_from_vector(a)
    cs = CStruct{A1}(p)
    @test cs.a == 0
    @test cs.b == 0
    @test cs.c == 0
    @test cs.d == ""
    v1 = true
    v2 = 0x12345678
    v3 = 47.11
    v4 = "hallo"
    cs.a = v1
    cs.b = v2
    cs.c = v3
    cs.d = v4
    @test cs.a == v1
    @test cs.b == v2
    @test cs.c == v3
    @test cs.d == v4
end

@testset "index access" begin
    a = fill(UInt64(0), 100)
    p = pointer_from_vector(a)
    cv = CVector{Int}(p, 3)
    @test length(cv) == 3
    cv[1:3] .= (1, 2, 3)
    @test cv[2] == 2
    @test cv[[1,3]] == [1, 3]
end

struct A2 <: Layout
    a::Ptr{A2}
end

@testset "self-referencing" begin
    a = fill(UInt8(0), 100)
    p = pointer_from_vector(a)
    cs = CStruct{A2}(p)
    @test cs.a === nothing
    io = IOBuffer()
    show(io, cs)
    @test String(take!(io)) == "CStruct{A2}(nothing)"
    cs.a = cs
    @test cs.a === cs
    show(io, cs)
    @test String(take!(io)) == "CStruct{A2}(CStruct{A2}(#= circular reference @-1 =#))"
end

struct A3 <: Layout
    len::Int
    vec::LVarVector{Float64, (x) -> x.len}
end

@testset "variable vector at end of struct" begin
    a = fill(Int(0), 1024)
    p = pointer_from_vector(a)
    LEN = 25
    cs = CStruct{A3}(p)
    cs.len = LEN
    @test cs.vec isa CVector{Float64}
    @test length(cs.vec) == cs.len == LEN
end

struct A4 <: Layout
    len::Int
    vec::Ptr{LVarVector{Float64, (x) -> x.len}}
end

@testset "pointer to variable vector" begin
    a = fill(Int(0), 1024)
    p = pointer_from_vector(a)
    a[2] = p + 32
    LEN = 25
    cs = CStruct{A4}(p)
    cs.len = LEN
    @test cs.vec isa CVector{Float64}
    @test length(cs.vec) == cs.len == LEN
end


text:
mutable struct TestDev <: AbstractDAQ
    devname::String
    nchans::Int
    channames::Vector{String}
    reading::Bool
    rate::Float64
    nsamples::Int
    freq::Float64
    noise::Float64
    t1::UInt64
    E::Matrix{Float64}
    tsk::Task
    time::DateTime
    chanidx::Dict{String,Int}
end

"""
`TestDev(devname, nchans; channames="E")`

Creates a test device, useful for testing stuff.
"""
function TestDev(devname, nchans; channames="E")

    nz = ceil(Int, log10(nchans + 1000*eps(Float64(nchans))))
    nz = max(nz,1)
    if isa(channames, AbstractString) || isa(channames, Symbol)
        chn = string(channames) .* numstring.(1:nchans, nz)
    elseif length(channames != nchans)
        throw(DomainError(nchans, "`channames` should have length 1 or equal to `nchans`"))
    else
        chn = string.(channames)
    end

    chanidx = Dict{String,Int}()
    for (i,c) in enumerate(chn)
        chanidx[c] = i
    end

    return TestDev(devname, nchans, chn, false, 1000.0,
                   1, 1.0, 0.01, 0, zeros(0,0), Task(_->1), now(), chanidx)
end

numchannels(dev::TestDev) = dev.nchans
daqchannels(dev::TestDev) = dev.channames

function daqaddinput(dev::TestDev, chans; names="E")
    nchans = length(chans)

    nz = ceil(Int, log10(nchans + 1000*eps(Float64(nchans))))
    nz = max(nz,1)
    if isa(names, AbstractString) || isa(names, Symbol)
        chn = string(names) .* numstring.(1:nchans, nz)
    elseif length(names != nchans)
        throw(DomainError(nchans, "`names` should have length 1 or equal to `nchans`"))
    else
        chn = string.(names)
    end

    dev.nchans = nchans
    dev.channames = chn
    chanidx = Dict{String,Int}()
    for (i,c) in enumerate(chn)
        chanidx[c] = i
    end
    dev.chanidx = chanidx
    return
end

function daqconfig(dev::TestDev; kw...)

    if haskey(kw, :rate) && haskey(kw, :dt)
        error("Parameters `rate` and `dt` can not be specified simultaneously!")
    elseif haskey(kw, :rate) || haskey(kw, :dt)
        if haskey(kw, :rate)
            rate = kw[:rate]
        else
            dt = kw[:dt]
            rate = 1.0 / dt
        end
    else
        error("Either `rate` or `dt` should be specified!")
    end

    if haskey(kw, :nsamples) && haskey(kw, :time)
        error("Parameters `nsamples` and `time` can not be specified simultaneously!")
    elseif haskey(kw, :nsamples) || haskey(kw, :time)
        if haskey(kw, :nsamples)
            nsamples = kw[:nsamples]
        else
            tt = kw[:time]
            nsamples = round(Int, tt * rate)
        end
    else
        error("Either `nsamples` or `time` should be specified")
    end

    dev.rate = rate
    dev.nsamples = nsamples

    # Allocate the buffer:
    dev.E = zeros(dev.nchans, nsamples)
    return
end

function filldata!(dev)
    nchans = dev.nchans
    nt = dev.nsamples

    dev.E .= randn(nchans, nt) .* dev.noise
    f = dev.freq
    fs = dev.rate
    for i in 1:nt
        x = sin(2π*f * (i-1)/fs)
        for k in 1:nchans
            dev.E[k,i] += x
        end
    end

end

function readtestsamples(dev::TestDev)

    dev.reading = true
    ttot = dev.nsamples / dev.rate
    dev.t1 = time_ns()
    sleep(ttot)
    filldata!(dev)
    dev.reading = false

end

function daqacquire(dev::TestDev)
    dev.time = now()
    readtestsamples(dev)
    return MeasData{Matrix{Float64}}(dev.devname, "TestDev", dev.time, dev.rate, dev.E,
                                     dev.chanidx)
end

function daqstart(dev::TestDev, usethread=true)
    dev.time = now()
    t =  @async readtestsamples(dev)
    dev.tsk = t
    return t
end

function daqread(dev::TestDev)

    #if !istaskdone(dev.tsk) && istaskstarted(dev.tsk)
    #    wait(dev.tsk)
    #end
    wait(dev.tsk)
    #println("waiting")

    return MeasData{Matrix{Float64}}(dev.devname, "TestDev", dev.time, dev.rate, dev.E,
                                     dev.chanidx)
end

isreading(dev::TestDev) = dev.reading
samplesread(dev::TestDev) =
    min(dev.nsamples, round(Int, (time_ns()-dev.t1)*1e-9 * dev.rate))
isdaqfinished(dev::TestDev) = !dev.reading


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_deps = joinpath(_pkg_root,"deps")
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

download("https://api.plot.ly/v2/plot-schema?sha1", joinpath(_pkg_deps,"plotschema.json"))
download("https://cdn.plot.ly/plotly-latest.min.js", joinpath(_pkg_assets,"plotly-latest.min.js"))


text:
# Public interface for AbstractConstraints

"""
    constraints(c::AbstractConstraints, f::AbstractObjective, x)

Return a values of constraints for a variable `x` given the set of constraints object `c`.
"""
value(c::AbstractConstraints, x) = nothing

"""
    penalty(constraints, x)

Calculate a penalty for the variable `x` given the set of `constraints`.
"""
penalty(c::AbstractConstraints, x) = 0

"""
    penalty!(fitness, constraints, population)

Set a penalty to the `fitness` given the set of `constraints` and `population`.
"""
penalty!(fitness, c::AbstractConstraints, population) = fitness

"""
    clip!(c::AbstractConstraints, x)

Clip a variable `x` values to satisfy the bound constraints `c`.
"""
apply!(c::AbstractConstraints, x) = x

# Auxillary functions

"""
    isfeasible(constraints, x) -> Bool
    isfeasible(bounds, x) -> Bool

Return `true` if point `x` is feasible, given the `constraints` which
specify bounds `lx`, `ux`, `lc`, and `uc`. `x` is feasible if

    lx[i] <= x[i] <= ux[i]
    lc[i] <= c(x)[i] <= uc[i]

for all possible `i`.
"""
function isfeasible(bounds::ConstraintBounds, x, c::Union{Vector,Nothing}=nothing)
    isf = true
    for (i,j) in enumerate(bounds.eqx)
        isf &= x[j] == bounds.valx[i]
    end
    for (i,j) in enumerate(bounds.ineqx)
        isf &= bounds.σx[i]*(x[j] - bounds.bx[i]) >= 0
    end
    if c !== nothing
        for (i,j) in enumerate(bounds.eqc)
            isf &= c[j] == bounds.valc[i]
        end
        for (i,j) in enumerate(bounds.ineqc)
            isf &= bounds.σc[i]*(c[j] - bounds.bc[i]) >= 0
        end
    end
    isf
end
isfeasible(c::AbstractConstraints, x) = isfeasible(c.bounds, x, value(c, x))

# Implementations

"""Type for an empty set of constratins"""
struct NoConstraints <: AbstractConstraints end
isfeasible(c::NoConstraints, x)  = true
getproperty(c::NoConstraints, s::Symbol) =
    s == :bounds ? ConstraintBounds(Float64[], Float64[], Float64[], Float64[]) : nothing

"""Box constraints"""
struct BoxConstraints{T} <: AbstractConstraints
    bounds::ConstraintBounds{T}
end
BoxConstraints(lower::AbstractVector{T}, upper::AbstractVector{T}) where {T} =
    BoxConstraints(ConstraintBounds(lower, upper, [], []))
BoxConstraints(lower::T, upper::T, size) where {T} =
    BoxConstraints(fill(lower, size), fill(upper, size))
apply!(c::BoxConstraints, x) = clip!(c.bounds, x)

"""
This type encodes constraints as a following additional penalty for an objective function:

``p(x) = \\sum^n_{i=1} r_i max(0, g_i(x))^2``

where ``g_i(x)`` is an inequality constraint. The equality constraints ``h_i(x) = 0`` are transformed
to inequality constraints as follows:

``h(x) - \\epsilon  \\leq 0``
"""
struct PenaltyConstraints{T,F} <: AbstractConstraints
    coef::Vector{T}
    constraints::F   # constraints(x) returns the value of the constraint-functions at x
    bounds::ConstraintBounds{T}
    PenaltyConstraints(penalty::Vector{T}, bounds::ConstraintBounds{T}, c::F) where {T,F} =
        new{T,F}(penalty, c, bounds)
end
PenaltyConstraints(penalty::T, bounds::ConstraintBounds{T}, cf=(x)->nothing) where {T} =
    PenaltyConstraints(fill(penalty, nconstraints(bounds)+nconstraints_x(bounds)), bounds, cf)
PenaltyConstraints(penalty::AbstractVector{T}, lx::AbstractVector{T}, ux::AbstractVector{T},
                   lc::AbstractVector{T}, uc::AbstractVector{T}, cf=(x)->nothing) where {T} =
    PenaltyConstraints(penalty, ConstraintBounds(lx, ux, lc, uc), cf)
PenaltyConstraints(penalty::T, lx::AbstractVector{T}, ux::AbstractVector{T},
                   lc::AbstractVector{T}=T[], uc::AbstractVector{T}=T[],
                   cf=(x)->nothing) where {T} =
    PenaltyConstraints(penalty, ConstraintBounds(lx, ux, lc, uc), cf)
PenaltyConstraints(lx::AbstractVector{T}, ux::AbstractVector{T},
                   lc::AbstractVector{T}=T[], uc::AbstractVector{T}=T[],
                   cf=(x)->nothing) where {T} =
    PenaltyConstraints(one(T), ConstraintBounds(lx, ux, lc, uc), cf)
value(c::PenaltyConstraints, x) = c.constraints(x)
penalty(c::PenaltyConstraints, x) = penalty(c.bounds, c.coef, x, value(c, x))
function penalty!(fitness::AbstractVector{T}, c::PenaltyConstraints{T,F}, population) where {T,F}
    for (i,x) in enumerate(population)
        fitness[i] += penalty(c, x)
    end
    return fitness
end

"""
This type encodes constraints as a following additional penalty for an objective function:

``p(x) = f_{worst} + \\sum^n_{i=1} |g_i(x)|``

if ``x`` is not feasible, otherwise no penalty is applied.
"""
struct WorstFitnessConstraints{T,F} <: AbstractConstraints
    bounds::ConstraintBounds{T}
    constraints::F   # constraints(x) returns the value of the constraint-functions at x
end
value(c::WorstFitnessConstraints, x) = c.constraints(x)
apply!(c::WorstFitnessConstraints, x) = clip!(c.bounds, x)
function penalty!(fitness::AbstractVector{T}, c::WorstFitnessConstraints{T,F}, population) where {T,F}
    worst = maximum(fitness)
    for (i,x) in enumerate(population)
        cv = value(c, x)
        p = zeros(size(cv))
        if !isfeasible(c.bounds, x, cv)
            for (i,j) in enumerate(c.bounds.eqc)
                p[j] = cv[j] - c.bounds.valc[i] - eps()
            end
            for (i,j) in enumerate(c.bounds.ineqc)
                p[j] = c.bounds.σc[i]*(c.bounds.bc[i]-cv[j])
            end
            fitness[i] = worst + sum(abs, p)
        end
    end
    return fitness
end

"""
This type provides an additional type constraints on the varaibles required for mixed integer optimization problmes.
"""
struct MixedTypePenaltyConstraints{C<:AbstractConstraints} <: AbstractConstraints
    penalty::C
    types::Vector{DataType}
end
value(c::MixedTypePenaltyConstraints, x) = value(c.penalty, x)
penalty!(fitness, c::MixedTypePenaltyConstraints, population) = penalty!(fitness, c.penalty, population)
function apply!(c::MixedTypePenaltyConstraints, x)
    y = apply!(c.penalty, x)
    for (i,t) in enumerate(c.types)
        if t <: Integer
            y[i] = round(t, y[i])
        end
    end
    return y
end

# Utilities

function clip!(bounds::ConstraintBounds{T}, x) where {T}
    for (i,j) in enumerate(bounds.eqx)
        x[j] = bounds.valx[i]
    end
    for (i,j) in enumerate(bounds.ineqx)
        if bounds.σx[i]*(x[j] - bounds.bx[i]) < 0
            x[j] = bounds.bx[i]
        end
    end
    x
end

function penalty(bounds::ConstraintBounds{T}, coeff::AbstractVector{T},
                 x, c::Union{AbstractVector,Nothing}=nothing) where {T}
    penalty = 0
    xc = nconstraints_x(bounds)
    for (i,j) in enumerate(bounds.eqx)
        penalty += coeff[j]*(x[j] - bounds.valx[i] - eps())^2
    end
    for (i,j) in enumerate(bounds.ineqx)
        penalty += coeff[j]*max(0, bounds.σx[i]*(bounds.bx[i]-x[j]))^2
    end
    if c !== nothing
        for (i,j) in enumerate(bounds.eqc)
            penalty += coeff[xc+j]*max(0, c[j] - bounds.valc[i] - eps())^2
        end
        for (i,j) in enumerate(bounds.ineqc)
            penalty += coeff[xc+j]*max(0, bounds.σc[i]*(bounds.bc[i]-c[j]))^2
        end
    end
    return penalty
end


text:
using Distributed
# add processes to workspace
addprocs(2)

@everywhere using ModelingToolkit, OrdinaryDiffEq

# create the Lorenz system
@everywhere @parameters t σ ρ β
@everywhere @variables x(t) y(t) z(t)
@everywhere D = Differential(t)

@everywhere eqs = [D(x) ~ σ * (y - x),
    D(y) ~ x * (ρ - z) - y,
    D(z) ~ x * y - β * z]

@everywhere @named de = ODESystem(eqs)
@everywhere ode_func = ODEFunction(de, [x, y, z], [σ, ρ, β])

@everywhere u0 = [19.0, 20.0, 50.0]
@everywhere params = [16.0, 45.92, 4]

@everywhere ode_prob = ODEProblem(ode_func, u0, (0.0, 10.0), params)

@everywhere begin
    using OrdinaryDiffEq
    using ModelingToolkit

    function solve_lorenz(ode_problem)
        print(solve(ode_problem, Tsit5()))
    end
end

solve_lorenz(ode_prob)

future = @spawn solve_lorenz(ode_prob)
@test_broken fetch(future)


text:
using Test
using LinearAlgebra
using Ferrite
using Ferrite: cellid, getcoordinates, CellIterator

using TopOpt
using TopOpt.TopOptProblems: boundingbox, nnodespercell, getgeomorder, getmetadata, getdh, getE, getdim
using TopOpt.TrussTopOptProblems: getA, default_quad_order
# if get(ENV, "CI", nothing) != "true"
#     import Makie
#     using TopOpt.TrussTopOptProblems.TrussVisualization: visualize
# end

include("utils.jl")

problem_json = ["mgz_truss1.json", "mgz_truss2.json", "mgz_truss3.json"]
u_solutions = [
    1e-3 * vcat([2.41, 0.72], zeros(2*2)),
    1e-3 * vcat([0.1783, 2.7222, -0.4863], zeros(4*3)),
    1e-3 * vcat([0.871, 1.244], [0,0], [-0.193, 0]),
]
ins_dir = joinpath(@__DIR__, "instances", "fea_examples");

@testset "Truss problem solve - $(problem_json[i])" for i in 1:length(problem_json)
    # i = 3
    file_name = problem_json[i]
    problem_file = joinpath(ins_dir, file_name)

    node_points, elements, mats, crosssecs, fixities, load_cases = parse_truss_json(problem_file);
    ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
    loads = load_cases["0"]

    problem = TrussProblem(Val{:Linear}, node_points, elements, loads, fixities, mats, crosssecs);

    @test getdim(problem) == ndim
    @test Ferrite.getncells(problem) == ncells
    @test getE(problem) == [m.E for m in mats]
    @test problem.black == problem.white == falses(ncells)
    @test problem.force == loads
    @test problem.varind == 1:ncells
    grid = problem.ch.dh.grid
    @test length(grid.cells) == ncells

    @test getgeomorder(problem) == 1
    @test nnodespercell(problem) == 2

    quad_order=default_quad_order(problem)
    elementinfo = ElementFEAInfo(problem, quad_order, Val{:Static});
    As = getA(problem)
    Es = getE(problem)
    for cell in CellIterator(getdh(problem))
        cellidx = cellid(cell)
        coords = getcoordinates(cell)
        L = norm(coords[1] - coords[2])
        A = As[cellidx]
        E = Es[cellidx]
        @test elementinfo.cellvolumes[cellidx] ≈ L * A

        Γ = global2local_transf_matrix(coords...)
        Ke_m = (A*E/L)*Γ'*[1 -1; -1 1]*Γ
        Ke = elementinfo.Kes[cellidx]
        @test Ke_m ≈ Ke
    end

    solver = FEASolver(Displacement, Direct, problem)
    solver()

    # if get(ENV, "CI", nothing) != "true"
    #     ## TODO plot analysis result with
    #     fig = visualize(problem, solver.u)
    #     Makie.display(fig)
    # end

    # we use kN for force and m for length
    # thus, pressure/modulus is in kN/m
    # the textbook uses a quite rough rounding scheme...
    # 0.3 mm error
    to_K_full = solver.globalinfo.K.data
    @assert norm(solver.u - u_solutions[i]) < 3e-4

end # end test set


text:
using Gurobi, MathProgBase, Compat.Test

using Compat: undef

@testset "Large Coefficients" begin

    @testset "Large Objective Coefficients" begin
        # Min   1.1e100x
        #       x >= 1
        m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))
        MathProgBase.loadproblem!(m, Array{Float64}(undef, 0, 1), [1],[Inf], [1.1e100], Float64[], Float64[], :Min)
        MathProgBase.optimize!(m)
        @test MathProgBase.getsolution(m) == [1.0]
        @test MathProgBase.getobjval(m) == 1e100

        # Max   -1.1e100x
        #       x >= 1
        m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))
        MathProgBase.loadproblem!(m, Array{Float64}(undef, 0, 1), [1],[Inf], [-1.1e100], Float64[], Float64[], :Max)
        MathProgBase.optimize!(m)
        @test MathProgBase.getsolution(m) == [1.0]
        @test MathProgBase.getobjval(m) == -1e100
    end

    @testset "Large Variable Bounds" begin
        # Min   x
        #       x >= 1.1e30
        m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))
        MathProgBase.loadproblem!(m, Array{Float64}(undef, 0, 1), [1.1e30],[Inf], [1], Float64[], Float64[], :Min)
        MathProgBase.optimize!(m)
        @test MathProgBase.status(m) == :Infeasible

        # Max   x
        #       x <= -1.1e30
        m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))
        MathProgBase.loadproblem!(m, Array{Float64}(undef, 0, 1), [-Inf],[-1.1e30], [1], Float64[], Float64[], :Max)
        MathProgBase.optimize!(m)
        @test MathProgBase.status(m) == :Infeasible

        # Min   x
        #       x >= -1.1e30
        m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))
        MathProgBase.loadproblem!(m, Array{Float64}(undef, 0, 1), [-1.1e30],[Inf], [1], Float64[], Float64[], :Min)
        MathProgBase.optimize!(m)
        @test MathProgBase.status(m) == :Unbounded

        # Max   x
        #       x <= 1.1e30
        m = MathProgBase.LinearQuadraticModel(GurobiSolver(OutputFlag=0))
        MathProgBase.loadproblem!(m, Array{Float64}(undef, 0, 1), [-Inf],[1.1e30], [1], Float64[], Float64[], :Max)
        MathProgBase.optimize!(m)
        @test MathProgBase.status(m) == :Unbounded

    end
end


text:
"""
    NOCB(; limit=nothing)

Next observation carried backward (NOCB) iterates backwards through the `data` and fills
missing data with the next existing observation.

See also:
- [`Impute.LOCF`](@ref): Last Observation Carried Forward

!!! Missing elements at the tail of the array may not be imputed if there is no
existing observation to carry backward. As a result, this method does not guarantee
that all missing values will be imputed.

# Keyword Arguments
* `limit::Union{UInt, Nothing}`: Optionally limits the amount of consecutive missing values to replace.

# Example
```jldoctest
julia> using Impute: NOCB, impute

julia> M = [1.0 2.0 missing missing missing 6.0; 1.1 missing missing 4.4 5.5 6.6]
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0       missing   missing   missing  6.0
 1.1   missing  missing  4.4       5.5       6.6

julia> impute(M, NOCB(); dims=:rows)
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0  6.0  6.0  6.0  6.0
 1.1  4.4  4.4  4.4  5.5  6.6

julia> impute(M,  NOCB(; limit=2); dims=:rows)
2×6 Matrix{Union{Missing, Float64}}:
 1.0  2.0   missing  6.0  6.0  6.0
 1.1  4.4  4.4       4.4  5.5  6.6
```
"""
struct NOCB <: Imputor
    limit::Union{UInt, Nothing}
end

NOCB(; limit=nothing) = NOCB(limit)

function _impute!(data::AbstractVector{Union{T, Missing}}, imp::NOCB) where T
    @assert !all(ismissing, data)
    end_idx = findlast(!ismissing, data)
    count = 1

    for i in end_idx - 1:-1:firstindex(data)
        if ismissing(data[i])
            if imp.limit === nothing
                data[i] = data[i+1]
            elseif count <= imp.limit
                data[i] = data[end_idx]
                count += 1
            end
        else
            end_idx = i
            count = 1
        end
    end

    return data
end


text:
using CSV, DataFrames
using Plots
source_data = DataFrame(CSV.File("./DataSum.csv"))

scatter(1:137, source_data.case, xlabel = "Days after March 1", ylabel = "Daily Cases", label = "Real Data", ms = 5)
plot!(1:137, source_data.DeepLearningIt, label = "Deep Learning Method", w = 2)
plot!(1:137, source_data.EpiEstimIt, label = "EpiEstim Method", w = 2)
plot!(1:137, source_data.EpiNow2It, label = "EpiNow2 Method", w = 2)
plot!(1:137, source_data.EmsembleIt, label = "Ensemble Method", w = 2)

savefig("Figures/comparisonensembleIt.png")

plot(1:137, source_data.DeepLearningRt, label = "Deep Learning Method", w = 2)
plot!(1:137, source_data.EpiEstimRt, label = "EpiEstim Method", w = 2)
plot!(1:137, source_data.EpiNow2Rt, label = "EpiNow2 Method", w = 2)
plot!(1:137, source_data.EmsembleRt, label = "Ensemble Method", w = 2)
plot!(1:137, source_data.KalmanRt, label = "Kalman Method", w = 2)
xlabel!("Days after March 1")
ylabel!("Effective Reproduction Number")

savefig("Figures/comparisonensembleRt.png")


text:
@testset "Pegasus-like instance has the correct ground state energy" begin
    ground_energy = -23.301855

    m = 3
    n = 4
    t = 3

    β = 1.

    schedule = 1.

    L = n * m * t
    states_to_keep = 20

    instance = "$(@__DIR__)/instances/pathological/cross_$(m)_$(n)_dd.txt"

    ig = ising_graph(instance)

    fg = factor_graph(
        ig,
        spectrum=full_spectrum,
        cluster_assignment_rule=super_square_lattice((m, n, t))
    )

    for transform ∈ all_lattice_transformations
        peps = PEPSNetwork{true}(m, n, fg, transform, β=β)
        update_gauges!(peps, :rand)
        sol = low_energy_spectrum(peps, states_to_keep)
        @test first(sol.energies) ≈ ground_energy
    end
end


text:
abstract type AbstractResult end

"""
    Result(total, duration, docs)

Formatted response object from redis.

# Fields
- `total::Integer`: Number of results found for the
- `duration::Float64`: Length of time taken to find all results.
- `docs::Vector{Document}`: Documents objects that are returned from  a search.
"""
mutable struct Result <: AbstractResult
    total::Integer
    duration::Float64
    docs::Vector{Document}
end

"""
    Results(res, has_content, duration; has_payload=false, with_scores=false)

Create the Result object from a Redis search response.

# Arguments
- `res`: Response from Redis
- `has_content`: Bool flag for if there is content returned
- `duration`: Time taken for query to return results

# Keywords
- `has_payload::Bool`: Flag for payload bring returned.
- `with_scores::Bool`: Flag for scores being returned.
"""
function Results(
    res,
    has_content::Bool,
    duration::Number;
    has_payload::Bool=false,
    with_scores::Bool=false,
)
    total = res[1]
    docs = Vector{Document}()

    step = 1
    if has_content
        step += 1
    end

    if has_payload
        step += 1
    end

    if with_scores
        step += 1
    end

    offset = 1
    if with_scores
        offset = 2
    end

    for i in 2:step:length(res)
        id = res[i]
        payload = Dict()
        field_offset = offset
        score = nothing

        if has_payload
            payload = res[i + offset]
            field_offset = offset + 1
        end

        if with_scores
            score = parse(Float64, res[i + 1])
        end

        fields = Dict()

        if has_content
            hset_fields = res[i + field_offset]
            fields = Dict(
                Symbol(hset_fields[j]) => hset_fields[j + 1] for
                j in 1:2:length(hset_fields)
            )
        end

        doc = Document(id; score=score, payload=payload, fields...)
        push!(docs, doc)
    end

    @info "Result: $total total"
    return Result(total, duration, docs)
end


text:
"""
    siegel_transform(τ)

Compute Siegel transformation on input Riemann matrix.
"""
function siegel_transform(τ::Array{<:Number})
    M = copy(τ);
    g = size(τ)[1];
    Γ = Matrix(I, 2*g, 2*g); # transformation
    A = [0 zeros(1, g-1); zeros(g-1, 1) Matrix(I, g-1, g-1)];
    B = [-1 zeros(1, g-1); zeros(g-1, 1) zeros(g-1, g-1)];
    C = [1 zeros(1, g-1); zeros(g-1, 1) zeros(g-1, g-1)];
    D = [0 zeros(1, g-1); zeros(g-1, 1) Matrix(I, g-1, g-1)];
    Q = [A B; C D];
    e = 0; # first entry of riemann matrix
    while e < 1
        # transform T such that the shortest vector of T is at the first column
        Y = 0.5*(imag(M) + transpose(imag(M)));
        T = (size(Y)[1] > 1 ? convert(Array{Float64, 2}, cholesky(Symmetric(Y)).U) : sqrt.(Y));
        U = hkz_transform!(T);
        M = transpose(U)*real(M)*U + im*transpose(T)*T;
        Γ = [transpose(U) zeros(g, g); zeros(g, g) inv(U)]*Γ;
        # set the real part of τ to have entries of absolute value at most 1/2
        b = -round.(real(M));
        M += b;
        Γ = [Matrix(I, g, g) b; zeros(g, g) Matrix(I, g, g)]*Γ;
        e = abs(M[1]);
        if e < 1
            M = (A*M + B)*inv(C*M + D); # Y is no longer symmetric after this operation
            Γ = Q*Γ;
        end
    end
    b = -round.(real(M));
    M += b;
    Γ = [Matrix(I, g, g) b; zeros(g, g) Matrix(I, g, g)]*Γ;
    return [Γ, M];
end

"""
    symplectic_transform(Γ, τ)

Compute symplectic group action of Γ on τ.
"""
function symplectic_transform(Γ::Array{<:Number}, τ::Array{<:Number})
    g = size(τ)[1];
    return (Γ[1:g, 1:g]*τ + Γ[1:g, g+1:2*g])*inv(Γ[g+1:2*g, 1:g]*τ + Γ[g+1:2*g, g+1:2*g]);
end

"""
    hkz_transform!(M)

Compute a unimodular matrix Z such that the first column of M*Z is the shortest vector of the lattice generated by M, where M is an upper triangular matrix. Transform the input matrix accordingly.
"""
function hkz_transform!(M::Array{<:Real})
    n = size(M)[1];
    z = svp(M)[3]; # M*z is shortest vector of M
    Z = Matrix{Float64}(I, n, n); # Z is the transform, iteratively computed such that inv(Z)*z = [1; 0; ...; 0], or first column of Z is z
    z2 = z[n];
    for j = n:-1:2
        # each iteration computes a 2x2 transform U such that inv(U) transforms [z[j-1]; z[j]] into a unit vector
        z1 = z[j-1];
        d, a, b = gcdx(z1, z2); # gcd of z[j-1], z[j]
        if d != 0
            U = [z1/d -b; z2/d a]; #  inv(U)*[z1; z2] = [d; 0]
        else
            U = Matrix(I, 2, 2);
        end
        z2 = d;
        Z[:, j-1:j] = Z[:, j-1:j]*U;
        M[:, j-1:j] = M[:, j-1:j]*U;
        # compute Givens rotation so that M remains upper triangular
        G = givens_rotation(M[j-1,j-1], M[j,j-1]);
        M[j-1:j, j-1:n] = G*M[j-1:j, j-1:n];
    end
    return Z;
end

"""
    givens_rotation(x, y)

Compute a Givens rotation, which is a 2x2 rotation matrix G such that G[x; y] has second entry zero.
"""
function givens_rotation(x::Real, y::Real)
    if abs(y) >= abs(x) # split into cases to reduce numerical error
        t = x/y;
        s = -1/sqrt(1+t^2);
        c = -t*s;
    else
        t = y/x;
        c = 1/sqrt(1+t^2);
        s = -t*c;
    end
    G = [c -s; s c];
    return G;
end


text:
_magic = "⛾ gt"
_version = 1

function minutype(n::Integer)
    @assert n ≥ 0
    if n < 2^8
        return UInt8
    elseif n < 2^16
        return UInt16
    elseif n < 2^32
        return UInt32
    elseif n < 2^64
        return UInt64
    end
    error("No type big enough")
end

function writegt(io::IO, g::AGraphOrDiGraph)
    write(io, _magic)
    write(io, UInt8(_version))
    write(io, false) # endiannes, false=little endian
    write(io, UInt64(0)) # length of comment string
    write(io, is_directed(g))
    write(io, UInt64(nv(g)))

    T = minutype(nv(g))
    writegt_adj(io, g, T)
    write(io, UInt64(0)) # num of property maps

    return 1
end

function writenetgt(io::IO, g::ANetOrDiNet)
    write(io, _magic)
    write(io, UInt8(_version))
    write(io, false) # endiannes, false=little endian
    write(io, UInt64(0)) # length of comment string
    write(io, is_directed(g))
    write(io, UInt64(nv(g)))

    T = minutype(nv(g))
    writegt_adj(io, g, T)
    writegt_props(io, g)
    return 1
end

function writegt_adj(io::IO, g::ADiGraph, ::Type{T}) where T
    for i=1:nv(g)
        write(io, UInt64(out_degree(g, i)))
        for v in out_neighbors(g, i)
            write(io, T(v-1))
        end
    end
end

function writegt_adj(io::IO, g::AGraph, ::Type{T}) where T
    for i=1:nv(g)
        neigs = out_neighbors(g, i)
        k = count(j -> j>=i, neigs)
        write(io, UInt64(k))
        for j in neigs
            if j >= i
                write(io, T(j-1))
            end
        end
    end
end

function writegt_props(io::IO, g::ANetOrDiNet)
    nprop = length(gprop(g)) + length(vprop(g)) + length(eprop(g))
    write(io, nprop) # num of property maps
    # @show nprop
    #graph props
    for (pname, p) in gprop(g)
        write(io, UInt8(0)) #property type (graph/edge/vertex)
        writegt_prop(io, pname)
        pvaln = findfirst(isequal(typeof(p)), gtpropmap)
        write(io, UInt8(pvaln-1))
        writegt_prop(io, p)
    end
    #vertex props
    for (pname, p) in vprop(g)
        write(io, UInt8(1)) #property type (graph/edge/vertex)
        writegt_prop(io, pname)
        pvaln = findfirst(isequal(valtype(p)), gtpropmap)
        write(io, UInt8(pvaln-1))

        for i=1:nv(g)
            writegt_prop(io, p[i])
        end
    end
    #edge props
    for (pname, p) in eprop(g)
        write(io, UInt8(2)) #property type (graph/edge/vertex)
        writegt_prop(io, pname)
        pvaln = findfirst(isequal(valtype(p)), gtpropmap)
        write(io, UInt8(pvaln-1))
        if is_directed(g)
            for i=1:nv(g)
                for e in out_edges(g, i)
                    writegt_prop(io, p[e])
                end
            end
        else
            for i=1:nv(g)
                for e in out_edges(g, i)
                    dst(e) < i && continue
                    writegt_prop(io, p[e])
                end
            end
        end
    end
end

function readgt(io::IO, ::Type{G}) where G<:AGraphOrDiGraph
    @assert String(read(io, 6)) == _magic "gt file not correctly formatted"
    ver = read(io, UInt8)  ## version
    indian = read(io, Bool)
    @assert indian == false
    lencomment = read(io, UInt64)
    read(io, lencomment)
    isdir = read(io, Bool)
    n = read(io, UInt64)
    T = minutype(n)
    g = isdir ? digraph(G(Int(n))) : graph(G(Int(n)))

    readgt_adj!(io, g, T)
    return g
end

function readnetgt(io::IO, ::Type{G}) where G<:ANetOrDiNet
    @assert String(read(io, 6)) == _magic "gt file not correctly formatted"
    ver = read(io, UInt8)  ## version
    indian = read(io, Bool)
    @assert indian == false
    lencomment = read(io, UInt64)
    read(io, lencomment) #skip comments
    isdir = read(io, Bool)
    n = read(io, UInt64)
    T = minutype(n)
    g = isdir ? digraph(G(Int(n))) : graph(G(Int(n)))

    readgt_adj!(io, g, T)
    readgt_props!(io, g)
    return g
end

function readgt_props!(io::IO, g::ANetOrDiNet)
    nprops = read(io, UInt64)
    for i=1:nprops
        ptype = read(io, UInt8)
        num = ptype == 0 ? 1 :
              ptype == 1 ? nv(g) : ne(g)
        namelen = read(io, UInt64)
        pname = String(read(io, namelen))
        pvaln = read(io, UInt8) + 1
        # @show num namelen pname pvaln
        pvaln ∈ [6,13,14] && error("not implemented") #TODO
        T = gtpropmap[pvaln]
        m = readgt_prop(io, T, num)

        if ptype == 0
            set_graph_property!(g, pname, m[1])
        elseif ptype == 1
            add_vertex_property!(g, pname, m)
        elseif ptype == 2
            add_edge_property!(g, pname, EdgeMap(g, m))
        end
    end

end

const Num = Union{Bool,Int16,Int32,Int64,Float64}

const gtpropmap = DataType[   Bool,                   #0x00
                              Int16,                  #0x01
                              Int32,                  #0x02
                              Int64,                  #0x03
                              Float64,                #0x04
                              Float64,               #0x05 #TODO => should be float128
                              String,                 #0x06
                              Vector{Bool},           #0x07
                              Vector{Int16},          #0x08
                              Vector{Int32},          #0x09
                              Vector{Int64},          #0x0a
                              Vector{Float64},        #0x0b
                              Vector{Float64},       #0x0c #TODO => should be float128
                              Vector{String}          #0x0d
                        ]

readgt_prop(io, ::Type{T}, num) where {T<:Num} = read!(io, Vector{T}(undef, num))
function readgt_prop(io, ::Type{Vector{T}}, num) where {T <: Num}
    [(l = read(io, UInt64); read!(io, Vector{T}(undef, l))) for _=1:num]
end
readgt_prop(io, ::Type{String}, num) = [(l = read(io, UInt64); String(read(io, l))) for _=1:num]

writegt_prop(io, x::T) where {T<:Num} = write(io, x)
writegt_prop(io, x::Vector{T}) where {T<:Num} =
    (write(io, length(x)); write(io, x))
writegt_prop(io, x::String) = (write(io, sizeof(x)); write(io, x))

function readgt_adj!(io::IO, g::ADiGraph, ::Type{T}) where T
    for i=1:nv(g)
        k = read(io, UInt64)
        for _=1:k
            j = read(io, T) + 1
            # unsafe_add_edge!(g, i, j)
            add_edge!(g, i, j)
        end
    end
    # rebuild!(g)
end

function readgt_adj!(io::IO, g::AGraph, ::Type{T}) where T
    for i=1:nv(g)
        k = read(io, UInt64)
        for _=1:k
            j = read(io, T) + 1
            # @assert j >= i
            # unsafe_add_edge!(g, i, j)
            add_edge!(g, i, j)
        end
    end
    # rebuild!(g)
end

filemap[:gt] = (readgt, writegt, readnetgt, writenetgt)


text:
"""Linear interpolation of a variable u in the x-direction.
m,n = size(ux) must be m+1,n = size(u)."""
function Ix!(ux::AbstractMatrix,u::AbstractMatrix)
    m, n = size(ux)
    @boundscheck (m+1,n) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            ux[i,j] = one_half*(u[i+1,j] + u[i,j])
        end
    end
end

""" Linear interpolation a variable u in the y-direction.
    m,n = size(uy) must be m,n+1 = size(u)."""
function Iy!(uy::AbstractMatrix,u::AbstractMatrix)
    m,n = size(uy)
    @boundscheck (m,n+1) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            uy[i,j] = one_half*(u[i,j+1] + u[i,j])
        end
    end
end

""" Bilinear interpolation a variable u in x and y-direction.
m,n = size(uxy) must be m+1,n+1 = size(u). """
function Ixy!(uxy::AbstractMatrix,u::AbstractMatrix)
    m,n = size(uxy)
    @boundscheck (m+1,n+1) == size(u) || throw(BoundsError())

    @inbounds for j ∈ 1:n
        for i ∈ 1:m
            uxy[i,j] = one_quart*(u[i,j] + u[i+1,j] + u[i,j+1] + u[i+1,j+1])
        end
    end
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
call = [0 0 0.0001 0.0011 0.0055 ;
        0 0 0.005 0.0412 0.1399 ;
        0 0.0403 0.4195 1.1159 1.9965 ;
        0 4.7469 7.1559 9.1754 10.9895 ;
        22.3035 23.4719 24.7110 25.9912 27.2682 ;
        39.6325 40.5313 41.4170 42.2922 43.1598 ;
        53.0966 53.7949 54.4828 55.1604 55.8281]
put = [113.2040 110.0298 106.903 103.8236 100.7935 ;
       65.6521 63.1859 60.7613 58.4042 56.1451 ;
       28.7059 26.83 25.3216 24.1584 23.2072 ;
       0 3.2581 4.2004 4.7752 5.1660 ;
       0 0.0117 0.1112 0.2689 0.44 ;
       0 0 0.0003 0.0033 0.0118 ;
       0 0 0 0 0.0001]

put[1,1]

T, S, r, N, dx = 1, 100, 0.06, 4, 0.2524

dt = T/N
edx = exp(dx)

## Initialise asset prices at maturity

St = zeros(Float64, 2N+1)
Q = zeros(Float64, (2N+1, N+1))

  ## creating a shifted ranged for addressing the Array in the usual i,j notation
J = sub(1:(2N+1), (1+(N+1):((2N+1)+(N+1))))

## Asset price at node N,-N
St[2N+1] = S * exp(-N * dx)

## Asset prices from bottom -N, to N, at time step N
for j = 2N:-1:1
St[j] = St[j+1] * edx
end

for i = 1:N

    for j = i:-1:0

            sum = 0

            for k = i:-1:(j+1)

                 sum = sum + Q[(N+1)-k, i]*(St[(N+1)-k] - St[(N+1)-j])

            end

        Q[(N+1)-j, i] = (call[(N+1)-j, i] - sum) / (St[(N+1)-j] - St[(N+1)-j+1])

    end

    for j = -i:1:-1

            sum = 0
            println("I am here $j")

            for k = -i:-1:j
                     println("k is $k")
                sum = sum + Q[(N+1)-k, i]*(St[(N+1)+j-1] - St[(N+1)-k])
                #println("$sum")
            end

         Q[N+1-j, i] = (put[N-j-1, i] - sum) / (St[N+1+j-1] - St[N+1+j])

    end

end

Q

Q = zeros(Float64, (2N+1, N+1))

for j in -4:1:-1
    x = put[N-j-1,1]
    println("x is $x")
end
call[]

  ## initialise Array for the tree

  C = zeros(Float64, (2N+1, N+1))

  ## Option Values at expiry, time step = N, in rows -N to
  for j = -N:N
    C[J[j], N+1] = max(0, St[J[j]] - K)
  end

  ## Stepping back through the lattice

  for i in range(N-1, -1, N)
    for j = -i:i
      C[J[j], i+1] = disc * (pu * C[J[j]+1, i+2] + pm * C[J[j], i+2] + pd * C[J[j]-1, i+2])
    end
  end

C


text:
@testset "trunc" begin

@test trunc(Int64, 3.14) == 3
@test trunc(Int64, 3.41) == 3
@test trunc(Int64, 3.51) == 3
@test trunc(Int64, 3.91) == 3

end

@testset "round" begin
@test round(1.4) == 1
@test round(1.5) == 2
end

@testset "parse" begin
@test parse(Int64, "13") == 13
# Why would this work???
@test parse(Float64, "13.1415926") == 13.1415926
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Operate <: Edit
        operation::String
        from::Array{Symbol,1}
        to::Array{Symbol,1}
        input::Array{Symbol,1}
        output::Symbol
    end

Perform an arithmetic operation across multiple DataFrame columns.

# Arguments
- `operation::String`: operation to perform (+, -, *, /)
- `from::Array{Symbol,1}`: name of original comment column (ex. units)
- `to::Array{Symbol,1}`: name of new comment column (ex. units)
- `input::Array{Symbol,1}`: names of columns on which to operate
- `output::Symbol`: name of result column
"""
mutable struct Operate <: Edit
    "operation to perform (+, -, *, /)"
    operation::String
    "name of original comment column (ex. units)"
    from::Array{Symbol,1}
    "name of new comment column (ex. units)"
    to::Array{Symbol,1}
    "names of columns on which to operate"
    input::Array{Symbol,1}
    "name of result column"
    output::Symbol
end

function Operate(; operation, from, to, input, output, )
    Operate(operation, from, to, input, output, )
end

"""Get [`Operate`](@ref) `operation`."""
get_operation(value::Operate) = value.operation
"""Get [`Operate`](@ref) `from`."""
get_from(value::Operate) = value.from
"""Get [`Operate`](@ref) `to`."""
get_to(value::Operate) = value.to
"""Get [`Operate`](@ref) `input`."""
get_input(value::Operate) = value.input
"""Get [`Operate`](@ref) `output`."""
get_output(value::Operate) = value.output

"""Set [`Operate`](@ref) `operation`."""
set_operation!(value::Operate, val) = value.operation = val
"""Set [`Operate`](@ref) `from`."""
set_from!(value::Operate, val) = value.from = val
"""Set [`Operate`](@ref) `to`."""
set_to!(value::Operate, val) = value.to = val
"""Set [`Operate`](@ref) `input`."""
set_input!(value::Operate, val) = value.input = val
"""Set [`Operate`](@ref) `output`."""
set_output!(value::Operate, val) = value.output = val

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
abstract AbstractOp{T,N}

##################################
#  Matrix Operators
##################################

type MatrixOp{N}<:AbstractOp{AbstractMatrix,N}
    name::AbstractString
    mat::AbstractMatrix
end

MatrixOp(num::Integer,name::AbstractString,mat::AbstractMatrix) = MatrixOp{num}(name,mat)

function call{N}(matop::MatrixOp{N})
    return [matop.mat[:,i] for i=1:size(matop.mat,2)]
end

function call{N}(matop::MatrixOp{N},state::AbstractVector)
    return matop.mat*state
end

function show{N}(io::IO,matop::MatrixOp{N})
    println("$N bits matrix operator $(matop.name):")
    show(matop.mat)
end

OP_Hadamard = MatrixOp(1,"Hadamard",hadamard)
# Pauli Groups
OP_sigmax   = MatrixOp(1,"Pauli Sigma X",σ₁)
OP_sigmay   = MatrixOp(1,"Pauli Sigma Y",σ₂)
OP_sigmaz   = MatrixOp(1,"Pauli Sigma Z",σ₃)

# TODO
# this comes from linalg/uniformscaling.jl
# same operators should be overloaded
immutable IdentityOp{T<:Number}<:AbstractOp{T}
    λ::T
end

OP_I = IdentityOp(1)

##################################
# Function Operators
##################################

const FUNCTION_OP_PARA_INF = -1

"""
Function operators should be able to accept AbstractSparseVector inputs

The member `f` should be an one input only function
"""
type FunctionOp{N}<:AbstractOp{Function,N}
    name::AbstractString
    f::Function
end

function basis(n::Int)
    return [i->sparsevec(Dict[i=>1]) for i=1:2^n]
end

function call{N}(funcop::FunctionOp{N})
    res = Array(Array{Complex,1},0)
    for i in basis(N)
        push!(res,funcop.f(i))
    end
    return res
end

function call{N}(funcop::FunctionOp{N},state::AbstractVector)
    return funcop(state)
end

function TimeOp{N}(state::ComputState{N};Hamiltonian=I,dt=1e-6)
    return expmv(-im*dt,Hamiltonian,coeffs(state))
end


text:
using HomotopyContinuation
using Base.Test
import DynamicPolynomials
const PolyImpl = DynamicPolynomials

include("base_utilities_test.jl")

include("straight_line_test.jl")
include("gamma_trick_test.jl")
include("fixedhomotopy_test.jl")
include("spherical_test.jl")
include("affine_test.jl")
include("test_systems_test.jl")
include("cauchyendgame_test.jl")


text:
using BinaryBuilder

name = "Musl"
version = v"1.1.24"

# sources to build, such as mingw32, our patches, etc....
sources = [
    ArchiveSource("https://www.musl-libc.org/releases/musl-$(version).tar.gz",
                  "1370c9a812b2cf2a7d92802510cca0058cc37e66a7bedd70051f0a34015022a3"),
]

# Bash recipe for building across all platforms
script = raw"""
mkdir ${WORKSPACE}/srcdir/musl_build
cd ${WORKSPACE}/srcdir/musl_build
musl_arch()
{
    case "${target}" in
        i686*)
            echo i386 ;;
        arm*)
            echo armhf ;;
        *)
            echo ${target%%-*} ;;
    esac
}

export LDFLAGS="${LDFLAGS} -Wl,-soname,libc.musl-$(musl_target).so.1"
${WORKSPACE}/srcdir/musl-*/configure --prefix=/usr \
    --build=${MACHTYPE} \
    --host=${target} \
    --disable-multilib \
    --disable-werror \
    --enable-optimize \
    --enable-debug

make -j${nproc}

# Install to fake directory
make install DESTDIR=${WORKSPACE}/srcdir/musl_deploy

# In reality, all we want is what we need to run Musl programs, so
# we just copy the dynamic loader directly over.
mkdir -p ${prefix}/lib
loader=$(echo ${WORKSPACE}/srcdir/musl_deploy/lib/*.so*)
mv ${WORKSPACE}/srcdir/musl_deploy/usr/lib/libc.so ${prefix}/lib/$(basename "${loader}")
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [p for p in supported_platforms() if libc(p) == :musl]

# The products that we will ensure are always built
products = [
    LibraryProduct(["ld-musl-x86_64", "ld-musl-i386", "ld-musl-armhf", "ld-musl-aarch64"], :libc; dont_dlopen=true),
]

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using ..ITensorAutoHOOT
using ..ITensorAutoHOOT: SubNetwork, neighboring_tensors

function tree(sub_peps_bra::Vector, sub_peps_ket::Vector, mpo::Vector)
  @assert length(sub_peps_bra) == length(sub_peps_ket)
  @assert length(sub_peps_bra) == length(mpo)
  center = div(length(sub_peps_bra), 2)
  len = length(sub_peps_bra)
  front_tree = SubNetwork(sub_peps_bra[1], sub_peps_ket[1], mpo[1])
  for i in 2:center
    front_tree = SubNetwork(sub_peps_bra[i], sub_peps_ket[i], mpo[i], front_tree)
  end
  back_tree = SubNetwork(sub_peps_bra[len], sub_peps_ket[len], mpo[len])
  for i in (len - 1):-1:(center + 1)
    back_tree = SubNetwork(sub_peps_bra[i], sub_peps_ket[i], mpo[i], back_tree)
  end
  return SubNetwork(front_tree, back_tree)
end

function tree(sub_peps_bra::Vector, sub_peps_ket::Vector)
  @assert length(sub_peps_bra) == length(sub_peps_ket)
  center = div(length(sub_peps_bra), 2)
  len = length(sub_peps_bra)
  front_tree = SubNetwork(sub_peps_bra[1], sub_peps_ket[1])
  for i in 2:center
    front_tree = SubNetwork(sub_peps_bra[i], sub_peps_ket[i], front_tree)
  end
  back_tree = SubNetwork(sub_peps_bra[len], sub_peps_ket[len])
  for i in (len - 1):-1:(center + 1)
    back_tree = SubNetwork(sub_peps_bra[i], sub_peps_ket[i], back_tree)
  end
  return SubNetwork(front_tree, back_tree)
end

"""Returns a tree structure for a line of tensors with projectors
Parameters
----------
line_size: size of the line structure
center_index: the center index of the line
site_tensors: a function, site_tensors(i) returns a list of tensors at position i
projectors: the projectors of the line structure
Returns
-------
Two trees, one in front of the center_index and another one after center_index
Example
-------
   |  |   |
p1-p2 |  p3-p4
 | |  |   | |
 | |  |   | |
s1-s2-s3-s4-s5
here line_size=5, center_index=3, si represents the list of tensors returned by site_tensors(i),
projectors are [p1, p2, p3, p4].
Returns two trees: [[s1, p1], s2, p2] and [[s5, p4], s4, p3]
"""
function tree_w_projectors(
  line_size, center_index, site_tensors, projectors::Vector{ITensor}
)
  front_tree, back_tree = nothing, nothing
  for i in 1:(center_index - 1)
    connect_projectors = neighboring_tensors(SubNetwork(site_tensors(i)), projectors)
    if front_tree == nothing
      inputs = vcat(site_tensors(i), connect_projectors)
    else
      inputs = vcat(site_tensors(i), connect_projectors, [front_tree])
    end
    front_tree = SubNetwork(inputs)
  end
  for i in line_size:-1:(center_index + 1)
    connect_projectors = neighboring_tensors(SubNetwork(site_tensors(i)), projectors)
    if back_tree == nothing
      inputs = vcat(site_tensors(i), connect_projectors)
    else
      inputs = vcat(site_tensors(i), connect_projectors, [back_tree])
    end
    back_tree = SubNetwork(inputs)
  end
  return front_tree, back_tree
end

function tree_w_projectors(
  sub_peps_bra::Vector, sub_peps_ket::Vector, projectors::Vector{ITensor}
)
  out_inds = inds(SubNetwork(vcat(sub_peps_bra, sub_peps_ket, projectors)))
  is_neighbor(t) = length(intersect(out_inds, inds(t))) > 0
  center_index = [i for (i, t) in enumerate(sub_peps_bra) if is_neighbor(t)]
  @assert length(center_index) == 1
  center_index = center_index[1]
  @assert is_neighbor(sub_peps_ket[center_index])
  site_tensors(i) = [sub_peps_bra[i], sub_peps_ket[i]]
  front_tree, back_tree = tree_w_projectors(
    length(sub_peps_bra), center_index, site_tensors, projectors::Vector{ITensor}
  )
  return SubNetwork([
    sub_peps_bra[center_index], sub_peps_ket[center_index], front_tree, back_tree
  ])
end

function tree_w_projectors(
  sub_peps_bra::Vector, sub_peps_ket::Vector, mpo::Vector, projectors::Vector{ITensor}
)
  out_inds = inds(SubNetwork(vcat(sub_peps_bra, sub_peps_ket, mpo, projectors)))
  is_neighbor(t) = length(intersect(out_inds, inds(t))) > 0
  center_index = [i for (i, t) in enumerate(sub_peps_bra) if is_neighbor(t)]
  @assert length(center_index) == 1
  center_index = center_index[1]
  @assert is_neighbor(sub_peps_ket[center_index])
  site_tensors(i) = [sub_peps_bra[i], sub_peps_ket[i], mpo[i]]
  front_tree, back_tree = tree_w_projectors(
    length(sub_peps_bra), center_index, site_tensors, projectors::Vector{ITensor}
  )
  return SubNetwork([
    sub_peps_bra[center_index],
    sub_peps_ket[center_index],
    mpo[center_index],
    front_tree,
    back_tree,
  ])
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    AbstractInterpreter

An abstract base class that allows multiple dispatch to determine the method of
executing Julia code.  The native Julia LLVM pipeline is enabled by using the
`NativeInterpreter` concrete instantiation of this abstract class, others can be
swapped in as long as they follow the `AbstractInterpreter` API.

If `interp` is an `AbstractInterpreter`, it is expected to provide at least the following methods:
- `InferenceParams(interp)` - return an `InferenceParams` instance
- `OptimizationParams(interp)` - return an `OptimizationParams` instance
- `get_world_counter(interp)` - return the world age for this interpreter
- `get_inference_cache(interp)` - return the runtime inference cache
- `code_cache(interp)` - return the global inference cache
"""
abstract type AbstractInterpreter end

struct ArgInfo
    fargs::Union{Nothing,Vector{Any}}
    argtypes::Vector{Any}
end

struct TriState; state::UInt8; end
const ALWAYS_FALSE     = TriState(0x00)
const ALWAYS_TRUE      = TriState(0x01)
const TRISTATE_UNKNOWN = TriState(0x02)

function tristate_merge(old::TriState, new::TriState)
    (old === ALWAYS_FALSE || new === ALWAYS_FALSE) && return ALWAYS_FALSE
    old === TRISTATE_UNKNOWN && return old
    return new
end

"""
    effects::Effects

Represents computational effects of a method call.

The effects are composed of the following set of different properties:
- `effects.consistent::TriState`: this method is guaranteed to return or terminate consistently
- `effect_free::TriState`: this method is free from externally semantically visible side effects
- `nothrow::TriState`: this method is guaranteed to not throw an exception
- `terminates::TriState`: this method is guaranteed to terminate
- `nonoverlayed::Bool`: indicates that any methods that may be called within this method
  are not defined in an [overlayed method table](@ref OverlayMethodTable)
- `notaskstate::TriState`: this method does not access any state bound to the current
  task and may thus be moved to a different task without changing observable
  behavior. Note that this currently implies that `noyield` as well, since
  yielding modifies the state of the current task, though this may be split
  in the future.
See [`Base.@assume_effects`](@ref) for more detailed explanation on the definitions of these properties.

Along the abstract interpretation, `Effects` at each statement are analyzed locally and
they are merged into the single global `Effects` that represents the entire effects of
the analyzed method (see `tristate_merge!`).
Each effect property is represented as tri-state and managed separately.
The tri-state consists of `ALWAYS_TRUE`, `TRISTATE_UNKNOWN` and `ALWAYS_FALSE`.
An effect property is initialized with `ALWAYS_TRUE` and then transitioned towards
`TRISTATE_UNKNOWN` or `ALWAYS_FALSE`. When we find a statement that has some effect,
`ALWAYS_TRUE` is propagated if that effect is known to _always_ happen, otherwise
`TRISTATE_UNKNOWN` is propagated. If a property is known to be `ALWAYS_FALSE`,
there is no need to do additional analysis as it can not be refined anyway.
Note that however, within the current data-flow analysis design, it is hard to derive a global
conclusion from a local analysis on each statement, and as a result, the effect analysis
usually propagates `TRISTATE_UNKNOWN` currently.
"""
struct Effects
    consistent::TriState
    effect_free::TriState
    nothrow::TriState
    terminates::TriState
    nonoverlayed::Bool
    notaskstate::TriState
    # This effect is currently only tracked in inference and modified
    # :consistent before caching. We may want to track it in the future.
    inbounds_taints_consistency::Bool
end
function Effects(
    consistent::TriState,
    effect_free::TriState,
    nothrow::TriState,
    terminates::TriState,
    nonoverlayed::Bool,
    notaskstate::TriState)
    return Effects(
        consistent,
        effect_free,
        nothrow,
        terminates,
        nonoverlayed,
        notaskstate,
        false)
end

const EFFECTS_TOTAL    = Effects(ALWAYS_TRUE,      ALWAYS_TRUE,      ALWAYS_TRUE,      ALWAYS_TRUE,      true, ALWAYS_TRUE)
const EFFECTS_THROWS   = Effects(ALWAYS_TRUE,      ALWAYS_TRUE,      TRISTATE_UNKNOWN, ALWAYS_TRUE,      true, ALWAYS_TRUE)
const EFFECTS_UNKNOWN  = Effects(TRISTATE_UNKNOWN, TRISTATE_UNKNOWN, TRISTATE_UNKNOWN, TRISTATE_UNKNOWN, true, TRISTATE_UNKNOWN)  # mostly unknown, but it's not overlayed at least (e.g. it's not a call)
const EFFECTS_UNKNOWN′ = Effects(TRISTATE_UNKNOWN, TRISTATE_UNKNOWN, TRISTATE_UNKNOWN, TRISTATE_UNKNOWN, false, TRISTATE_UNKNOWN) # unknown, really

function Effects(e::Effects = EFFECTS_UNKNOWN′;
    consistent::TriState = e.consistent,
    effect_free::TriState = e.effect_free,
    nothrow::TriState = e.nothrow,
    terminates::TriState = e.terminates,
    nonoverlayed::Bool = e.nonoverlayed,
    notaskstate::TriState = e.notaskstate,
    inbounds_taints_consistency::Bool = e.inbounds_taints_consistency)
    return Effects(
        consistent,
        effect_free,
        nothrow,
        terminates,
        nonoverlayed,
        notaskstate,
        inbounds_taints_consistency)
end

is_consistent(effects::Effects)   = effects.consistent === ALWAYS_TRUE
is_effect_free(effects::Effects)  = effects.effect_free === ALWAYS_TRUE
is_nothrow(effects::Effects)      = effects.nothrow === ALWAYS_TRUE
is_terminates(effects::Effects)   = effects.terminates === ALWAYS_TRUE
is_notaskstate(effects::Effects)  = effects.notaskstate === ALWAYS_TRUE
is_nonoverlayed(effects::Effects) = effects.nonoverlayed

# implies :notaskstate, but not explicitly checked here
is_foldable(effects::Effects) =
    is_consistent(effects) &&
    is_effect_free(effects) &&
    is_terminates(effects)

is_total(effects::Effects) =
    is_foldable(effects) &&
    is_nothrow(effects)

is_removable_if_unused(effects::Effects) =
    is_effect_free(effects) &&
    is_terminates(effects) &&
    is_nothrow(effects)

function encode_effects(e::Effects)
    return (e.consistent.state << 0) |
           (e.effect_free.state << 2) |
           (e.nothrow.state << 4) |
           (e.terminates.state << 6) |
           (UInt32(e.nonoverlayed) << 8) |
           (UInt32(e.notaskstate.state) << 9)
end
function decode_effects(e::UInt32)
    return Effects(
        TriState((e >> 0) & 0x03),
        TriState((e >> 2) & 0x03),
        TriState((e >> 4) & 0x03),
        TriState((e >> 6) & 0x03),
        _Bool(   (e >> 8) & 0x01),
        TriState((e >> 9) & 0x03),
        false)
end

function tristate_merge(old::Effects, new::Effects)
    return Effects(
        tristate_merge(
            old.consistent, new.consistent),
        tristate_merge(
            old.effect_free, new.effect_free),
        tristate_merge(
            old.nothrow, new.nothrow),
        tristate_merge(
            old.terminates, new.terminates),
        old.nonoverlayed & new.nonoverlayed,
        tristate_merge(
            old.notaskstate, new.notaskstate),
        old.inbounds_taints_consistency | new.inbounds_taints_consistency)
end

struct EffectsOverride
    consistent::Bool
    effect_free::Bool
    nothrow::Bool
    terminates_globally::Bool
    terminates_locally::Bool
    notaskstate::Bool
end

function encode_effects_override(eo::EffectsOverride)
    e = 0x00
    eo.consistent && (e |= 0x01)
    eo.effect_free && (e |= 0x02)
    eo.nothrow && (e |= 0x04)
    eo.terminates_globally && (e |= 0x08)
    eo.terminates_locally && (e |= 0x10)
    eo.notaskstate && (e |= 0x20)
    return e
end

function decode_effects_override(e::UInt8)
    return EffectsOverride(
        (e & 0x01) != 0x00,
        (e & 0x02) != 0x00,
        (e & 0x04) != 0x00,
        (e & 0x08) != 0x00,
        (e & 0x10) != 0x00,
        (e & 0x20) != 0x00)
end

"""
    InferenceResult

A type that represents the result of running type inference on a chunk of code.
"""
mutable struct InferenceResult
    linfo::MethodInstance
    argtypes::Vector{Any}
    overridden_by_const::BitVector
    result                   # ::Type, or InferenceState if WIP
    src                      # ::Union{CodeInfo, OptimizationState} if inferred copy is available, nothing otherwise
    valid_worlds::WorldRange # if inference and optimization is finished
    ipo_effects::Effects     # if inference is finished
    effects::Effects         # if optimization is finished
    argescapes               # ::ArgEscapeCache if optimized, nothing otherwise
    function InferenceResult(linfo::MethodInstance,
                             arginfo#=::Union{Nothing,Tuple{ArgInfo,InferenceState}}=# = nothing)
        argtypes, overridden_by_const = matching_cache_argtypes(linfo, arginfo)
        return new(linfo, argtypes, overridden_by_const, Any, nothing,
            WorldRange(), Effects(), Effects(), nothing)
    end
end

"""
    OptimizationParams

Parameters that control optimizer operation.
"""
struct OptimizationParams
    inlining::Bool              # whether inlining is enabled
    inline_cost_threshold::Int  # number of CPU cycles beyond which it's not worth inlining
    inline_nonleaf_penalty::Int # penalty for dynamic dispatch
    inline_tupleret_bonus::Int  # extra inlining willingness for non-concrete tuple return types (in hopes of splitting it up)
    inline_error_path_cost::Int # cost of (un-optimized) calls in blocks that throw

    trust_inference::Bool

    # Duplicating for now because optimizer inlining requires it.
    # Keno assures me this will be removed in the near future
    MAX_METHODS::Int
    MAX_TUPLE_SPLAT::Int
    MAX_UNION_SPLITTING::Int

    function OptimizationParams(;
            inlining::Bool = inlining_enabled(),
            inline_cost_threshold::Int = 100,
            inline_nonleaf_penalty::Int = 1000,
            inline_tupleret_bonus::Int = 250,
            inline_error_path_cost::Int = 20,
            max_methods::Int = 3,
            tuple_splat::Int = 32,
            union_splitting::Int = 4,
            trust_inference::Bool = false
        )
        return new(
            inlining,
            inline_cost_threshold,
            inline_nonleaf_penalty,
            inline_tupleret_bonus,
            inline_error_path_cost,
            trust_inference,
            max_methods,
            tuple_splat,
            union_splitting
        )
    end
end

"""
    InferenceParams

Parameters that control type inference operation.
"""
struct InferenceParams
    ipo_constant_propagation::Bool
    aggressive_constant_propagation::Bool
    unoptimize_throw_blocks::Bool

    # don't consider more than N methods. this trades off between
    # compiler performance and generated code performance.
    # typically, considering many methods means spending lots of time
    # obtaining poor type information.
    # It is important for N to be >= the number of methods in the error()
    # function, so we can still know that error() is always Bottom.
    MAX_METHODS::Int
    # the maximum number of union-tuples to swap / expand
    # before computing the set of matching methods
    MAX_UNION_SPLITTING::Int
    # the maximum number of union-tuples to swap / expand
    # when inferring a call to _apply_iterate
    MAX_APPLY_UNION_ENUM::Int

    # parameters limiting large (tuple) types
    TUPLE_COMPLEXITY_LIMIT_DEPTH::Int

    # when attempting to inline _apply_iterate, abort the optimization if the
    # tuple contains more than this many elements
    MAX_TUPLE_SPLAT::Int

    function InferenceParams(;
            ipo_constant_propagation::Bool = true,
            aggressive_constant_propagation::Bool = false,
            unoptimize_throw_blocks::Bool = true,
            max_methods::Int = 3,
            union_splitting::Int = 4,
            apply_union_enum::Int = 8,
            tupletype_depth::Int = 3,
            tuple_splat::Int = 32,
        )
        return new(
            ipo_constant_propagation,
            aggressive_constant_propagation,
            unoptimize_throw_blocks,
            max_methods,
            union_splitting,
            apply_union_enum,
            tupletype_depth,
            tuple_splat,
        )
    end
end

"""
    NativeInterpreter

This represents Julia's native type inference algorithm and codegen backend.
It contains many parameters used by the compilation pipeline.
"""
struct NativeInterpreter <: AbstractInterpreter
    # Cache of inference results for this particular interpreter
    cache::Vector{InferenceResult}
    # The world age we're working inside of
    world::UInt

    # Parameters for inference and optimization
    inf_params::InferenceParams
    opt_params::OptimizationParams

    function NativeInterpreter(world::UInt = get_world_counter();
                               inf_params = InferenceParams(),
                               opt_params = OptimizationParams(),
                               )
        # Sometimes the caller is lazy and passes typemax(UInt).
        # we cap it to the current world age
        if world == typemax(UInt)
            world = get_world_counter()
        end

        # If they didn't pass typemax(UInt) but passed something more subtly
        # incorrect, fail out loudly.
        @assert world <= get_world_counter()

        return new(
            # Initially empty cache
            Vector{InferenceResult}(),

            # world age counter
            world,

            # parameters for inference and optimization
            inf_params,
            opt_params,
        )
    end
end

# Quickly and easily satisfy the AbstractInterpreter API contract
InferenceParams(ni::NativeInterpreter) = ni.inf_params
OptimizationParams(ni::NativeInterpreter) = ni.opt_params
get_world_counter(ni::NativeInterpreter) = ni.world
get_inference_cache(ni::NativeInterpreter) = ni.cache
code_cache(ni::NativeInterpreter) = WorldView(GLOBAL_CI_CACHE, get_world_counter(ni))

"""
    lock_mi_inference(ni::NativeInterpreter, mi::MethodInstance)

Hint that `mi` is in inference to help accelerate bootstrapping.
This helps us limit the amount of wasted work we might do when inference is working on initially inferring itself
by letting us detect when inference is already in progress and not running a second copy on it.
This creates a data-race, but the entry point into this code from C (`jl_type_infer`) already includes detection and restriction on recursion,
so it is hopefully mostly a benign problem (since it should really only happen during the first phase of bootstrapping that we encounter this flag).
"""
lock_mi_inference(::NativeInterpreter, mi::MethodInstance) = (mi.inInference = true; nothing)
lock_mi_inference(::AbstractInterpreter, ::MethodInstance) = return

"""
See `lock_mi_inference`.
"""
unlock_mi_inference(::NativeInterpreter, mi::MethodInstance) = (mi.inInference = false; nothing)
unlock_mi_inference(::AbstractInterpreter, ::MethodInstance) = return

"""
Emit an analysis remark during inference for the current line (`sv.pc`).
These annotations are ignored by the native interpreter, but can be used by external tooling
to annotate inference results.
"""
add_remark!(::AbstractInterpreter, sv#=::InferenceState=#, s) = return

may_optimize(::AbstractInterpreter) = true
may_compress(::AbstractInterpreter) = true
may_discard_trees(::AbstractInterpreter) = true
verbose_stmt_info(::AbstractInterpreter) = false

"""
    method_table(interp::AbstractInterpreter) -> MethodTableView

Returns a method table this `interp` uses for method lookup.
External `AbstractInterpreter` can optionally return `OverlayMethodTable` here
to incorporate customized dispatches for the overridden methods.
"""
method_table(interp::AbstractInterpreter) = InternalMethodTable(get_world_counter(interp))

"""
By default `AbstractInterpreter` implements the following inference bail out logic:
- `bail_out_toplevel_call(::AbstractInterpreter, sig, ::InferenceState)`: bail out from inter-procedural inference when inferring top-level and non-concrete call site `callsig`
- `bail_out_call(::AbstractInterpreter, rt, ::InferenceState)`: bail out from inter-procedural inference when return type `rt` grows up to `Any`
- `bail_out_apply(::AbstractInterpreter, rt, ::InferenceState)`: bail out from `_apply_iterate` inference when return type `rt` grows up to `Any`

It also bails out from local statement/frame inference when any lattice element gets down to `Bottom`,
but `AbstractInterpreter` doesn't provide a specific interface for configuring it.
"""
bail_out_toplevel_call(::AbstractInterpreter, @nospecialize(callsig), sv#=::InferenceState=#) =
    return sv.restrict_abstract_call_sites && !isdispatchtuple(callsig)
bail_out_call(::AbstractInterpreter, @nospecialize(rt), sv#=::InferenceState=#) =
    return rt === Any
bail_out_apply(::AbstractInterpreter, @nospecialize(rt), sv#=::InferenceState=#) =
    return rt === Any

"""
    infer_compilation_signature(::AbstractInterpreter)::Bool

For some call sites (for example calls to varargs methods), the signature to be compiled
and executed at run time can differ from the argument types known at the call site.
This flag controls whether we should always infer the compilation signature in addition
to the call site signature.
"""
infer_compilation_signature(::AbstractInterpreter) = false
infer_compilation_signature(::NativeInterpreter) = true

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia