text:
"""
    DAMMfit(x::VecOrMat{<: Real}, Rₛ::Vector{Float64}, poro_val::Float64)

Fit the DAMM model parameters to data.

# Examples
```julia-repl
julia> df = DAMMfdata(100) # generates a fake dataset
100×3 DataFrame
 Row │ Tₛ       θ        Rₛ
     │ Float64  Float64  Float64
─────┼─────────────────────────────
   1 │    27.1      0.3   4.345
   2 │    38.7      0.6  12.0106
  ⋮  │    ⋮        ⋮         ⋮
  99 │    18.6      0.5   0.894257
 100 │    19.4      0.4   3.79532
julia> p = DAMMfit(hcat(df.Tₛ, df.θ), df.Rₛ, 0.7)
(2.034002955272664e10, 71.65411256289629, 9.903541279858033e-8, 0.003688664956456453, 0.7, 0.02, 1.0)
julia> DAMM(hcat(df.Tₛ, df.θ), p)
100-element Vector{Float64}:
  4.233540174412755
 10.41149919818871
  ⋮
  1.746141124513421
  1.9599317903590014
```
"""
function DAMMfit(x::VecOrMat{<: Real}, Rₛ::Vector{Float64}, poro_val::Float64)
  lb = [0.0, 0.0, 0.0, 0.0] # params can't be negative
  p_ini = [1.0, 6.4, 1.0, 1.0] # initial parameters
  p_fact = (1e9, 1e1, 1e-6, 1e-4, 1.0, 1.0, 1.0) # factor of param, better fit
  DAMMfix = (x, p) -> DAMM(x, p_fact.*(p[1], p[2], p[3], p[4], poro_val, 0.02, 1.0))
  fit_bounds = curve_fit(DAMMfix, x, Rₛ, p_ini, lower=lb)
  p = fit_bounds.param # fitted param result
  fparam = p_fact[1:4].*p # multiply them by p_fact
  params = (fparam[1], fparam[2], fparam[3], fparam[4], poro_val, 0.02, 1.0)
  return params
end


text:
@testset "21.merge-two-sorted-lists.jl" begin
    list1 = ListNode{Int}([1, 2, 4])
    list2 = ListNode{Int}([1, 3, 4])
    expected = ListNode{Int}([1, 1, 2, 3, 4, 4])
    @test merge_two_lists(list1, list2) == expected
    @test merge_two_lists(nothing, ListNode(0)) == ListNode(0)
    @test merge_two_lists(nothing, nothing) === nothing
end

ext: jl


text:
@testset "arb_types" begin
    mag = Arblib.mag_struct()
    arf = Arblib.arf_struct()
    arb = Arblib.arb_struct()
    acb = Arblib.acb_struct()

    prec = 256
    for x in (arf, arb, acb)
        @test precision(x) == prec
        @test precision(Ptr{typeof(x)}()) == prec
        @test precision(typeof(x)) == prec
        @test precision(Ptr{typeof(x)}) == prec
    end
end


text:
"""
This module provides APIs for handling documentation metadata in modules.

The implementation is similar to how docstrings are handled in `Base` by the `Base.Docs`
module — a special variable is created in each module that has documentation metadata.

# Public API

* [`DocMeta.getdocmeta`](@ref)
* [`DocMeta.setdocmeta!`](@ref)

# Supported metadata

* `DocTestSetup`: contains the doctest setup code for doctests in the module.
"""
module DocMeta
using ..Documenter: Utilities
using DocStringExtensions

"The unique `Symbol` that is used to store the metadata dictionary in each module."
const META = gensym(:docmeta)

"List of modules that have the metadata dictionary added."
const METAMODULES = Module[]

"Type of the metadata dictionary."
const METATYPE = Dict{Symbol,Any}

"Dictionary of all valid metadata keys and their types."
const VALIDMETA = Dict{Symbol,Type}(:DocTestSetup => Union{Expr,Symbol})

"""
"""
function initdocmeta!(m::Module)
    if !isdefined(m, META)
        @debug "Creating documentation metadata dictionary (META=$META) in $m"
        Core.eval(m, :(const $META = $(METATYPE())))
        push!(METAMODULES, m)
    else
        @warn "Existing documentation metadata dictionary (META=$META) in $m. Ignoring."
    end
    return getfield(m, META)
end

"""
    getdocmeta(m::Module)

Returns the documentation metadata dictionary for the module `m`. The dictionary should be
considered immutable and assigning values to it is not well-defined. To set documentation
metadata values, [`DocMeta.setdocmeta!`](@ref) should be used instead.
"""
getdocmeta(m::Module) = isdefined(m, META) ? getfield(m, META) : METATYPE()

"""
    getdocmeta(m::Module, key::Symbol, default=nothing)

Return the `key` entry from the documentation metadata for module `m`, or `default` if the
value is unset.
"""
getdocmeta(m::Module, key::Symbol, default=nothing) = get(getdocmeta(m), key, default)

"""
    setdocmeta!(m::Module, key::Symbol, value; recursive=false, warn=true)

Set the documentation metadata value `key` for module `m` to `value`.

If `recursive` is set to `true`, it sets the same metadata value for all the submodules too.
If `warn` is `true`, it prints a warning when `key` already exists and is gets rewritten.
"""
function setdocmeta!(m::Module, key::Symbol, value; warn=true, recursive=false)
    key in keys(VALIDMETA) || throw(ArgumentError("Invalid metadata key\nValid keys are: $(join(keys(VALIDMETA), ", "))"))
    isa(value, VALIDMETA[key]) || throw(ArgumentError("Bad value type ($(typeof(value))) for metadata key $(key). Must be <: $(VALIDMETA[key])"))
    if recursive
        for mod in Utilities.submodules(m)
            setdocmeta!(mod, key, value; warn=warn, recursive=false)
        end
    else
        isdefined(m, META) || initdocmeta!(m)
        meta = getdocmeta(m)
        if warn && haskey(meta, key)
            @warn "$(key) already set for module $m. Overwriting."
        end
        meta[key] = value
    end
    return nothing
end

end


text:
#################################################################################################
# AIM 4: uniform distribution of null p-values
#################################################################################################
## Deps
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate("./") # <=
@everywhere using CpelNano
using StatsBase
using StatsPlots
using Distributions

## Source code
root_dir = dirname(dirname(pathof(CpelNano)))
include("$(root_dir)/calls/PaperSimulations/Aim-4/Aim-4-Source.jl")

## Constants

# IO
const aim_dir = "/Users/jordiabante/OneDrive - Johns Hopkins/CpelNano/Data/Simulations/Aim-4/" # <=
const blind_friend_col = ["#999999","#E69F00","#56B4E9","#009E73","#F0E442","#0072B2","#D55E00","#CC79A7"]

# Set up
@everywhere const M = 5
@everywhere const POBS = 1.0
@everywhere const LEVEL_STD_M = 2.0
@everywhere const LEVEL_STD_U = 2.0
@everywhere const LEVEL_MEAN_M = 84.0
@everywhere const LEVEL_MEAN_U = 80.0
@everywhere const EF_PERC_RNG = 0.1:0.1:0.1

## Default attributes
default(titlefont=(14, "arial"),guidefont=(16, "arial"),tickfont=(12, "arial"))

#################################################################################################
# Calls
#################################################################################################

# Args
matched = parse(Bool, ARGS[1])
n_sim_reps = parse(Int64, ARGS[2])
n_grp_reps = parse(Int64, ARGS[3])

# Generate null p-values
pvals_tmml, pvals_tnme, pvals_tcmd = run_null_sim(n_sim_reps, n_grp_reps, matched)

## Plot histogram

# Plot
plt = plt_null_pvals_hist(pvals_tmml, pvals_tnme, pvals_tcmd, matched)

# Store
if matched
    savefig("$(aim_dir)/Null-Pvals/Histogram-Null-Pvals-Matched.pdf")
else
    savefig("$(aim_dir)/Null-Pvals/Histogram-Null-Pvals-Unmatched.pdf")
end

## Plot ECDF

# Plot
plt = plt_null_pvals_ecdf(pvals_tmml, pvals_tnme, pvals_tcmd, matched)

# Store
if matched
    savefig("$(aim_dir)/Null-Pvals/ECDF-Null-Pvals-Matched.pdf")
else
    savefig("$(aim_dir)/Null-Pvals/ECDF-Null-Pvals-Unmatched.pdf")
end


text:
"""
    ScheduledOptim{T<:AbstractSchedule, O, F}
    ScheduledOptim(schedule::AbstractSchedule, opt, update_func)
    ScheduledOptim(schedule::AbstractSchedule, opt; update_func = (o, s) -> (o.eta = s))
    (::Type{<:AbstractSchedule})(opt; update_func = (o, s) -> (o.eta = s), kwargs...)

Wrap a `schedule` and `opt` together into a `ScheduledOptim`.
The `schedule` is iterated statefully using [`ScheduleIterator`](#) on every
call to [`Flux.apply!`](https://github.com/FluxML/Flux.jl/blob/master/src/optimise/optimisers.jl).
The `ScheduledOptim` can be used anywhere a Flux optimizer is used.

The keyword argument constructor sets `update_func` to schedule the learning rate of `opt`.

Instead of constructing `schedule` and `opt` separately,
you can use the `::Type{<:AbstractSchedule}` constructor (e.g. [`Exp`](#) below):
```julia
# create a Flux.Momentum optimizer
# where the learning rate is adjusted with
# Exp(λ = 0.1, γ = 0.8) decay schedule
opt = Exp(Momentum(); λ = 0.1, γ = 0.8)
```

# Arguments
- `schedule::AbstractSchedule`: the schedule to use
- `opt`: a Flux optimizer
- `update_func`: a mutating function of with inputs `(optim, param)`
                 that updates `optim` based on the current `param` value
"""
mutable struct ScheduledOptim{T<:AbstractSchedule, O, F}
    schedule::ScheduleIterator{T}
    optim::O
    update_func::F
end
ScheduledOptim(schedule::AbstractSchedule, opt, update_func) =
    ScheduledOptim(ScheduleIterator(schedule), opt, update_func)

function Flux.Optimise.apply!(opt::ScheduledOptim, x, Δ)
    # set param
    opt.update_func(opt.optim, next!(opt.schedule))

    # do normal apply
    return Flux.Optimise.apply!(opt.optim, x, Δ)
end

for Opt in (Descent, Momentum, Nesterov, RMSProp,
            ADAM, RADAM, AdaMax, OADAM, ADAGrad,
            ADADelta, AMSGrad, NADAM, AdaBelief)
    @eval begin
        ScheduledOptim(schedule::AbstractSchedule, opt::$Opt; update_func = (o, s) -> (o.eta = s)) =
            ScheduledOptim(schedule, opt, update_func)
        (::Type{T})(opt::$Opt; update_func = (o, s) -> (o.eta = s), kwargs...) where T<:AbstractSchedule =
            ScheduledOptim(T(kwargs...), opt, update_func)
    end
end


text:
# ComplexDual encompasses all complex dual types
const ComplexDual{T,V,N} = Complex{Dual{T,V,N}}

# Constructor for ComplexDual based on real Partials for dreal and dimag
@inline function ComplexDual{T,V,N}(z::Number,dr::Partials{N,V},di::Partials{N,V}) where {T,N,V<:Real}
    zr, zi = reim(z)
    return Dual{T}(zr,dr) + im*Dual{T}(zi,di)
end

@inline ComplexDual{T}(z::Number,u::Partials{N,V},v::Partials{N,V}) where {T,N,V} =
    ComplexDual{T,V,N}(z,u,v)

# Constructor for ComplexDual based on complex Partials for dz and dz*
@inline function ComplexDual{T}(z::Number,dz::Partials{N,V},dzstar::Partials{N,V}) where {T,N,V<:Complex}
    dr, di = reim_partials(dz,dzstar)
    return ComplexDual{T}(z,dr,di)
end

@inline ComplexDual{T}(z::Number,dz::Complex,dzstar::Complex) where {T} =
    ComplexDual{T}(z,Partials(tuple(dz)),Partials(tuple(dzstar)))

@inline Base.one(::Type{<:ComplexDual{T}},z::Number) where T = ComplexDual{T}(z,one(z),zero(z))

@inline valtype(::ComplexDual{T,V,N}) where {T,V,N} = V
@inline valtype(::Type{ComplexDual{T,V,N}}) where {T,V,N} = V
@inline valtype(::Array{<:ComplexDual{T,V,N}}) where {T,V,N} = V
@inline valtype(::Type{<:Array{<:ComplexDual{T,V,N}}}) where {T,V,N} = V

# ==== extensions of Partials ==== #
@inline Base.conj(v::Partials{N,V}) where {N,V} = Partials{N,V}(conj.(v.values))
@inline Base.real(v::Partials{N,Complex{V}}) where {N,V} = Partials{N,V}(real.(v.values))
@inline Base.imag(v::Partials{N,Complex{V}}) where {N,V} = Partials{N,V}(imag.(v.values))

@inline Base.:*(x::Complex, partials::Partials) = partials*x

@inline function Base.:*(partials::Partials, x::Complex)
        return Partials(ForwardDiff.scale_tuple(partials.values, x))
end

# extract the partial derivatives from a complex dual number
@inline function ForwardDiff.partials(d::ComplexDual)
    dr, di = reim(d)
    return partials(dr), partials(di)
end

# assemble the derivatives d/dz and d/dz* from the partials
@inline function _derivs(dvrx,dvry,dvix,dviy)
    return 0.5*(dvrx+dviy) + 0.5*im*(dvix-dvry), 0.5*(dvrx-dviy) + 0.5*im*(dvix+dvry)
end

# Given a complex dual number, extract Partials of dz and dz*
@inline function dz_partials(d::ComplexDual{T,V,N}) where {T,V,N}
    dr, di = partials(d)
    tmp  = [_derivs(dr[i], dr[i+1], di[i],di[i+1]) for i in 1:2:N]
    return Partials(tuple(first.(tmp)...)), Partials(tuple(last.(tmp)...))
end

# For the dz and dz* partials with respect to the ith element
@inline Base.@propagate_inbounds function dz_partials(d::ComplexDual, i)
    dr, di = partials(d)
    ix = 2*i.-1
    iy = 2*i
    return _derivs(dr[ix], dr[iy], di[ix],di[iy])
end

@inline Base.@propagate_inbounds dz_partials(::Type{T}, d::ComplexDual{T}, i...) where T = dz_partials(d, i...)

# For dealing with arrays of complex duals, return the dz and dz*
# partials as separate arrays
@inline dz_partials(d::AbstractVector{ComplexDual{T,V,N}}) where {T,V,N} =
  (tmp = dz_partials.(d); return first.(tmp), last.(tmp))

# For the dz and dz* partials with respect to the ith element
@inline Base.@propagate_inbounds dz_partials(d::AbstractVector{ComplexDual{T,V,N}},i) where {T,V,N} =
  (tmp = dz_partials.(d,i); return first.(tmp), last.(tmp))

@inline Base.@propagate_inbounds dz_partials(::Type{T}, d::AbstractVector{ComplexDual{T,V,N}}, i...) where {T,V,N} = dz_partials(d, i...)

# Given complex Partials for d/dz and d/dz*, return partials for dr/dx,y and di/dx,y
@inline function reim_partials(dz::Partials{N,V},dzstar::Partials{N,V}) where {T,V<:Complex,N}
    drx, dix = reim(dz+dzstar)
    dry, diy = reim(im*dz-im*dzstar)
    return _splice_xy(drx,dry), _splice_xy(dix,diy)
end

# Given partials for the Create separate partials for the real part and the imaginary part.
# These are organized as drx[1],dry[1],drx[2],dry[2],... and dix[1],diy[1],dix[2],diy[2],...
@inline function _splice_xy(dx::Partials{N},dy::Partials{N}) where {N}
    return Partials(tuple([ifelse(mod(i,2) == 0,dy[max(i÷2,1)],dx[(i+1)÷2]) for i in 1:2N]...))
end

# Return the value of a complex dual
@inline function value(d::ComplexDual)
    dr, di = reim(d)
    return value(dr)+im*value(di)
end


text:
abstract type AbstractChargeDriftModel{T <: SSDFloat} end
abstract type AbstractTemperatureModel{T <: SSDFloat} end

include("Vacuum/Vacuum.jl")
include("ADL/ADL.jl")


text:
# Utilities:
function draw_fullscreen(vao_id)
    glBindVertexArray(vao_id)
    glDrawArrays(GL_TRIANGLES, 0, 3)
    glBindVertexArray(0)
    return
end

struct PostprocessPrerender end

function (sp::PostprocessPrerender)()
    glDepthMask(GL_TRUE)
    glDisable(GL_DEPTH_TEST)
    glDisable(GL_BLEND)
    glDisable(GL_CULL_FACE)
    return
end

rcpframe(x) = 1f0 ./ Vec2f(x[1], x[2])

struct PostProcessor{F}
    robjs::Vector{RenderObject}
    render::F
end

function empty_postprocessor(args...; kwargs...)
    PostProcessor(RenderObject[], screen -> nothing)
end

function OIT_postprocessor(framebuffer)
    # Based on https://jcgt.org/published/0002/02/09/, see #1390
    # OIT setup
    shader = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/OIT_blend.frag")
    )
    data = Dict{Symbol, Any}(
        # :opaque_color => framebuffer[:color][2],
        :sum_color => framebuffer[:HDR_color][2],
        :prod_alpha => framebuffer[:OIT_weight][2],
    )
    pass = RenderObject(
        data, shader,
        () -> begin
            glDepthMask(GL_TRUE)
            glDisable(GL_DEPTH_TEST)
            glDisable(GL_CULL_FACE)
            glEnable(GL_BLEND)
            # shader computes:
            # src.rgb = sum_color / sum_weight * (1 - prod_alpha)
            # src.a = prod_alpha
            # blending: (assumes opaque.a = 1)
            # opaque.rgb = 1 * src.rgb + src.a * opaque.rgb
            # opaque.a   = 0 * src.a   + 1 * opaque.a
            glBlendFuncSeparate(GL_ONE, GL_SRC_ALPHA, GL_ZERO, GL_ONE)
        end,
        nothing
    )
    pass.postrenderfunction = () -> draw_fullscreen(pass.vertexarray.id)

    color_id = framebuffer[:color][1]
    full_render = screen -> begin
        fb = screen.framebuffer
        w, h = size(fb)

        # Blend transparent onto opaque
        glDrawBuffer(color_id)
        glViewport(0, 0, w, h)
        glDisable(GL_STENCIL_TEST)
        GLAbstraction.render(pass)
    end

    PostProcessor(RenderObject[pass], full_render)
end

function ssao_postprocessor(framebuffer)
    # Add missing buffers
    if !haskey(framebuffer, :position)
        glBindFramebuffer(GL_FRAMEBUFFER, framebuffer.id[1])
        position_buffer = Texture(
            Vec3f, size(framebuffer), minfilter = :nearest, x_repeat = :clamp_to_edge
        )
        pos_id = attach_colorbuffer!(framebuffer, :position, position_buffer)
        push!(framebuffer.render_buffer_ids, pos_id)
    end
    if !haskey(framebuffer, :normal)
        if !haskey(framebuffer, :HDR_color)
            glBindFramebuffer(GL_FRAMEBUFFER, framebuffer.id[1])
            normal_occlusion_buffer = Texture(
                Vec4{Float16}, size(framebuffer), minfilter = :nearest, x_repeat = :clamp_to_edge
            )
            normal_occ_id = attach_colorbuffer!(framebuffer, :normal_occlusion, normal_occlusion_buffer)
        else
            normal_occ_id = framebuffer[:HDR_color][1]
        end
        push!(framebuffer.render_buffer_ids, normal_occ_id)
    end

    # SSAO setup
    N_samples = 64
    lerp_min = 0.1f0
    lerp_max = 1.0f0
    kernel = map(1:N_samples) do i
        n = normalize([2.0rand() .- 1.0, 2.0rand() .- 1.0, rand()])
        scale = lerp_min + (lerp_max - lerp_min) * (i / N_samples)^2
        v = Vec3f(scale * rand() * n)
    end

    # compute occlusion
    shader1 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/SSAO.frag"),
        view = Dict(
            "N_samples" => "$N_samples"
        )
    )
    data1 = Dict{Symbol, Any}(
        :position_buffer => framebuffer[:position][2],
        :normal_occlusion_buffer => getfallback(framebuffer, :normal_occlusion, :HDR_color)[2],
        :kernel => kernel,
        :noise => Texture(
            [normalize(Vec2f(2.0rand(2) .- 1.0)) for _ in 1:4, __ in 1:4],
            minfilter = :nearest, x_repeat = :repeat
        ),
        :noise_scale => map(s -> Vec2f(s ./ 4.0), framebuffer.resolution),
        :projection => Observable(Mat4f(I)),
        :bias => 0.025f0,
        :radius => 0.5f0
    )
    pass1 = RenderObject(data1, shader1, PostprocessPrerender(), nothing)
    pass1.postrenderfunction = () -> draw_fullscreen(pass1.vertexarray.id)

    # blur occlusion and combine with color
    shader2 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/SSAO_blur.frag")
    )
    data2 = Dict{Symbol, Any}(
        :normal_occlusion => getfallback(framebuffer, :normal_occlusion, :HDR_color)[2],
        :color_texture => framebuffer[:color][2],
        :ids => framebuffer[:objectid][2],
        :inv_texel_size => lift(rcpframe, framebuffer.resolution),
        :blur_range => Int32(2)
    )
    pass2 = RenderObject(data2, shader2, PostprocessPrerender(), nothing)
    pass2.postrenderfunction = () -> draw_fullscreen(pass2.vertexarray.id)
    color_id = framebuffer[:color][1]

    full_render = screen -> begin
        fb = screen.framebuffer
        w, h = size(fb)

        # Setup rendering
        # SSAO - calculate occlusion
        glDrawBuffer(normal_occ_id)  # occlusion buffer
        glViewport(0, 0, w, h)
        glDisable(GL_STENCIL_TEST)
        glEnable(GL_SCISSOR_TEST)

        for (screenid, scene) in screen.screens
            # Select the area of one leaf scene
            # This should be per scene because projection may vary between
            # scenes. It should be a leaf scene to avoid repeatedly shading
            # the same region (though this is not guaranteed...)
            isempty(scene.children) || continue
            a = pixelarea(scene)[]
            glScissor(minimum(a)..., widths(a)...)
            # update uniforms
            data1[:projection] = scene.camera.projection[]
            data1[:bias] = scene.ssao.bias[]
            data1[:radius] = scene.ssao.radius[]
            GLAbstraction.render(pass1)
        end

        # SSAO - blur occlusion and apply to color
        glDrawBuffer(color_id)  # color buffer
        for (screenid, scene) in screen.screens
            # Select the area of one leaf scene
            isempty(scene.children) || continue
            a = pixelarea(scene)[]
            glScissor(minimum(a)..., widths(a)...)
            # update uniforms
            data2[:blur_range] = scene.ssao.blur
            GLAbstraction.render(pass2)
        end
        glDisable(GL_SCISSOR_TEST)
    end

    PostProcessor(RenderObject[pass1, pass2], full_render)
end

"""
    fxaa_postprocessor(framebuffer)

Returns a PostProcessor that handles fxaa.
"""
function fxaa_postprocessor(framebuffer)
    # Add missing buffers
    if !haskey(framebuffer, :color_luma)
        if !haskey(framebuffer, :HDR_color)
            glBindFramebuffer(GL_FRAMEBUFFER, framebuffer.id[1])
            color_luma_buffer = Texture(
                RGBA{N0f8}, size(framebuffer), minfilter=:linear, x_repeat=:clamp_to_edge
            )
            luma_id = attach_colorbuffer!(framebuffer, :color_luma, color_luma_buffer)
        else
            luma_id = framebuffer[:HDR_color][1]
        end
    end

    # calculate luma for FXAA
    shader1 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/postprocess.frag")
    )
    data1 = Dict{Symbol, Any}(
        :color_texture => framebuffer[:color][2],
        :object_ids => framebuffer[:objectid][2]
    )
    pass1 = RenderObject(data1, shader1, PostprocessPrerender(), nothing)
    pass1.postrenderfunction = () -> draw_fullscreen(pass1.vertexarray.id)

    # perform FXAA
    shader2 = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/fxaa.frag")
    )
    data2 = Dict{Symbol, Any}(
        :color_texture => getfallback(framebuffer, :color_luma, :HDR_color)[2],
        :RCPFrame => lift(rcpframe, framebuffer.resolution),
    )
    pass2 = RenderObject(data2, shader2, PostprocessPrerender(), nothing)
    pass2.postrenderfunction = () -> draw_fullscreen(pass2.vertexarray.id)

    color_id = framebuffer[:color][1]
    full_render = screen -> begin
        fb = screen.framebuffer
        w, h = size(fb)

        # FXAA - calculate LUMA
        glDrawBuffer(luma_id)
        glViewport(0, 0, w, h)
        # necessary with negative SSAO bias...
        glClearColor(1, 1, 1, 1)
        glClear(GL_COLOR_BUFFER_BIT)
        GLAbstraction.render(pass1)

        # FXAA - perform anti-aliasing
        glDrawBuffer(color_id)  # color buffer
        GLAbstraction.render(pass2)
    end

    PostProcessor(RenderObject[pass1, pass2], full_render)
end

"""
    to_screen_postprocessor(framebuffer)

Sets up a Postprocessor which copies the color buffer to the screen. Used as a
final step for displaying the screen.
"""
function to_screen_postprocessor(framebuffer)
    # draw color buffer
    shader = LazyShader(
        loadshader("postprocessing/fullscreen.vert"),
        loadshader("postprocessing/copy.frag")
    )
    data = Dict{Symbol, Any}(
        :color_texture => framebuffer[:color][2]
    )
    pass = RenderObject(data, shader, PostprocessPrerender(), nothing)
    pass.postrenderfunction = () -> draw_fullscreen(pass.vertexarray.id)

    full_render = screen -> begin
        fb = screen.framebuffer
        w, h = size(fb)

        # transfer everything to the screen
        glBindFramebuffer(GL_FRAMEBUFFER, 0)
        glViewport(0, 0, w, h)
        glClear(GL_COLOR_BUFFER_BIT)
        GLAbstraction.render(pass) # copy postprocess
    end

    PostProcessor(RenderObject[pass], full_render)
end


text:
solve(i::ElementaryPathInstance{T}, ::BellmanFordAlgorithm; kwargs...) where T = lp_dp(i; kwargs...)

function lp_dp(i::ElementaryPathInstance{T}) where T # I.e. Bellman-Ford algorithm. Assumption: no positive-cost cycle in the graph.
  V = Dict{T, Float64}()
  S = Dict{T, Vector{Edge{T}}}()

  # Initialise.
  for v in vertices(graph(i))
    V[v] = -Inf
    S[v] = Edge{T}[]
  end
  V[src(i)] = 0.0

  # Dynamic part.
  for _ in 1:ne(graph(i))
    changes = false # Stop the algorithm as soon as it no more makes progress.

    for e in edges(graph(i))
      u, v = src(e), dst(e)
      w = cost(i, u, v)

      # If using the solution to the currently explored subproblem would
      # lead to a cycle, skip it.
      if any(src(e) == v for e in S[u])
        continue
      end

      # Compute the maximum: is passing through w advantageous?
      if V[u] + w > V[v]
        V[v] = V[u] + w
        S[v] = vcat(S[u], Edge(u, v))

        changes = true
      end
    end

    if ! changes
      break
    end
  end

  # Checking existence of negative-cost cycle.
  for e in edges(graph(i))
    u, v = src(e), dst(e)
    w = cost(i, u, v)

    if V[u] + w > V[v]
      @warn("The graph contains a positive-cost cycle around edge $(u) -> $(v).")
    end
  end

  return ElementaryPathSolution(i, S[dst(i)], V, S)
end


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, MATLAB # Requires a matlab license and the DDS package to be in the right subfolder to run this code.

function solveBasicRelax(A, theK, eps)
    model=Model(Mosek.Optimizer)
    n=size(A,1)
    m=size(A,2)

    @variable(model, z[1:m]>=0.0)
    @constraint(model, z.<=1.0)
    @constraint(model, sum(z)<=theK)
    @variable(model, t)
    E = A * LinearAlgebra.diagm(z) * A'+eps*Matrix(I,n,n)
    @constraint(model, [t, 1, (E[i, j] for i in 1:n for j in 1:i)...] in MOI.LogDetConeTriangle(n))

    @objective(model, Max, t)

    @suppress optimize!(model)

    @show ofv_dual=objective_value(model)

    z_rounded=zeros(m)
    ind_max=sortperm(value.(z),rev=true)[1:theK]
    z_rounded[ind_max].=1.0
    @show ofv_primal=logdet(sum(z_rounded[i]*A[:,i]*A[:,i]' for i in 1:m)+eps*Matrix(I,n,n))
    @show ofv_dual=logdet(sum(value.(z[i])*A[:,i]*A[:,i]' for i in 1:m)+eps*Matrix(I,n,n))

    return value.(z), ofv_dual, z_rounded, ofv_primal
end

function solveBasicGD(A, k, eps)
    n=size(A,1)
    m=size(A,2)
    z0=zeros(m)
    ofv_primal=0.0
    for t=1:k
        z0, ofv_primal=getkplusone(A, k, m, n, z0, eps)
    end

    return z0, ofv_primal
end

function getkplusone(A, k, m, n, z0, eps) # Remark: introducing a +epsilon I term here so that the logdet term isn't negative
    candidateIndices=findall(z0.==0)
    iBest=-1
    ofvBest=-1e6
    Sigma_Gram_current=sum(z0[i]*A[:,i]*A[:,i]' for i in 1:m)
    for i in candidateIndices
        ofv_test=logdet(Sigma_Gram_current+A[:,i]*A[:,i]'+eps*Matrix(I,n,n)) # Use det rather than logdet here since logdet is undefined in the rank-deficient case... Note that greedy doesn't do anything when k<n anyway.
        if ofv_test>ofvBest
            ofvBest=ofv_test
            iBest=i
        end
    end
    z0[iBest]=1.0
    ofv_best=logdet(sum(z0[i]*A[:,i]*A[:,i]' for i in 1:m)+eps*Matrix(I,n,n))
    return z0, ofv_best
end

function solveQRERelax(Adata, theK, eps)

    n=size(Adata,1)
    m=size(Adata,2)
    # Call Matlab here, run via DDS solver, return results, round greedily and return an upper bound
    # eps_mat=mxarray(eps)
    # n_mat=mxarray(n)
    # m_mat=mxarray(m)
    # k_mat=mxarray(theK)
    # A_mat=mxarray(Adata)
    # Plan: build it in Matlab with random data, debug there, then move it here once it works!
    mat"""
    clear; A_mat=[[0.541299839388014,0.449426591726197,0.564473921782696,0.408709992443545,0.855835900703399,0.876145347529076,0.365236483082946,0.104519529950634,0.232254753050780,0.295421839485955,0.927064140870894,0.416765701751856,0.709240953533567,0.899426437101584,0.0990172329213107,0.241554770085929,0.955437531793689,0.228764772549632,0.517723514394718,0.490445656962020;0.868968871664997,0.659616768761805,0.540981282819544,0.947512642788601,0.724406283688824,0.420428701240150,0.598590500684565,0.00997765604969036,0.739832227063707,0.622028790367341,0.0877717606376327,0.280292866619113,0.641341360210975,0.452946115671495,0.570988144269953,0.912720162652781,0.142650346985794,0.423524871222756,0.0618247144123362,0.599441720175534;0.557047467058915,0.753208345950567,0.0689214476617044,0.919280969984633,0.199110301818659,0.487663607210372,0.668487146264776,0.0591533040753343,0.888992345152850,0.0475341941914539,0.332398561915437,0.598100317141488,0.174060860147872,0.0580261764965854,0.325875094193265,0.825734269699093,0.512563384780399,0.273596220894793,0.231367220579307,0.0902161331201864;0.0213981169649879,0.804739013812769,0.988430867357282,0.121243677375576,0.157285874413643,0.460325494339181,0.894564090918275,0.322652370507717,0.859807625912606,0.994609895249091,0.526180895421029,0.0364748567011317,0.0621519426336853,0.106269376654518,0.450492336647337,0.444545883918906,0.971925137586938,0.444565843866667,0.118485991561029,0.978224034127668;0.482680551891677,0.0291555821799541,0.251095440428348,0.591944834715854,0.370475519779316,0.515677249607505,0.0873355239070054,0.779476506848221,0.597060788610351,0.206764149130446,0.246642124058313,0.0636865030513627,0.406659420297776,0.998425291506397,0.577843795117458,0.982062563214574,0.648320621265659,0.627515046901424,0.0987800358735032,0.653018772578069;0.807985579141521,0.779817267477441,0.315472991527951,0.359652173164073,0.862272004954744,0.271994390193763,0.539010349313424,0.335478586307354,0.654752143480193,0.607363956107084,0.542899456586946,0.322882505170560,0.463062925909165,0.866327403667170,0.0748443075998355,0.578267561462169,0.614671003416509,0.534641253536576,0.890275455546976,0.461132126549800;0.736009366155118,0.567350498584087,0.300656948802795,0.719306932475474,0.684762628110137,0.231579838051303,0.428449936007297,0.619568612502970,0.915013993756726,0.347630325692841,0.780867850604077,0.0983785303848749,0.202674299437306,0.615205796593277,0.0573426693754819,0.234423496393930,0.469650371915959,0.385442159974304,0.0333890773864002,0.863758110483555;0.572304395272699,0.0761166268515728,0.0420390381653749,0.523551091500771,0.634197879289704,0.899534076872039,0.617153494459035,0.992885717464124,0.433183080639157,0.717737067921327,0.521882381409282,0.170044626268754,0.869549890049143,0.0269444844838657,0.300956435005654,0.810590442462569,0.577780697646748,0.873450099080692,0.839034991901461,0.262829072110785;0.00898485336788568,0.251621824309164,0.527901300708439,0.260843351982082,0.141323291186943,0.908697326644413,0.558875876956837,0.648005694952403,0.289760755842240,0.0279932951713195,0.931949933886006,0.371163798766144,0.597940369107138,0.322522635620462,0.521721248782275,0.451274164510932,0.911313227169849,0.300348734897046,0.507266724816296,0.823962942846020;0.718280463661178,0.133454528618960,0.256022333768320,0.493079306041446,0.0793022450665238,0.603642852234244,0.225850688199737,0.539775820318468,0.631882300545926,0.0668417611778476,0.147111980322787,0.0397616226212385,0.0230135370914409,0.463776233943276,0.561880235062444,0.249969169410276,0.376220283460707,0.400029715619638,0.113716441044558,0.328975305064346]];

    n_mat=10;
    m_mat=20;
    k_mat=7;
    eps_mat=1e-6;
    disp(n_mat)

    cvx_begin

        variable Y(n_mat,n_mat) semidefinite
        variable z(m_mat)

        eye(n_mat)-Y==semidefinite(n_mat);
        trace(Y)<=min(k_mat, n_mat);
        sum(z)==k_mat;
        z>=0;
        z<=ones(m_mat,1);

        minimize quantum_rel_entr(Y, A_mat*diag(z)*A_mat'+eps_mat*eye(n_mat))-(n_mat-trace(Y))*log(eps_mat)
    cvx_end

      $ofv_dual=-quantum_rel_entr(Y, A_mat*diag(z)*A_mat'+eps_mat*eye(n_mat))+(n_mat-trace(Y))*log(eps_mat);
      disp(ofv_dual)
      $z_relax=z
      [~,inds]=maxk(z, k_mat);
      $z_rounded=zeros(m_mat,1);
      $z_rounded(inds)=1.0;
      $ofv_primal=log(det(A_mat*diag(z_rounded)*A_mat'+eps_mat*eye(n_mat)));

    """

    @show ofv_primal, ofv_dual, z_relax, z_rounded

    return z_relax, ofv_dual, z_rounded, ofv_primal
end

max_stars_repo_path: D-opt Design/convexpenalties.jl

max_issues_repo_path: D-opt Design/convexpenalties.jl

max_forks_repo_path: D-opt Design/convexpenalties.jl


text:
export sympToGraph,  getColors,  getColorSetsFromSymp, contract

# Generate a Graph based on the given Symplectic Matrix
const Graph = AbstractMatrix{Bool}
const Colors = Vector
const Vertex = Int
const Colored = Bool
const NofModes = Int

function sympToGraph( S::Symp )::Graph
    n = nModes(S)
    return [ norm(S.S[2*i-1:2*i, 2*j-1:2*j]) > tolerance for i in 1:n, j in 1:n ]
end

function colorSuccessors(sucs, Cs::Colors)::Colors
    firstColor = Cs[sucs[1]]
    Cs[sucs] .= firstColor
    return Cs
end

function  getColors( G::Graph )::Colors
    colored = false
    n = isSquare(G) # = length(Cs)
    Cs = [1:n;]
    for v in 1:n
        sameColorVs = G[ findall(x->x == Cs[v], Cs ), :]
        sucs = findall(vi->any(sameColorVs[:, vi] ),  [1:n;] )
        nColors = length( Set(Cs[sucs]) )
        if nColors > 1
            Cs = colorSuccessors( sucs, Cs )
            colored = true
        end
    end
    return colored ? getColors(G) : Cs
end

# Is i a successor of j (both colored)?
function areColoredVeticesConnected(ci, cj, colorSets::Dict, G::Graph)::Bool
    vis = colorSets[ci]
    vjs = colorSets[cj]
    return any(G[vis, vjs])
end

getAllColors(Cs::Colors)::Vector = [Set(Cs)...]

function getColorSets(Cs::Colors)::Dict
    allColors = getAllColors(CS)
    return Dict(c => findall(x -> x == c, Cs) for c in allColors)
end

getColorSetsFromSymp = getColorSets ∘ getColors ∘ sympToGraph

function contract( G::Graph, Cs::Colors )
    colorSets = getColorSets(Cs)
    return [areColoredVeticesConnected(ci, cj, colorSets, G) for  ci in allColors, cj in allColors]
end


text:
#---
#hidden: true
#---

# This is a hidden julia card


text:
#=
Physics-inspired process simulation and data analysis related
to the albedo
=#
export surface_atmosphere_contributions, total_toa_albedo

"""
    surface_atmosphere_contributions(I, F_toa_⬆, F_s_⬆, F_s_⬇) → α_ATM, α_SFC
Calculate the atmospheric and surface **contributions** of the planetary albedo, so that
the TOA albedo is `α = α_ATM + α_SFC`, using the
simple 1-layer radiative transfer model by Donohoe & Battisti (2011) or G. Stephens (2015).
Stephens' formulas are incorrect and I have corrected them!
"""
function surface_atmosphere_contributions(I, F_toa_⬆, F_s_⬆, F_s_⬇)
    R = F_toa_⬆ ./ I   # planetary albedo (system reflectance)
    T = F_s_⬇ ./ I     # system transmisttance
    α = F_s_⬆ ./ F_s_⬇ # surface albedo

    # Formulas by Graeme, which are wrong!
    # τ = @. (1 - α*R) / (1 - (α^2) * (T^2)) # atmosphere transmittance
    # r = @. R - τ*α*T # atmospheric contribution to albedo

    # My calculations:
    r = @. (R - α*T^2) / (1 - (α*T)^2) # atmospheric contribution to albedo
    t = @. T*(1 - r*α)                 # atmospheric transmittance
    s = @. (α*t^2) / (1 - r*α)         # surface contribution to albedo
    return r, s
end

"""
    total_toa_albedo(a, s, t) = a + s*t^2/(1-a*s)
Combine given atmosphere albedo `a`, surface albedo `s` and atmosphere transmittance `t`
into a total top-of-the-atmosphere albedo `α` according to the model of Donohoe & Battisti (2011).
"""
total_toa_albedo(a, s, t) = a + s*t^2/(1-a*s)


text:
"""
    Variable{D <: AbstractDomain}
A structure containing the necessary information for a solver's variables: `name`, `domain`, and `constraints` it belongs.
```
struct Variable{D <: AbstractDomain}
    domain::D
    constraints::Indices{Int}
end
```
"""
mutable struct Variable{D <: AbstractDomain}
    domain::D
    constraints::Indices{Int}
end

function Variable(D, x::Variable{D2}) where {D2 <: AbstractDomain}
    return Variable{D}(x.domain, x.constraints)
end

# Methods: lazy forwarding from ConstraintDomains.domain.jl
@forward Variable.domain Base.length, Base.rand, Base.delete!, Base.isempty
# @forward Variable.domain get_domain, domain_size

add!(x::Variable, value) = ConstraintDomains.add!(x.domain, value)
get_domain(x::Variable) = ConstraintDomains.get_domain(x.domain)
domain_size(x::Variable) = ConstraintDomains.domain_size(x.domain)

"""
    _get_constraints(x::Variable)
Access the list of `constraints` of `x`.
"""
_get_constraints(x::Variable) = x.constraints

"""
    _add_to_constraint!(x::Variable, id)
Add a constraint `id` to the list of contraints of `x`.
"""
_add_to_constraint!(x::Variable, id) = set!(_get_constraints(x), id)

"""
    _delete_from_constraint!(x::Variable, id)
Delete a constraint `id` from the list of contraints of `x`.
"""
_delete_from_constraint!(x::Variable, id) = delete!(x.constraints, id)

"""
    _constriction(x::Variable)
Return the `cosntriction` of `x`, i.e. the number of constraints restricting `x`.
"""
_constriction(x::Variable) = length(x.constraints)

"""
    x::Variable ∈ constraint
    value ∈ x::Variable
Check if a variable `x` is restricted by a `constraint::Int`, or if a `value` belongs to the domain of `x`.
"""
Base.in(x::Variable, constraint) = constraint ∈ x.constraints
Base.in(value, x::Variable) = value ∈ x.domain

"""
    variable(values::AbstractVector{T}, name::AbstractString; domain = :set) where T <: Number
    variable(domain::AbstractDomain, name::AbstractString) where D <: AbstractDomain
Construct a variable with discrete domain. See the `domain` method for other options.

```julia
d = domain([1,2,3,4], types = :indices)
x1 = variable(d, "x1")
x2 = variable([-89,56,28], "x2", domain = :indices)
```
"""
variable() = Variable(domain(), Indices{Int}())
variable(domain::D) where {D <: AbstractDomain} = Variable(domain, Indices{Int}())
variable(vals) = isempty(vals) ? variable() : variable(domain(vals))


text:
abstract type RBFParameterSelection end

struct RBFParLinear <: RBFParameterSelection
    c
end
struct RBFParSquareRoot <: RBFParameterSelection
    c
end
struct RBFParConstant <: RBFParameterSelection
    c
end

RBFparameter(p::RBFParLinear, n::Int) = p.c*n
RBFparameter(p::RBFParSquareRoot, n::Int) = p.c*sqrt(convert(typeof(p.c),n))
RBFparameter(p::RBFParConstant, n::Int) = p.c

name(p::RBFParLinear) = "$(@sprintf("%.2f", p.c)) * N"
name(p::RBFParSquareRoot) = "$(@sprintf("%.2f", p.c)) * sqrt(N)"
name(p::RBFParConstant) = "$(@sprintf("%.2f", p.c))"

rbfselect(rbf::RBF, n, select) = similar(rbf, RBFparameter(select, n))

struct RBFPlatform{R <: RBF} <: FrameFun.Platform
    rbf     ::  R
    domain  ::  Domain
    box     ::  Domain
    parameterstyle
    gridgen
end

default_selection(rbf::SmoothRBF{T}) where {T} = RBFParLinear(one(T))
default_selection(rbf::PiecewiseSmoothRBF) = RBFParConstant(rbf.p)

default_gridgenerator(rbf::RBF, domain::Domain{T}, box::Domain{T}) where {T} =
    EquispacedGrids{T}(box)

RBFPlatform(rbf::SmoothRBF{T}; domain = zero(T)..one(T),
        box = 2*boundingbox(domain),
        select = default_select(rbf),
        gridgen = default_gridgenerator(rbf, domain, box)) where {T} =
    RBFPlatform(rbf, domain, box, select, gridgen)

dictionary(pf::RBFPlatform, ppar) = dictionary(pf, ppar, pf.rbf, pf.domain, pf.box, pf.parameterstyle, pf.gridgen)

function dictionary(pf::RBFPlatform, n, rbf, domain, box, paramstyle, gridgen)
    rbfargs = RBFparameter(paramstyle, n)
    # centers = equispaced_centers(pf.domain, pf.box, n)
    centers = collect(generate(gridgen, n))
    RBFDictionary(centers[:], similar(rbf, rbfargs...), domain)
end

SamplingStyle(pf::RBFPlatform) = OversamplingStyle()
SolverStyle(pf::RBFPlatform) = DirectStyle()

# function equispaced_centers(domain::Domain, box, n)
#     g = equispaced_grid(box, n)
#     collect(g)
# end
#
# equispaced_grid(box::AbstractInterval{T}, n) where {T} =
#     EquispacedGrid{T}(n, infimum(box), supremum(box))
#
# equispaced_grid(box::ProductDomain, n) =
#     ProductGrid(equispaced_grid.(elements(box), n)...)
#
# periodic_equispaced_grid(box::AbstractInterval{T}, n) where {T} =
#     EquispacedGrid{T}(n, infimum(box), supremum(box))
#
# periodic_equispaced_grid(box::ProductDomain, n) =
#     ProductGrid(map(periodic_equispaced_grid, elements(box), n)...)


text:
# This file was generated, do not modify it. # hide
@assert length(G) == 20_000

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
# An IntervalCollection is an efficiently stored and indexed set of annotated genomic intervals.
# It looks something like this.
#
#                                      ┌─────┐
#                                      │trees│
#                                      └─────┘
#                                         │
#                              ┌──────────┴──────────┬────────────┐
#                              ▼                     ▼            │
#  Each sequence has       ┌──────┐              ┌──────┐         ▼
#    an associated         │ chr1 │              │ chr2 │
#    IntervalTree     ┌────┴──────┴────┐    ┌────┴──────┴────┐    ...
#                     │ IntervalTree 1 │    │ IntervalTree 2 │
#                     └────────────────┘    └────────────────┘
#                              │
#                    ┌─────────┴─────────────────────────┬────────────────────────┐
#                    │                                   │                        │
#                    ▼                                   ▼                        ▼
#   ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓  ┏━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
#   ┃ Interval{T}(10000, 20000, ...) ┃  ┃ Interval{T}(35000, 40000, ...) ┃      ...
#   ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛  ┗━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┛
#
#
#
#
#                               ┌────────────────┐
#       ordered_trees holds     │ IntervalTree 1 │
#         an array of the       ├────────────────┤
#       some IntervalTrees,     │ IntervalTree 2 │
#           ordered by          └────────────────┘
#       chromosome for fast
#       ordered iteration.             ...
#

# Aliases for types of IntervalTrees.jl (IC: Interval Collection).
const ICTree{T}                               = IntervalTrees.IntervalBTree{Int64,Interval{T},64}
const ICTreeIteratorState{T}                  = IntervalTrees.IntervalBTreeIteratorState{Int64,Interval{T},64}
const ICTreeIntersection{T}                   = IntervalTrees.Intersection{Int64,Interval{T},64}
const ICTreeIntersectionIterator{F,S,T}       = IntervalTrees.IntersectionIterator{F,Int64,Interval{S},64,Interval{T},64}
const ICTreeIntervalIntersectionIterator{F,T} = IntervalTrees.IntervalIntersectionIterator{F, Int64,Interval{T},64}

"An IntervalCollection is an efficiently stored and indexed set of annotated genomic intervals."
mutable struct IntervalCollection{T}
    # Sequence name mapped to IntervalTree, which in turn maps intervals to a list of metadata.
    trees::Dict{String,ICTree{T}}

    # Keep track of the number of stored intervals.
    length::Int

    # A vector of values(trees) sorted on sequence name.
    # This is used to iterate intervals as efficiently as possible, but is only updated as needed, indicated by the ordered_trees_outdated flag.
    ordered_trees::Vector{ICTree{T}}
    ordered_trees_outdated::Bool

    "Empty initaialzation."
    function IntervalCollection{T}() where T
        return new{T}(Dict{String,ICTree{T}}(), 0, ICTree{T}[], false)
    end

    "Bulk insertion."
    function IntervalCollection{T}(intervals::AbstractVector{Interval{T}}, sort::Bool=false) where T
        if sort
            sort!(intervals)
        else
            if !issorted(intervals)
                error("Intervals must be sorted, or `sort=true` set, to construct an IntervalCollection")
            end
        end

        n = length(intervals)
        trees = Dict{String,ICTree{T}}()
        i = 1
        while i <= n
            j = i
            while j <= n && seqname(intervals[i]) == seqname(intervals[j])
                j += 1
            end
            trees[seqname(intervals[i])] = ICTree{T}(view(intervals, i:j-1))
            i = j
        end
        return new{T}(trees, n, ICTree{T}[], true)
    end
end

"""
    IntervalCollection(intervals::AbstractVector{Interval{T}}, sort::Bool=false) where T
Shorthand constructor.
"""
function IntervalCollection(intervals::AbstractVector{Interval{T}}, sort::Bool=false) where T
    return IntervalCollection{T}(intervals, sort)
end

"""
    IntervalCollection{T}(data, sort::Bool=false) where T
Constructor that offers conversion through collection.
"""
function IntervalCollection{T}(data, sort::Bool=false) where T
    return IntervalCollection(collect(Interval{T}, data), sort)
end

"""
    IntervalCollection(data, sort::Bool=false)
Constructor that guesses metadatatype, and offers conversion through collection.
"""
function IntervalCollection(data, sort::Bool=false)
    return IntervalCollection(collect(Interval{metadatatype(data)}, data), sort)
end

function update_ordered_trees!(ic::IntervalCollection{T}) where T
    if ic.ordered_trees_outdated
        ic.ordered_trees = collect(ICTree{T}, values(ic.trees))
        p = sortperm(collect(AbstractString, keys(ic.trees)), lt = isless)
        ic.ordered_trees = ic.ordered_trees[p]
        ic.ordered_trees_outdated = false
    end
end

function Base.push!(ic::IntervalCollection{T}, i::Interval{T}) where T
    tree = get!(ic.trees, seqname(i)) do
        # Setup empty tree for new seqname key.
        ic.ordered_trees_outdated = true
        return ICTree{T}()
    end
    push!(tree, i)
    ic.length += 1
    return ic
end

function Base.show(io::IO, ic::IntervalCollection{T}) where T
    n_entries = length(ic)
    println(io, "IntervalCollection{$(T)} with $(n_entries) intervals:")
    if n_entries > 0
        for (k, i) in enumerate(ic)
            if k > 8
                break
            end
            println(IOContext(io, :compact=>true), "  ", i)
        end
        if n_entries > 8
            print(io, "  ⋮")
        end
    end
end

function Base.length(ic::IntervalCollection)
    return ic.length
end

function Base.eltype(::Type{IntervalCollection{T}}) where T
    return Interval{T}
end

function Base.:(==)(a::IntervalCollection, b::IntervalCollection)
    if length(a) != length(b)
        return false
    end

    # Return false if any of the intervals are not equivalent.
    return all(zip(a,b)) do (ia,ib)
        ia == ib
    end
end

# Iterators
# ---------

#=
mutable struct IntervalCollectionIteratorState{T}
    i::Int # index into ordered_trees
    tree_state::ICTreeIteratorState{T}

    function IntervalCollectionIteratorState{T}(i::Int) where T
        return new{T}(i)
    end

    function IntervalCollectionIteratorState{T}(i::Int, tree_state) where T
        return new{T}(i, tree_state)
    end
end

function iterinitstate(ic::IntervalCollection{T}) where T
    update_ordered_trees!(ic)
    i = 1
    while i <= length(ic.ordered_trees)
        tree_state = IntervalTrees.iterinitstate(ic.ordered_trees[i])
        if !(tree_state.leaf === nothing || isempty(tree_state.leaf))
            return IntervalCollectionIteratorState{T}(i, tree_state)
        end
        i += 1
    end
    return IntervalCollectionIteratorState{T}(i)
end

function Base.start(ic::IntervalCollection{T}) where T
    update_ordered_trees!(ic)
    i = 1
    while i <= length(ic.ordered_trees)
        tree_state = start(ic.ordered_trees[i])
        if !done(ic.ordered_trees[i], tree_state)
            return IntervalCollectionIteratorState{T}(i, tree_state)
        end
        i += 1
    end
    return IntervalCollectionIteratorState{T}(i)
end
=#

#=
function Base.iterate(ic::IntervalCollection, state=iterinitstate(ic))
    i = state.i
    treeit = iterate(ic.ordered_trees[i], state.tree_state)
    if treeit === nothing
        i += 1
        while i <= length(ic.ordered_trees)
            tree_state = IntervalTrees.iterinitstate(ic.ordered_trees[i])
            if !(tree_state.leaf === nothing || isempty(tree_state.leaf))
                break
            end
            i += 1
        end
    end
    state.i, state.tree_state = i, tree_state
    return value, state
end
=#

function iterprep(ic::IntervalCollection)
    update_ordered_trees!(ic)
    return ()
end

# State is a tuple:
# (ot iteration state, current ot element, current ot element state)
# "ot" is shorthand for "ordered_trees"
# This iterate method basically works like the Iterators.Flatten iterator.
@propagate_inbounds function Base.iterate(ic::IntervalCollection, state = iterprep(ic))
    if state !== ()
        # Iterate over each interval in current ordered tree.
        tree_it = iterate(Base.tail(state)...)
        tree_it !== nothing && return (tree_it[1], (state[1], state[2], tree_it[2]))
    end
    # Iterate over ordered_trees.
    ot_it = (state === () ? iterate(ic.ordered_trees) : iterate(ic.ordered_trees, state[1]))
    ot_it === nothing && return nothing
    iterate(ic, (ot_it[2], ot_it[1]))
end

#=
function Base.next(ic::IntervalCollection, state)
    i = state.i
    value, tree_state = next(ic.ordered_trees[i], state.tree_state)
    if done(ic.ordered_trees[i], tree_state)
        i += 1
        while i <= length(ic.ordered_trees)
            tree_state = start(ic.ordered_trees[i])
            if !done(ic.ordered_trees[i], tree_state)
                break
            end
            i += 1
        end
    end
    state.i, state.tree_state = i, tree_state
    return value, state
end

function Base.done(ic::IntervalCollection, state)
    return state.i > length(ic.ordered_trees)
end
=#

# Filter predicates
# -----------------

true_cmp(a, b) = true

eachoverlap_equal_filter(a, b) = first(a) == first(b) && last(a) == last(b)

# TODO: other prebaked filter predicates

"""
Find a the first interval with matching start and end points.

Returns that interval, or 'nothing' if no interval was found.
"""
function Base.findfirst(a::IntervalCollection, b::Interval; filter=true_cmp)
    if !haskey(a.trees, seqname(b))
        return nothing
    end

    return findfirst(a.trees[seqname(b)], b, filter)
end

# Overlaps
# --------

function eachoverlap(a::IntervalCollection{T}, query::Interval; filter::F = true_cmp) where {F,T}
    if haskey(a.trees, seqname(query))
        return ICTreeIntervalIntersectionIterator{F,T}(filter, ICTreeIntersection{T}(), a.trees[seqname(query)], query)
    end

    return ICTreeIntervalIntersectionIterator{F,T}(filter, ICTreeIntersection{T}(), ICTree{T}(), query)
end

function eachoverlap(a::IntervalCollection, b::IntervalCollection; filter = true_cmp)
    seqnames = collect(AbstractString, keys(a.trees) ∩ keys(b.trees))
    sort!(seqnames, lt = isless)
    a_trees = [a.trees[seqname] for seqname in seqnames]
    b_trees = [b.trees[seqname] for seqname in seqnames]
    return IntersectIterator(filter, a_trees, b_trees)
end

struct IntersectIterator{F, S, T}
    filter::F
    a_trees::Vector{ICTree{S}}
    b_trees::Vector{ICTree{T}}
end

#=
mutable struct IntersectIteratorState{F,S,T}
    i::Int  # index into a_trees/b_trees.
    intersect_iterator::ICTreeIntersectionIterator{F,S,T}

    function IntersectIteratorState{F,S,T}(i) where {F,S,T}
        return new{F,S,T}(i)
    end

    function IntersectIteratorState{F,S,T}(i, iter) where {F,S,T}
        return new{F,S,T}(i, iter)
    end
end
=#

function Base.eltype(::Type{IntersectIterator{F,S,T}}) where {F,S,T}
    return Tuple{Interval{S},Interval{T}}
end

function Base.IteratorSize(::Type{IntersectIterator{F,S,T}}) where {F,S,T}
    return Base.SizeUnknown()
end

#=
function Base.start(it::IntersectIterator{F,S,T}) where {F,S,T}
    i = 1
    while i <= length(it.a_trees)
        intersect_iterator = intersect(it.a_trees[i], it.b_trees[i], it.filter)
        intersect_iterator_state = start(intersect_iterator)
        if !done(intersect_iterator, intersect_iterator_state)
            return IntersectIteratorState{F,S,T}(i, intersect_iterator)
        end
        i += 1
    end
    return IntersectIteratorState{F,S,T}(i)
end
=#

# State is a tuple:
# (tree pair iteration state, current intersect iterator, current intersect iterator state)
function Base.iterate(it::IntersectIterator{F, S, T}, state = ()) where {F,S,T}
    if state !== ()
        # Iterate over each intersection between two trees.
        isect = iterate(Base.tail(state)...)
        isect !== nothing && return (isect[1], (state[1], state[2], isect[2]))
    end
    # Iterate over tree pairs.
    treeA = (state === () ? iterate(it.a_trees) : iterate(it.a_trees, state[1]))
    treeB = (state === () ? iterate(it.b_trees) : iterate(it.b_trees, state[1]))
    treeA === nothing && return nothing
    intersect_iterator = intersect(treeA[1], treeB[1], it.filter)
    iterate(it, (treeA[2], intersect_iterator))
end

#=
function Base.next(it::IntersectIterator{F, S, T}, state) where {F,S,T}
    i, intersect_iterator = state.i, state.intersect_iterator
    value, intersect_iterator_state = next(intersect_iterator, nothing)
    if done(intersect_iterator, intersect_iterator_state)
        i += 1
        while i <= length(it.a_trees)
            intersect_iterator = intersect(it.a_trees[i], it.b_trees[i],
                                           it.filter)
            intersect_iterator_state = start(intersect_iterator)
            if !done(intersect_iterator, intersect_iterator_state)
                break
            end
            i += 1
        end
    end
    state.i, state.intersect_iterator = i, intersect_iterator
    return value, state
end

function Base.done(it::IntersectIterator{F, S, T}, state) where {F, S, T}
    return state.i > length(it.a_trees)
end
=#

function eachoverlap(a, b::IntervalCollection; filter=true_cmp)
    return IntervalCollectionStreamIterator(filter, a, b)
end

struct IntervalCollectionStreamIterator{F,S,T}
    filter::F
    stream::S
    collection::IntervalCollection{T}
end

function Base.eltype(::Type{IntervalCollectionStreamIterator{F,S,T}}) where {F,S,T}
    return Tuple{Interval{metadatatype(S)},Interval{T}}
end

function Base.IteratorSize(::Type{IntervalCollectionStreamIterator{F,S,T}}) where {F,S,T}
    return Base.SizeUnknown()
end

#= ### The old iteration protocol prior to version 0.7 / 1.0 of julia.

mutable struct IntervalCollectionStreamIteratorState{F,Ta,Tb,U}
    intersection::ICTreeIntersection{Tb}
    stream_value::Interval{Ta}
    stream_state::U

    function IntervalCollectionStreamIteratorState{F,Ta,Tb,U}(intersection, stream_value, stream_state) where {F,Ta,Tb,U}
        return new{F,Ta,Tb,U}(intersection, stream_value, stream_state)
    end

    function IntervalCollectionStreamIteratorState{F,Ta,Tb,U}() where {F,Ta,Tb,U}
        return new{F,Ta,Tb,U}(ICTreeIntersection{Tb}())
    end
end

# This mostly follows from SuccessiveTreeIntersectionIterator in IntervalTrees
function Base.start(it::IntervalCollectionStreamIterator{F,S,T}) where {F,S,T}
    stream_state = start(it.stream)
    intersection = ICTreeIntersection{T}()
    while !done(it.stream, stream_state)
        stream_value, stream_state = next(it.stream, stream_state)
        if haskey(it.collection.trees, seqname(stream_value))
            tree = it.collection.trees[seqname(stream_value)]
            IntervalTrees.firstintersection!(tree, stream_value, nothing, intersection, it.filter)
            if intersection.index != 0
                return IntervalCollectionStreamIteratorState{F,T,metadatatype(it.stream),typeof(stream_state)}(intersection, stream_value, stream_state)
            end
        end
    end
    return IntervalCollectionStreamIteratorState{S,metadatatype(it.stream),typeof(stream_state)}()
end

function Base.next(it::IntervalCollectionStreamIterator{F,S,T}, state) where {F,S,T}
    intersection = state.intersection
    entry = intersection.node.entries[intersection.index]
    return_value = (state.stream_value, entry)
    IntervalTrees.nextintersection!(intersection.node, intersection.index, state.stream_value, intersection, it.filter)
    while intersection.index == 0 && !done(it.stream, state.stream_state)
        state.stream_value, state.stream_state = next(it.stream, state.stream_state)
        if haskey(it.b.trees, seqname(state.stream_value))
            tree = it.b.trees[seqname(state.stream_value)]
            IntervalTrees.firstintersection!(tree, state.stream_value, nothing, intersection, it.filter)
        end
    end
    return return_value, state
end

function Base.done(it::IntervalCollectionStreamIterator, state)
    return state.intersection.index == 0
end
=#

# New julia 0.7 / 1.0 iteration protocol for collection stream iterator.
# State is a tuple:
# (current_query, stream_state, intersection_object)
function Base.iterate(it::IntervalCollectionStreamIterator{F,S,T}, state = ()) where {F,S,T}
    # If first iteration, make empty intersection, otherwise get it from the state.
    intersection = (state !== () ? state[3] : ICTreeIntersection{T}())

    # If this is not the first iteration, and there is an available intersection for the current query, return it and search for the next intersection.
    if state !== () && intersection.index != 0
        entry = intersection.node.entries[intersection.index]
        return_value = (state[1], entry)
        IntervalTrees.nextintersection!(intersection.node, intersection.index, state[1], intersection, it.filter)
        return return_value, (state[1], state[2], intersection)
    end

    # If code reaches this point, there is no valid intersection to return for the current query, so we get the next query and start looking for intersections.
    while intersection.index == 0
        # Get a new query from the stream, and its first intersection in the collection.
        stream_it = (state === () ? iterate(it.stream) : iterate(it.stream, state[2]))
        stream_it === nothing && return nothing # You have reached the end of the stream, stop iterating.
        stream_value = stream_it[1]
        if haskey(it.collection.trees, seqname(stream_value))
            tree = it.collection.trees[seqname(stream_value)]
            IntervalTrees.firstintersection!(tree, stream_value, nothing, intersection, it.filter)
            iterate_result = iterate(it, (stream_value, stream_it[2], intersection))
            if iterate_result !== nothing
                return iterate_result
            end
        end
    end
end


text:
module Yfjd

using Dates
using DataFrames
using TimeSeries
using HTTP
using JSON

export get_symbols_df, get_symbols_df_raw
export get_symbols, get_symbols_ta, get_symbols_ta_clo, get_symbols_ta_adj

# deleterows! becomes delete!  since DataFrames v0.20 (breaking)
if  isdefined(DataFrames,:deleterows!)
  Ydelrows! = DataFrames.deleterows!
elseif isdefined(DataFrames,:delete!)
  Ydelrows! = DataFrames.delete!
else
  error("Cant find \"Dataframes.delete[rows]!\" ")
end

"""

# Download raw stock market data from Yahoo Finance
# Returns a DataFrame

    get_symbols_df_raw(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols_df_raw("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols_df_raw("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols_df_raw("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols_df_raw("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols_df_raw(symbol::String, from::String, to::String)

    from = string(Integer(datetime2unix(DateTime(from * "T12:00:00"))))
    to = string(Integer(datetime2unix(DateTime(to * "T12:00:00"))))
    parse(Int, from) > parse(Int, to) ? error("in get_symbols: to must be older than from") : nothing

    url = "https://query1.finance.yahoo.com/v7/finance/chart/$symbol?&interval=1d&period1=$from&period2=$to"

    response = HTTP.get(url, cookies = true)
    body = JSON.parse(String(response.body))["chart"]["result"][1]
    values = body["indicators"]["quote"][1]

    df = DataFrame(
        Open = values["open"],
        High =  values["high"],
        Low = values["low"],
        Close = values["close"],
        Volume = values["volume"],
        Adjusted = body["indicators"]["adjclose"][1]["adjclose"],
        Time = Dates.Date.(unix2datetime.(body["timestamp"]))
    )

    return df
end

"""

# Download "cooked" stock market data from Yahoo Finance
# Returns a DataFrame
# Cooking is converting to Float64, rounded to 2 decimals
    get_symbols_df(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols_df("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols_df("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols_df("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols_df("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols_df(symbol::String, from::String, to::String)
    df = get_symbols_df_raw(symbol, from, to)
    Ydelrows!(df, isnothing.(df).Close)

	df.Open     = Float64.(df.Open)
	df.High     = Float64.(df.High)
	df.Low      = Float64.(df.Low)
	df.Close    = Float64.(df.Close)
	df.Volume   =   Int64.(df.Volume)
	df.Adjusted = Float64.(df.Adjusted)

	df.Open     = round.(df.Open    ,digits=2)
	df.High     = round.(df.High    ,digits=2)
	df.Low      = round.(df.Low     ,digits=2)
	df.Close    = round.(df.Close   ,digits=2)
	df.Adjusted = round.(df.Adjusted,digits=2)

    return df
end

"""

# Download stock market data (:Close) from Yahoo Finance
# Returns a TimeArray, so omit times without quotation

    get_symbols_ta_clo(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols_ta_clo("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols_ta_clo("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols_ta_clo("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols_ta_clo("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols_ta_clo(symbol::String, from::String, to::String)
    df = get_symbols_df(symbol, from, to)
    Ydelrows!(df, isnothing.(df).Close)
    x = TimeArray(df, timestamp = :Time)
    return x
end

"""

# Download stock market data (:Adjusted) from Yahoo Finance
# Returns a TimeArray, so omit times without quotation

    get_symbols_ta_adj(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols_ta_adj("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols_ta_adj("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols_ta_adj("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols_ta_adj("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols_ta_adj(symbol::String, from::String, to::String)
    df = get_symbols_df(symbol, from, to)
	Ydelrows!(df, isnothing.(df).Adjusted)
    x = TimeArray(df, timestamp = :Time)
    return x
end

"""

# Download stock market data (:Adjusted) from Yahoo Finance
# Returns a TimeArray, so omit times without quotation

    get_symbols_ta(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols_ta("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols_ta("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols_ta("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols_ta("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols_ta(symbol::String, from::String, to::String)
    return get_symbols_ta_adj(symbol, from, to)
end

"""

# Download stock market data (:Adjusted)from Yahoo Finance
# Returns a TimeArray, so omit times without quotation

    get_symbols(symbol, from, to)

## Arguments

    - `symbol`: Market symbol, e.g. "AAPL" or "GOOG"
    - `from` : Date character, e.g. "2018-12-13" in the form of "YYYY-MM-DD"
    - `to` : Date character, e.g. "2019-12-13" in the form of "YYYY-MM-DD"

## Examples

```jldoctest
julia> get_symbols("GOOG", "2018-12-26", "2019-12-20")
julia> get_symbols("^GDAXI", "2018-12-26", "2019-12-20")
julia> get_symbols("EURUSD=X", "2018-12-26", "2019-12-20")
julia> get_symbols("BTC-USD", "2018-12-26", "2019-02-20")
```
"""
function get_symbols(symbol::String, from::String, to::String)
    return get_symbols_ta_adj(symbol, from, to)
end

end


text:
using Test

function parseInput(input)
    chars = split(input, ",")
    chars = parse.(Int, chars)
    return chars
end

function part1(codes)
    i = 1
    while codes[i] != 99
        operation = codes[i]
        pos1 = codes[i+1] + 1
        pos2 = codes[i+2] + 1
        dest = codes[i+3] + 1
        # pos1, pos2, dest = codes[i+1:i+3]
        if operation == 1
            codes[dest] = codes[pos1] + codes[pos2]
        elseif operation == 2
            codes[dest] = codes[pos1] * codes[pos2]
        end
        i += 4
    end
    return codes
end

function run(s)
    for i in 1:99
        for j in 1:99
            input = parseInput(split(s, "\n")[1])
            input[2] = i
            input[3] = j
            output = part1(input)
            if output[1] == 19690720
                return 100 * i + j
            end
        end
    end
end

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time*1000)")
    println(res)
end

main()


text:
# InferenceErrorReport
# ====================

# abstract stack trace
# --------------------

"""
    Signature

Represents an expression signature.
`print_signature` implements a frontend functionality to show this type.
"""
struct Signature
    _sig::Vector{Any}
end

# define equality functions that avoid dynamic dispatches
function Base.:(==)(sig1::Signature, sig2::Signature)
    sig1 = sig1._sig
    sig2 = sig2._sig
    length(sig1) == length(sig2) || return false
    for (a1, a2) in zip(sig1, sig2)
        a1 === a2 || return false
    end
    return true
end
function Base.hash(sig::Signature, h::UInt)
    for a in sig._sig
        h = @invoke hash(a::Any, h::UInt)
    end
    return h
end
@inline Base.iterate(sig::Signature, state...) = iterate(sig._sig, state...)

"""
    VirtualFrame

Stack information representing virtual execution context:
- `file::Symbol`: the path to the file containing the virtual execution context
- `line::Int`: the line number in the file containing the virtual execution context
- [`sig::Signature`](@ref Signature): a signature of this frame
- `linfo::MethodInstance`: The `MethodInstance` containing the execution context

This type is very similar to `Base.StackTraces.StackFrame`, but its execution context is
collected during abstract interpration, not collected from actual execution.
"""
@withmixedhash struct VirtualFrame
    file::Symbol
    line::Int
    sig::Signature
    linfo::MethodInstance
end

"""
    VirtualStackTrace

Represents a virtual stack trace in the form of a vector of `VirtualFrame`.
The vector holds `VirtualFrame`s in order of "from entry call site to error point", i.e.
the first element is the `VirtualFrame` of the entry call site, and the last element is that
contains the error.
"""
const VirtualStackTrace = Vector{VirtualFrame}

const INFERENCE_ERROR_REPORT_FIELD_DECLS = [
    :(vst::VirtualStackTrace),
    :(msg::String),
    :(sig::Signature),
]

# get location information at the given program counter (or a current counter if not specified)
function get_virtual_frame(state::StateAtPC)
    sig = get_sig(state)
    file, line = get_file_line(state)
    linfo = isa(state, MethodInstance) ? state : first(state).linfo
    return VirtualFrame(file, line, sig, linfo)
end
get_virtual_frame(sv::InferenceState) = get_virtual_frame((sv, get_currpc(sv)))
function get_virtual_frame(linfo::MethodInstance)
    sig = get_sig(linfo)
    file, line = get_file_line(linfo)
    return VirtualFrame(file, line, sig, linfo)
end

get_file_line(s::StateAtPC)          = get_file_line(get_lin(s))
get_file_line(lin::LineInfoNode)     = lin.file, lin.line
get_file_line(linfo::MethodInstance) = begin
    def = linfo.def

    isa(def, Method) && return def.file, Int(def.line)

    # top-level
    src = linfo.uninferred::CodeInfo
    return get_file_line(first(src.linetable::LineTable)::LineInfoNode)
end

# signature
# ---------

# adapted from https://github.com/JuliaLang/julia/blob/0f11a7bb07d2d0d8413da05dadd47441705bf0dd/base/show.jl#L989-L1011
function get_sig(l::MethodInstance)
    def = l.def
    ret = if isa(def, Method)
        if isdefined(def, :generator) && l === def.generator
            # print(io, "MethodInstance generator for ")
            # show(io, def)
            sprint(show, def)
        else
            # print(io, "MethodInstance for ")
            sprint(show_tuple_as_call, def.name, l.specTypes)
        end
    else
        # print(io, "Toplevel MethodInstance thunk")
        # # `thunk` is not very much information to go on. If this
        # # MethodInstance is part of a stacktrace, it gets location info
        # # added by other means.  But if it isn't, then we should try
        # # to print a little more identifying information.
        # if !from_stackframe
        #     linetable = l.uninferred.linetable
        #     line = isempty(linetable) ? "unknown" : (lt = linetable[1]; string(lt.file) * ':' * string(lt.line))
        #     print(io, " from ", def, " starting at ", line)
        # end
        "toplevel"
    end
    return Signature(Any[ret])
end

@inline function show_tuple_as_call(io::IO, name::Symbol, @nospecialize(sig::Type))
    @static if hasmethod(Base.show_tuple_as_call, (IO, Symbol, Type), (:demangle, :kwargs, :argnames, :qualified))
        Base.show_tuple_as_call(io, name, sig; qualified = true)
    else
        Base.show_tuple_as_call(io, name, sig, false, nothing, nothing, true)
    end
end

@inline get_sig(s::StateAtPC) = get_sig(s, get_stmt(s))
@inline get_sig(s::StateAtPC, @nospecialize(x)) = Signature(_get_sig(s, x))
@inline _get_sig(@nospecialize args...) = first(_get_sig_type(args...))

# to help inference
_get_sig_type(@nospecialize args...) = __get_sig_type(args...)::Tuple{Vector{Any},Any}

function _get_callsig(s::StateAtPC, @nospecialize(f), args::Vector{Any};
                      splat::Bool = false)
    sig = _get_sig(s, f)
    push!(sig, '(')

    nargs = length(args)
    for (i, arg) in enumerate(args)
        arg_sig = _get_sig(s, arg)
        append!(sig, arg_sig)
        if i ≠ nargs
            push!(sig, ", ")
        else
            splat && push!(sig, "...")
        end
    end
    push!(sig, ')')

    return sig
end

function __get_sig_type(s::StateAtPC, expr::Expr)
    head = expr.head
    if head === :call
        f = first(expr.args)
        args = expr.args[2:end]

        # special case splat call signature
        if isa(f, GlobalRef) && f.name === :_apply_iterate && begin
                itf = first(args)
                isa(itf, GlobalRef) && itf.name === :iterate
            end
            f = args[2]
            args = args[3:end]
            return _get_callsig(s, f, args; splat = true), nothing
        else
            return _get_callsig(s, f, args), nothing
        end
    elseif head === :invoke
        f = expr.args[2]
        args = expr.args[3:end]
        return _get_callsig(s, f, args), nothing
    elseif head === :(=)
        sigtyp = _get_sig_type(s, last(expr.args))
        sv = first(s)
        if isa(sv, InferenceState)
            lhs = first(expr.args)
            if isa(lhs, SlotNumber)
                name = sv.src.slotnames[slot_id(lhs)]
                sig, typ = sigtyp
                pushfirst!(sig, string(name), " = ")
                return sig, typ
            end
        end
        return sigtyp
    elseif head === :static_parameter
        typ = widenconst(first(s).sptypes[first(expr.args)::Int])
        return Any['_', typ], typ
    else
        return Any[string(expr)], nothing
    end
end
function __get_sig_type((sv, _)::StateAtPC, ssa::SSAValue)
    news = (sv, ssa.id)
    if isa(sv, OptimizationState)
        # when working on `OptimizationState`, the SSA traverse could be really long because
        # of inlining, so just give up for such a case
        typ = widenconst(ignorelimited(ignorenotfound(get_ssavaluetype(news))))
        sig = Any["%$(ssa.id)", typ]
    else
        # XXX the same problem _may_ happen for `InferenceState` too ?
        sig, sig_typ = _get_sig_type(news, get_stmt(news))
        typ = widenconst(ignorelimited(ignorenotfound(get_ssavaluetype(news))))
        sig_typ == typ || push!(sig, typ) # XXX I forgot why I added this line ...
    end
    return sig, typ
end
function __get_sig_type(s::StateAtPC, slot::SlotNumber)
    sv = first(s)
    name = get_slotname(sv, slot)
    sig = string(name)
    if isempty(sig)
        sig = string(slot) # fallback if no explicit slotname
    end
    if istoplevel(sv)
        # this is a abstract global variable, form the global reference
        return _get_sig_type(s, GlobalRef(sv.linfo.def::Module, name))
    else
        # we can use per-program counter type after inference
        t = (isa(sv, InferenceState) && sv.inferred) ? get_slottype(sv, slot) : get_slottype(s, slot)
        typ = widenconst(ignorelimited(t))
        return Any[sig, typ], typ
    end
end
# NOTE `Argument` is introduced by optimization, and so we don't need to handle abstract global variable here, etc.
function __get_sig_type((sv, _)::StateAtPC, arg::Argument)
    name = get_slotname(sv, arg.n)
    sig = string(name)
    typ = widenconst(ignorelimited(get_slottype(sv, arg))) # after optimization we shouldn't use `get_slottype(::StateAtPC, ::Any)`
    return Any[sig, typ], typ
end
__get_sig_type(_::StateAtPC, gr::GlobalRef) = Any[string(gr.mod, '.', gr.name)], nothing
__get_sig_type(_::StateAtPC, name::Symbol) = Any[repr(name; context = :compact => true)], nothing
function __get_sig_type(s::StateAtPC, gotoifnot::GotoIfNot)
    sig  = Any[string("goto %", gotoifnot.dest, " if not "), _get_sig(s, gotoifnot.cond)...]
    return sig, nothing
end
function __get_sig_type(s::StateAtPC, rn::ReturnNode)
    sig = is_unreachable(rn) ? Any["unreachable"] : Any["return ", _get_sig(s, rn.val)...]
    return sig, nothing
end
function __get_sig_type(::StateAtPC, qn::QuoteNode)
    typ = typeof(qn.value)
    return Any[string(qn), typ], typ
end
__get_sig_type(::StateAtPC, @nospecialize(x)) = Any[repr(x; context = :compact => true)], nothing

# new report
# ----------

"""
    InferenceErrorReport

An interface type of error reports that JET collects by abstract interpration.
If `T` implements this interface, the following requirements should be satisfied:

---
- **Required fields** \\
  `T` should have the following fields, which explains _where_ and _why_ this error is reported:
  * `vst::VirtualStackTrace`: a virtual stack trace of the error
  * `msg::String`: explains why this error is reported
  * [`sig::Signature`](@ref Signature): a signature of the error point

  Note that `T` can still have additional, specific fields.
---
- **A constructor interface to create `T` from abstraction interpretation** \\
  `T<:InferenceErrorReport` has the default constructor

      T(::AbstractAnalyzer, sv::InferenceState, spec_args...)

  which works when `T` is reported when `sv`'s program counter (`sv.currpc`) points to that
  of statement where the error may happen. If so `T` just needs to overload

      JET.get_msg(::Type{T}, ::InferenceState, spec_args...) -> msg::String

  to provide the message that describes why this error is reported (otherwise the senseless
  default message will be used).

  ---

  If `T` is reported when `sv`'s program counter (`sv.currpc`) may not point to the error
  location or even `sv::InferenceState` isn't available, `T` can implement its own constructor method.
---
- **A contructor interface to create `T` from the global report cache** \\
  In order to be cached and restored from [`JET_CACHE`](@ref), `T` _**must**_ implement
  the following interfaces:
  * `JET.get_spec_args(::T) -> Tuple{...}`:
    returns fields that are specific to `T`, which is internally used by the caching logic
  * `T(vst::VirtualStackTrace, msg::String, sig::Signature spec_args::Tuple{...}) -> T`:
    constructor to create `T` from the cache, which should expand `spec_args` into each specific field
---

To satisfy these requirements manually will be very tedious.
JET internally uses `@reportdef` utility macro, which takes the `struct` definition of
`InferenceErrorReport` and automatically defines the `struct` itself and the cache interfaces.

See also: [`VirtualStackTrace`](@ref), [`VirtualFrame`](@ref)
"""
abstract type InferenceErrorReport end

# to help inference
function Base.getproperty(er::InferenceErrorReport, sym::Symbol)
    return if sym === :vst
        getfield(er, sym)::VirtualStackTrace
    elseif sym === :msg
        getfield(er, sym)::String
    elseif sym === :sig
        getfield(er, sym)::Signature
    else
        getfield(er, sym) # fallback
    end
end

function Base.show(io::IO, report::InferenceErrorReport)
    print(io, typeof(report).name.name, '(')
    for a in report.sig
        _print_signature(io, a, (; annotate_types = true); bold = true)
    end
    print(io, ')')
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", report::InferenceErrorReport) =
    return report

# the default constructor to create a report from abstract interpretation
function (T::Type{<:InferenceErrorReport})(state, @nospecialize(spec_args...))
    vf = get_virtual_frame(state)
    msg = get_msg(T, state, spec_args...)
    return T([vf], msg, vf.sig, spec_args...)
end

get_msg(T::Type{<:InferenceErrorReport}, @nospecialize(_...)) = error("`get_msg` is not implemented for $T")
get_spec_args(T::Type{<:InferenceErrorReport}) =                error("`get_spec_args` is not implemented for $T")

# cache
# -----

@eval struct InferenceErrorReportCache
    T::Type{<:InferenceErrorReport}
    $(INFERENCE_ERROR_REPORT_FIELD_DECLS...)
    spec_args::NTuple{N,Any} where N
end

function cache_report!(cache, @nospecialize(report::InferenceErrorReport))
    vst = copy(report.vst)
    new = InferenceErrorReportCache(typeof(report), vst, report.msg, report.sig, get_spec_args(report))
    return push!(cache, new)
end

restore_cached_report(cache::InferenceErrorReportCache) =
    return cache.T(copy(cache.vst), cache.msg, cache.sig, cache.spec_args)::InferenceErrorReport

# utility
# -------

# TODO parametric definition

# a simple utility macro to define `InferenceErrorReport` w/o code duplications
macro reportdef(ex)
    @assert @capture(ex, struct T_ <: S_; spec_sigs__; end)
    @assert Core.eval(__module__, S) <: InferenceErrorReport

    spec_decls = map(spec_sigs) do x
        if @isexpr(x, :macrocall) && x.args[1] === Symbol("@nospecialize")
            return x.args[3]
        end
        return x
    end
    spec_names = extract_decl_name.(spec_decls)
    spec_types = esc.(extract_decl_type.(spec_decls))

    T, S = esc(T), esc(S)

    # cache constructor
    cache_constructor_sig = :($T(vst::VirtualStackTrace,
                                 msg::String,
                                 sig::Signature,
                                 @nospecialize(spec_args::Tuple),
                                 ))
    cache_constructor_call = :($T(vst, msg, sig))
    for (i, spec_type) in enumerate(spec_types)
        push!(cache_constructor_call.args, :($(esc(:spec_args))[$i]::$spec_type)) # needs escape since `@nospecialize`d
    end
    cache_constructor = Expr(:function, cache_constructor_sig, Expr(:block, __source__,
        :(return @inbounds $cache_constructor_call),
    ))

    # cache helper
    spec_getter_sig = :($(GlobalRef(JET, :get_spec_args))(report::$T))
    spec_getter_tuple = Expr(:tuple)
    for spec_name in spec_names
        getter = Expr(:call, GlobalRef(Base, :getproperty), :report, QuoteNode(spec_name))
        push!(spec_getter_tuple.args, getter)
    end
    spec_getter = Expr(:function, spec_getter_sig, Expr(:block, __source__,
        :(return $spec_getter_tuple::Tuple{$(spec_types...)}),
    ))

    return quote
        Base.@__doc__ struct $T <: $S
            $(INFERENCE_ERROR_REPORT_FIELD_DECLS...)
            $(map(esc, spec_decls)...)
            # esc is needed here since signanture might be `@nospecialize`d
            function $T($(INFERENCE_ERROR_REPORT_FIELD_DECLS...), $(map(esc, spec_sigs)...))
                new($(extract_decl_name.(INFERENCE_ERROR_REPORT_FIELD_DECLS)...), $(map(esc, spec_names)...))
            end
        end

        $cache_constructor

        $spec_getter
    end
end

extract_decl_name(@nospecialize(x)) = (@isexpr(x, :(::)) ? first(x.args) : x)::Symbol
extract_decl_type(@nospecialize(x)) = @isexpr(x, :(::)) ? last(x.args) : GlobalRef(Core, :Any)


text:
# Deep Learning with Flux: A 60 Minute Blitz
# =====================

# This is a quick intro to [Flux](https://github.com/FluxML/Flux.jl) loosely
# based on [PyTorch's
# tutorial](https://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html).
# It introduces basic Julia programming, as well Zygote, a source-to-source
# automatic differentiation (AD) framework in Julia.
# We'll use these tools to build a very simple neural network.

# Arrays
# -------

# The starting point for all of our models is the `Array` (sometimes referred to
# as a `Tensor` in other frameworks). This is really just a list of numbers,
# which might be arranged into a shape like a square. Let's write down an array
# with three elements.

x = [1, 2, 3]

# Here's a matrix – a square array with four elements.

x = [1 2; 3 4]

# We often work with arrays of thousands of elements, and don't usually write
# them down by hand. Here's how we can create an array of 5×3 = 15 elements,
# each a random number from zero to one.

x = rand(5, 3)

# There's a few functions like this; try replacing `rand` with `ones`, `zeros`,
# or `randn` to see what they do.

# By default, Julia works stores numbers is a high-precision format called
# `Float64`. In ML we often don't need all those digits, and can ask Julia to
# work with `Float32` instead. We can even ask for more digits using `BigFloat`.

x = rand(BigFloat, 5, 3)
#-
x = rand(Float32, 5, 3)

# We can ask the array how many elements it has.

length(x)

# Or, more specifically, what size it has.

size(x)

# We sometimes want to see some elements of the array on their own.

x
#-
x[2, 3]

# This means get the second row and the third column. We can also get every row
# of the third column.

x[:, 3]

# We can add arrays, and subtract them, which adds or subtracts each element of
# the array.

x + x
#-
x - x

# Julia supports a feature called *broadcasting*, using the `.` syntax. This
# tiles small arrays (or single numbers) to fill bigger ones.

x .+ 1

# We can see Julia tile the column vector `1:5` across all rows of the larger
# array.

zeros(5,5) .+ (1:5)

# The x' syntax is used to transpose a column `1:5` into an equivalent row, and
# Julia will tile that across columns.

zeros(5,5) .+ (1:5)'

# We can use this to make a times table.

(1:5) .* (1:5)'

# Finally, and importantly for machine learning, we can conveniently do things like
# matrix multiply.

W = randn(5, 10)
x = rand(10)
W * x

# Julia's arrays are very powerful, and you can learn more about what they can
# do [here](https://docs.julialang.org/en/v1/manual/arrays/).

# ### CUDA Arrays

# CUDA functionality is provided separately by the [CUDA
# package](https://github.com/JuliaGPU/CUDA.jl). If you have a GPU and CUDA
# available, you can run `] add CUDA` in a REPL or IJulia to get it.

# Once CUDA is loaded you can move any array to the GPU with the `cu`
# function, and it supports all of the above operations with the same syntax.

## using CUDA
## x = cu(rand(5, 3))

# Automatic Differentiation
# -------------------------

# You probably learned to take derivatives in school. We start with a simple
# mathematical function like

f(x) = 3x^2 + 2x + 1

f(5)

# In simple cases it's pretty easy to work out the gradient by hand – here it's
# `6x+2`. But it's much easier to make Flux do the work for us!

using Flux: gradient

df(x) = gradient(f, x)[1]

df(5)

# You can try this with a few different inputs to make sure it's really the same
# as `6x+2`. We can even do this multiple times (but the second derivative is a
# fairly boring `6`).

ddf(x) = gradient(df, x)[1]

ddf(5)

# Flux's AD can handle any Julia code you throw at it, including loops,
# recursion and custom layers, so long as the mathematical functions you call
# are differentiable. For example, we can differentiate a Taylor approximation
# to the `sin` function.

mysin(x) = sum((-1)^k*x^(1+2k)/factorial(1+2k) for k in 0:5)

x = 0.5

mysin(x), gradient(mysin, x)
#-
sin(x), cos(x)

# You can see that the derivative we calculated is very close to `cos(x)`, as we
# expect.

# This gets more interesting when we consider functions that take *arrays* as
# inputs, rather than just a single number. For example, here's a function that
# takes a matrix and two vectors (the definition itself is arbitrary)

myloss(W, b, x) = sum(W * x .+ b)

W = randn(3, 5)
b = zeros(3)
x = rand(5)

gradient(myloss, W, b, x)

# Now we get gradients for each of the inputs `W`, `b` and `x`, which will come
# in handy when we want to train models.

# Because ML models can contain hundreds of parameters, Flux provides a slightly
# different way of writing `gradient`. We instead mark arrays with `param` to
# indicate that we want their derivatives. `W` and `b` represent the weight and
# bias respectively.

using Flux: params

W = randn(3, 5)
b = zeros(3)
x = rand(5)

y(x) = sum(W * x .+ b)

grads = gradient(()->y(x), params([W, b]))

grads[W], grads[b]

# We can now grab the gradients of `W` and `b` directly from those parameters.

# This comes in handy when working with *layers*. A layer is just a handy
# container for some parameters. For example, `Dense` does a linear transform
# for you.

using Flux

m = Dense(10, 5)

# We can easily get the parameters of any layer or model with params with
# `params`.

params(m)

# This makes it very easy to calculate the gradient for all
# parameters in a network, even if it has many parameters.
x = rand(Float32, 10)
m = Chain(Dense(10, 5, relu), Dense(5, 2), softmax)
l(x) = Flux.Losses.crossentropy(m(x), [0.5, 0.5])
grads = gradient(params(m)) do
    l(x)
end
for p in params(m)
    println(grads[p])
end

# You don't have to use layers, but they can be convenient for many simple kinds
# of models and fast iteration.

# The next step is to update our weights and perform optimisation. As you might be
# familiar, *Gradient Descent* is a simple algorithm that takes the weights and steps
# using a learning rate and the gradients. `weights = weights - learning_rate * gradient`
# (note that `Flux.Optimise.update!(x, x̄)` already updates with the negative of x̄`).
using Flux.Optimise: update!, Descent
η = 0.1
for p in params(m)
  update!(p, η * grads[p])
end

# While this is a valid way of updating our weights, it can get more complicated as the
# algorithms we use get more involved.

# Flux comes with a bunch of pre-defined optimisers and makes writing our own really simple.
# We just give it the learning rate η

opt = Descent(0.01)

# `Training` a network reduces down to iterating on a dataset multiple times, performing these
# steps in order. Just for a quick implementation, let’s train a network that learns to predict
# `0.5` for every input of 10 floats. `Flux` defines the `train!` function to do it for us.

data, labels = rand(10, 100), fill(0.5, 2, 100)
loss(x, y) = Flux.Losses.crossentropy(m(x), y)
Flux.train!(loss, params(m), [(data,labels)], opt)
# You don't have to use `train!`. In cases where arbitrary logic might be better suited,
# you could open up this training loop like so:

# ```julia
#   for d in training_set # assuming d looks like (data, labels)
#     # our super logic
#     gs = gradient(params(m)) do #m is our model
#       l = loss(d...)
#     end
#     update!(opt, params(m), gs)
#   end
# ```

# Training a Classifier
# ---------------------

# Getting a real classifier to work might help cement the workflow a bit more.
# [CIFAR10](url) is a dataset of 50k tiny training images split into 10 classes.

# We will do the following steps in order:

# * Load CIFAR10 training and test datasets
# * Define a Convolution Neural Network
# * Define a loss function
# * Train the network on the training data
# * Test the network on the test data

# Loading the Dataset

# [Metalhead.jl](https://github.com/FluxML/Metalhead.jl) is an excellent package
# that has a number of predefined and pretrained computer vision models.
# It also has a number of dataloaders that come in handy to load datasets.

using Statistics
using Flux, Flux.Optimise
using Images: channelview
using Metalhead
using Metalhead: trainimgs, valimgs
using Images.ImageCore
using Flux: onehotbatch, onecold, flatten
using Base.Iterators: partition
# using CUDA

# The image will give us an idea of what we are dealing with.
# ![title](https://pytorch.org/tutorials/_images/cifar10.png)

Metalhead.download(CIFAR10)
X = trainimgs(CIFAR10)
labels = onehotbatch([X[i].ground_truth.class for i in 1:50000],1:10)

# Let's take a look at a random image from the dataset

image(x) = x.img # handy for use later
ground_truth(x) = x.ground_truth
image.(X[rand(1:end, 10)])

# The images are simply 32 X 32 matrices of numbers in 3 channels (R,G,B). We can now
# arrange them in batches of say, 1000 and keep a validation set to track our progress.
# This process is called minibatch learning, which is a popular method of training
# large neural networks. Rather that sending the entire dataset at once, we break it
# down into smaller chunks (called minibatches) that are typically chosen at random,
# and train only on them. It is shown to help with escaping
# [saddle points](https://en.wikipedia.org/wiki/Saddle_point).

# Defining a `getarray` function would help in converting the matrices to `Float` type.

getarray(X) = float.(permutedims(channelview(X), (3, 2, 1)))
imgs = [getarray(X[i].img) for i in 1:50000]

# The first 49k images (in batches of 1000) will be our training set, and the rest is
# for validation. `partition` handily breaks down the set we give it in consecutive parts
# (1000 in this case). `cat` is a shorthand for concatenating multi-dimensional arrays along
# any dimension.

train = ([(cat(imgs[i]..., dims = 4), labels[:,i]) for i in partition(1:49000, 1000)]) |> gpu
valset = 49001:50000
valX = cat(imgs[valset]..., dims = 4) |> gpu
valY = labels[:, valset] |> gpu

# ## Defining the Classifier
# --------------------------
# Now we can define our Convolutional Neural Network (CNN).

# A convolutional neural network is one which defines a kernel and slides it across a matrix
# to create an intermediate representation to extract features from. It creates higher order
# features as it goes into deeper layers, making it suitable for images, where the structure of
# the subject is what will help us determine which class it belongs to.

m = Chain(
  Conv((5,5), 3=>16, relu),
  MaxPool((2,2)),
  Conv((5,5), 16=>8, relu),
  MaxPool((2,2)),
  flatten,
  Dense(200, 120),
  Dense(120, 84),
  Dense(84, 10),
  softmax) |> gpu

#-
# We will use a crossentropy loss and the Momentum optimiser here. Crossentropy will be a
# good option when it comes to working with multiple independent classes. Momentum smooths out
# the noisy gradients and helps towards a smooth convergence. Gradually lowering the
# learning rate along with momentum helps to maintain a bit of adaptivity in our optimisation,
# preventing us from overshooting our desired destination.
#-

using Flux: Momentum

loss(x, y) = Flux.Losses.crossentropy(m(x), y)
opt = Momentum(0.01)

# We can start writing our train loop where we will keep track of some basic accuracy
# numbers about our model. We can define an `accuracy` function for it like so.

accuracy(x, y) = mean(onecold(m(x), 1:10) .== onecold(y, 1:10))

# ## Training
# -----------

# Training is where we do a bunch of the interesting operations we defined earlier,
# and see what our net is capable of. We will loop over the dataset 10 times and
# feed the inputs to the neural network and optimise.

epochs = 10

for epoch = 1:epochs
  for d in train
    gs = gradient(params(m)) do
      l = loss(d...)
    end
    update!(opt, params(m), gs)
  end
  @show accuracy(valX, valY)
end

# Seeing our training routine unfold gives us an idea of how the network learnt the
# This is not bad for a small hand-written network, trained for a limited time.

# Training on a GPU
# -----------------

# The `gpu` functions you see sprinkled through this bit of the code tell Flux to move
# these entities to an available GPU, and subsequently train on it. No extra faffing
# about required! The same bit of code would work on any hardware with some small
# annotations like you saw here.

# ## Testing the Network
# ----------------------

# We have trained the network for 10 passes over the training dataset. But we need to
# check if the network has learnt anything at all.

# We will check this by predicting the class label that the neural network outputs, and
# checking it against the ground-truth. If the prediction is correct, we add the sample
# to the list of correct predictions. This will be done on a yet unseen section of data.

# Okay, first step. Let us perform the exact same preprocessing on this set, as we did
# on our training set.

valset = valimgs(CIFAR10)
valimg = [getarray(valset[i].img) for i in 1:10000]
labels = onehotbatch([valset[i].ground_truth.class for i in 1:10000],1:10)
test = gpu.([(cat(valimg[i]..., dims = 4), labels[:,i]) for i in partition(1:10000, 1000)])

# Next, display some of the images from the test set.

ids = rand(1:10000, 10)
image.(valset[ids])

# The outputs are energies for the 10 classes. Higher the energy for a class, the more the
# network thinks that the image is of the particular class. Every column corresponds to the
# output of one image, with the 10 floats in the column being the energies.

# Let's see how the model fared.

rand_test = getarray.(image.(valset[ids]))
rand_test = cat(rand_test..., dims = 4) |> gpu
rand_truth = ground_truth.(valset[ids])
m(rand_test)

# This looks similar to how we would expect the results to be. At this point, it's a good
# idea to see how our net actually performs on new data, that we have prepared.

accuracy(test[1]...)

# This is much better than random chance set at 10% (since we only have 10 classes), and
# not bad at all for a small hand written network like ours.

# Let's take a look at how the net performed on all the classes performed individually.

class_correct = zeros(10)
class_total = zeros(10)
for i in 1:10
  preds = m(test[i][1])
  lab = test[i][2]
  for j = 1:1000
    pred_class = findmax(preds[:, j])[2]
    actual_class = findmax(lab[:, j])[2]
    if pred_class == actual_class
      class_correct[pred_class] += 1
    end
    class_total[actual_class] += 1
  end
end

class_correct ./ class_total

# The spread seems pretty good, with certain classes performing significantly better than the others.
# Why should that be?

max_stars_repo_name: Ossifragus/model-zoo

max_issues_repo_name: Ossifragus/model-zoo

max_forks_repo_name: Ossifragus/model-zoo


text:
# Precompile Utilities

"""
    $SIGNATURES

Load and solve a canonical or user factor graph to warm up---precompile---several RoME/Caesar related functions.
"""
function warmUpSolverJIT(;fg::AbstractDFG=generateCanonicalFG_Hexagonal(),
                          drawtree::Bool=true )::Nothing
  #

  fcts = ls(fg, :x0)
  fcts = ls(fg)
  fcts = lsf(fg, :x0f1)
  fcts = lsf(fg)
  getSolverParams(fg).drawtree = drawtree
  tree = solveTree!(fg)
  nothing
end