text:
import PyPlot.plot

"""
    plot_field(field::Array{<:Real,2}, a1::Array{<:Real,1}, a2::Array{<:Real,1}; cmap="coolwarm", vmin=nothing, vmax=nothing)

Plot the real-valued `field` on a unit cell with lattice vectors `a1` and `a2`.
"""
function plot_field(field::Array{<:Real,2}, a1::Array{<:Real,1}, a2::Array{<:Real,1}; cmap="coolwarm", vmin=nothing, vmax=nothing, label=nothing)
    ps = range(-0.5, stop=0.5, length=size(field,1)+1)
    qs = range(-0.5, stop=0.5, length=size(field,2)+1)
    xs = [p*a1[1]+q*a2[1] for p in ps, q in qs]
    ys = [p*a1[2]+q*a2[2] for p in ps, q in qs]
    if cmap == "coolwarm" && vmin == nothing && vmax == nothing
        vlim = maximum(abs.(field))
        vmin = -vlim
        vmax = +vlim
    end
    if cmap == "viridis" && vmin == nothing && vmax == nothing
        vmin = 0
        vmax = maximum(abs.(field))
    end
    pcolormesh(xs, ys, field, cmap=cmap, vmin=vmin, vmax=vmax)
    gca().set_aspect("equal")
    axis("off")
    colorbar()
    if label != nothing
        title(label)
    end
end

"""
    plot_field(field::Array{<:Real,2}, a1::Array{<:Real,1}, a2::Array{<:Real,1}; cmap="coolwarm", vmin=nothing, vmax=nothing)

Plot the complex-valued `field` on a unit cell with lattice vectors `a1` and `a2`.
"""
function plot_field(field::Array{<:Complex,2}, a1::Array{<:Real,1}, a2::Array{<:Real,1}; cmap="coolwarm", vmin=nothing, vmax=nothing, label=nothing)
    # Plot real and imaginary parts with the same color scale
    if cmap == "coolwarm" && vmin == nothing && vmax == nothing
        vlim = maximum(abs.(field))
        vmin = -vlim
        vmax = +vlim
    end
    if cmap == "viridis" && vmin == nothing && vmax == nothing
        vmin = 0
        vmax = maximum(abs.(field))
    end
    figure(figsize=(4,2))
    subplot(1,2,1)
    plot_field(real(field), a1, a2, cmap=cmap, vmin=vmin, vmax=vmax, label="Re[$label]")
    subplot(1,2,2)
    plot_field(imag(field), a1, a2, cmap=cmap, vmin=vmin, vmax=vmax, label="Im[$label]")
end

"""
    get_field(data::AbstractVector{<:Complex}, basis::PlaneWaveBasis;
                        k0=[0,0], t1s=-0.5:0.01:0.5, t2s=-0.5:0.01:0.5)

Convert the `data` from a `PlaneWaveBasis` to a real space grid.
"""
function get_field(data::AbstractVector{<:Complex}, basis::PlaneWaveBasis;
                    k0=[0,0], t1s=-0.5:0.01:0.5, t2s=-0.5:0.01:0.5)
    # Generate real space coordinates
    b1, b2 = basis.b1, basis.b2
    a1, a2 = bs_to_as(b1, b2)
    xs = [t1*a1[1]+t2*a2[1] for t1 in t1s, t2 in t2s]
    ys = [t1*a1[2]+t2*a2[2] for t1 in t1s, t2 in t2s]
    out = zeros(ComplexF64, length(t1s), length(t2s))
    # Inverse Fourier transform
    kxs = k0[1] .+ basis.kxs
    kys = k0[2] .+ basis.kys
    for (kx,ky,u) in zip(kxs, kys, data)
        for c in CartesianIndices(out)
            phase = kx*xs[c] + ky*ys[c]
            out[c] += u * exp(1im*phase)
        end
    end
    return out
end

function get_field_FDFD(data::AbstractVector{<:Complex}, basis::FDFDBasis; k0 = [0, 0])
    dx, dy = basis.resolution[1], basis.resolution[2]
    a1, a2 = basis.a1, basis.a2
    Nx, Ny = Int(1/dx), Int(1/dy)
    out = zeros(ComplexF64, Nx, Ny)
    for m = 1:Nx
        for n = 1:Ny
            Lx = (m - 0.5) * dx
            Ly = (n - 0.5) * dy
            index = Int((n - 1)*Nx + m)
            out[m, n] = data[index] * exp(-1im * k0[1] * Lx) * exp(-1im * k0[2] *Ly) #u function
        end
    end
    return out
end

function plot_field(data::AbstractVector{<:Complex}, basis::FDFDBasis; k0 = [0, 0], cmap="coolwarm", vmin=nothing, vmax=nothing, label=nothing)
    field = get_field_FDFD(data, basis, k0 = k0)
    plot_field(field, basis.a1, basis.a2, cmap=cmap, vmin=vmin, vmax=vmax, label=label)
end

"""
    plot_field(data::AbstractVector{<:Complex}, basis::PlaneWaveBasis;
                        k0=[0,0], cmap="coolwarm", vmin=nothing, vmax=nothing, label=nothing)

Converts the plane-wave amplitudes to real space using [`get_field`](@ref)
and then plots the field.
"""
function plot_field(data::AbstractVector{<:Complex}, basis::PlaneWaveBasis; k0=[0,0], cmap="coolwarm", vmin=nothing, vmax=nothing, label=nothing)
    field = get_field(data, basis, k0=k0)
    a1, a2 = bs_to_as(basis.b1, basis.b2)
    plot_field(field, a1, a2, cmap=cmap, vmin=vmin, vmax=vmax, label=label)
end

"""
    plot(geometry::Geometry)

Plot the permittivity and permeability of the `geometry` in real space.
"""
function PyPlot.plot(geometry::Geometry)
    # Plot permittivity
    plot_field(geometry.ep.epsxx, geometry.a1, geometry.a2, cmap="viridis", label="ϵ")
    # Plot permeability
    plot_field(geometry.mu.muxx, geometry.a1, geometry.a2, cmap="viridis", label="μ")
end

"""
    plot(solver::Solver)

Plot the representation of the `solver`'s geometry in real space.

The `solver` approximates the geometry using a truncated basis of plane waves,
so `plot(solver)` lets us judge how accurately the geometry is represented.
"""
function PyPlot.plot(solver::Solver)
    col = findfirst((solver.basis.ps .==0) .& (solver.basis.qs .== 0))
    # Plot permittivity
    plot_field(solver.epc[:,col], solver.basis, cmap="viridis", label="ϵ")
    # Plot permeability
    plot_field(solver.muc[:,col], solver.basis, cmap="viridis", label="μ")
end

"""
    plot(mode::Mode, [bloch_phase=true])

Plot the mode in real space.

To plot only the cell-periodic part of the Bloch wave, set `bloch_phase=false`.
"""
function PyPlot.plot(mode::Mode; bloch_phase=true)
    if bloch_phase
        label = mode.label
        plot_field(mode.data, mode.basis, label=label, k0=mode.k0)
    else
        label = L"e^{-i k \cdot r} \cdot" * mode.label
        plot_field(mode.data, mode.basis, label=label)
    end
end

function PyPlot.plot(mode::Mode_FDFD; bloch_phase=true)
    if bloch_phase
        label = mode.label
        plot_field(mode.data, mode.basis, label=label)
    else
        label = L"e^{-i k \cdot r} \cdot" * mode.label
        plot_field(mode.data, mode.basis, label=label , k0=mode.k0)
    end
end


text:
"""
It stores set of functions that need to be composed in order to prepare GLFW window and
display verticies needed for texture  display
"""
module PrepareWindowHelpers

using  ..ModernGlUtil, GLFW,ModernGL

export createDAtaBuffer
export createElementBuffer
export createVertexBuffer
export encodeDataFromDataBuffer
export controllWindowInput
export initializeWindow

"""
data is loaded into a buffer which passes it into thw GPU for futher processing
    - here the data is just passing the positions of verticies
    GL_STREAM_DRAW the data is set only once and used by the GPU at most a few times.
    GL_STATIC_DRAW the data is set only once and used many times.
    GL_DYNAMIC_DRAW the data is changed a lot and used many times.
    """
function createDAtaBuffer(positions)
    vbo = Ref(GLuint(0))   # initial value is irrelevant, just allocate space
    glGenBuffers(1, vbo)
    glBindBuffer(GL_ARRAY_BUFFER, vbo[])
    glBufferData(GL_ARRAY_BUFFER, sizeof(positions), positions, GL_STATIC_DRAW)
return vbo
end #createDAtaBuffer

"""
Similar to the VBO we bind the EBO and copy the indices into the buffer with glBufferData.
    """
function createElementBuffer(elements)
    ebo = Ref(GLuint(0))
    glGenBuffers(1, ebo)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ebo[])
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(elements), elements, GL_STATIC_DRAW)
return ebo
end #createElementBuffer

"""
vertex buffer keeping things simpler
    """
function createVertexBuffer()
    vao = Ref(GLuint(0))
    glGenVertexArrays(1, vao)
    glBindVertexArray(vao[])
    return vao
end #createVertexBuffer

glVertexAttribSettingStr= """
showing how openGL should read data from buffer in GPU
in case of code like below it would mean:

first parameter specifies which vertex attribute we want to configure Remember that we specified the location of the position vertex attribute in the vertex shader
next argument specifies the size of the vertex attribute. The vertex attribute is a vec2 so it is composed of 2 values.
The third argument specifies the type of the data which is GL_FLOAT
next argument specifies if we want the data to be normalized. If we’re inputting integer data types like int, byte and we’ve set this to GL_TRUE
    The fifth argument is known as the stride and tells us the space between consecutive vertex attributes. Since the next set of position data is
     located exactly 2 times the size of a float we could’ve also specified the stride as 0 to let OpenGL determine the stride
    he last parameter is of type void* and thus requires that weird cast. This is the offset of where the position data begins in the buffer.

glVertexAttribPointer positionAttribute, 2, GL_FLOAT, false, 0, C_NULL

    The position data is stored as 32-bit  so 4 byte floating point values.
    Each position is composed of 2 of those values.
    """
    @doc glVertexAttribSettingStr
function glVertexAttribSetting(positionAttribute)
    glVertexAttribPointer(positionAttribute, 2, GL_FLOAT, false, 0, C_NULL)

end #glVertexAttribSetting

"""
how data should be read from data buffer
    """
function encodeDataFromDataBuffer()
    typee = Float32

    # position attribute
    glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 8 * sizeof(typee), C_NULL);
    glEnableVertexAttribArray(0);
    # color attribute
    glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, 8 * sizeof(typee),  Ptr{Nothing}(3 * sizeof(typee)));
    glEnableVertexAttribArray(1);
    # texture coord attribute
    glVertexAttribPointer(2, 2, GL_FLOAT, GL_FALSE, 8 * sizeof(typee),  Ptr{Nothing}(6 * sizeof(typee)));
    glEnableVertexAttribArray(2);

end #encodeDataFromDataBuffer

"""
it will generally be invoked on GLFW.PollEvents()  in event loop and now depending on
what will be pressed or clicked it will lead to diffrent actions
"""
function controllWindowInput(window)

    # Input callbacks
GLFW.SetKeyCallback(window, (_, key, scancode, action, mods) -> begin
	name = GLFW.GetKeyName(key, scancode)
	if name == nothing
		println("scancode $scancode  $(typeof(scancode))", action)
        println("action $action $(typeof(action))")

	else
		println("key $name $(typeof(name))", action)
		println("action $action $(typeof(action))")
	end
end)
end #controllWindowInputDoc

"""
modified from ModernGL.jl github page  and GLFW page
stores primary
    """
function initializeWindow(windowWidth::Int,windowHeight::Int)
	GLFW.Init()

	# Create a windowed mode window and its OpenGL context
	window = GLFW.CreateWindow(windowWidth, windowHeight, "Segmentation Visualization")
	# Make the window's context current
	GLFW.MakeContextCurrent(window)
	GLFW.ShowWindow(window)
	GLFW.SetWindowSize(window, windowWidth, windowHeight) # Seems to be necessary to guarantee that window > 0
	glViewport(0, 0, windowWidth, windowHeight)
	glDisable(GL_LIGHTING);
	glEnable(GL_TEXTURE_2D);
	println(createcontextinfo())
	return window
	end #initializeWindow

end # PreperWindowHelpers

ext: jl


text:
#========== a mutable, typeless, almost-namedtuple ==========#

struct DotDict
    store::Dict{Symbol,Any}
end
DotDict(;kw...) = DotDict(Dict(pairs(kw)...))

Base.parent(x::DotDict) = getfield(x, :store)

Base.propertynames(x::DotDict) = Tuple(sort(collect(keys(parent(x)))))
Base.getproperty(x::DotDict, s::Symbol) = getindex(parent(x), s)
function Base.setproperty!(x::DotDict, s::Symbol, v)
    s in propertynames(x) || throw("DotDict has no field $s")
    T = typeof(getproperty(x, s))
    if T == Nothing
        setindex!(parent(x), v, s)
    else
        setindex!(parent(x), convert(T, v), s)
    end
end

function Base.show(io::IO, x::DotDict)
    print(io, "DotDict(")
    strs = map(k -> string(k, " = ", getproperty(x, k)), propertynames(x))
    print(io, join(strs, ", "), ")")
end

verboseprint(store) = begin
    printstyled("┌ store.\n", color=:blue)
    foreach(propertynames(store)) do k
        printstyled("│   $k = ", color=:blue)
        r = getproperty(store, k)
        if k ∉ [:outpre, :outex]
            println(repr(r))
        else
            str = repr(verbosetidy.(r))
            println(first(str, 150), length(str)>150 ? " ..." : "")
        end
    end
    printstyled("└\n", color=:blue)
    if store.verbose == 3
        if haskey(store, :outpre)
            printstyled("store.outpre = \n", color=:blue)
            printstyled(verbosetidy(store.outpre) , "\n", color=:green)
        end
        if haskey(store, :outex)
            printstyled("\nstore.outex = \n", color=:blue)
            printstyled(verbosetidy(store.outpre) , "\n", color=:green)
        end
    end
end

#========== capture macro ==========#
# My faster, more limited, version:

"""
    @capture_(ex, A_[ijk__])

Faster drop-in replacement for `MacroTools.@capture`, for a few patterns only:
* `A_[ijk__]` and `A_{ijk__}`
* `[ijk__]`
* `A_.field_`
* `A_ := B_` and  `A_ = B_` and `A_ += B_` etc.
* `f_(x_)`
"""
macro capture_(ex, pat::Expr)

    H = QuoteNode(pat.head)

    A,B = if pat.head in [:ref, :curly] && length(pat.args)==2 &&
        _endswithone(pat.args[1]) && _endswithtwo(pat.args[2]) # :( A_[ijk__] )
        _symbolone(pat.args[1]), _symboltwo(pat.args[2])

    elseif pat.head == :. && pat.args[2] isa QuoteNode &&
        _endswithone(pat.args[1]) && _endswithone(pat.args[2].value) # :( A_.field_ )
        _symbolone(pat.args[1]), _symbolone(pat.args[2].value)

    elseif pat.head == :call  && length(pat.args)==2 &&
        _endswithone(pat.args[1]) && _endswithone(pat.args[2]) # :( f_(x_) )
        _symbolone(pat.args[1]), _symbolone(pat.args[2])

    elseif pat.head in [:call, :(=), :(:=), :+=, :-=, :*=, :/=, :^=] &&
        _endswithone(pat.args[1]) && _endswithone(pat.args[2]) # :( A_ += B_ )
        _symbolone(pat.args[1]), _symbolone(pat.args[2])

    # elseif pat.head == :call  && length(pat.args)==3 && pat.args[1] == :!= &&
    #     _endswithone(pat.args[2]) && _endswithone(pat.args[3]) # :( A_ != B_ )
    #     H = QuoteNode(pat.args[1])
    #     _symbolone(pat.args[2]), _symbolone(pat.args[3])

    elseif pat.head == :vect && _endswithtwo(pat.args[1]) # :( [ijk__] )
        _symboltwo(pat.args[1]), gensym(:ignore)

    else
        throw("@capture_ doesn't work on pattern $pat")
    end

    @gensym res
    quote
        $A, $B = nothing, nothing
        $res = $_trymatch($ex, Val($H))
        if $res === nothing
            false
        else
            $A, $B = $res
            true
        end
    end |> esc
end

_endswithone(ex) = endswith(string(ex), '_') && !_endswithtwo(ex)
_endswithtwo(ex) = endswith(string(ex), "__")

_symbolone(ex) = Symbol(string(ex)[1:end-1])
_symboltwo(ex) = Symbol(string(ex)[1:end-2])

_getvalue(::Val{val}) where {val} = val

_trymatch(s, v) = nothing # Symbol, or other Expr
_trymatch(ex::Expr, pat::Union{Val{:ref}, Val{:curly}}) = # A_[ijk__] or A_{ijk__}
    if ex.head === _getvalue(pat)
        ex.args[1], ex.args[2:end]
    else
        nothing
    end
_trymatch(ex::Expr, ::Val{:.}) = # A_.field_
    if ex.head === :. && ex.args[2] isa QuoteNode
        ex.args[1], ex.args[2].value
    else
        nothing
    end
_trymatch(ex::Expr, pat::Val{:call}) =
    if ex.head === _getvalue(pat) && length(ex.args) == 2
        ex.args[1], ex.args[2]
    else
        nothing
    end
_trymatch(ex::Expr, pat::Union{Val{:(=)}, Val{:(:=)}, Val{:(+=)}, Val{:(-=)}, Val{:(*=)}, Val{:(/=)}, Val{:(^=)}}) =
    if ex.head === _getvalue(pat)
        ex.args[1], ex.args[2]
    else
        nothing
    end
# _trymatch(ex::Expr, pat::Val{:!=}) =
#     if ex.head === :call && length(ex.args) == 3 && ex.args[1] == :!=
#         ex.args[2], ex.args[3]
#     else
#         nothing
#     end
_trymatch(ex::Expr, ::Val{:vect}) = # [ijk__]
    if ex.head === :vect
        ex.args, nothing
    else
        nothing
    end

# Cases for Tullio:
# @capture(ex, B_[inds__].field_) --> @capture_(ex, Binds_.field_) && @capture_(Binds, B_[inds__])

#========== postwalk ==========#
# Copied verbatim from here:
# https://github.com/MikeInnes/MacroTools.jl/blob/master/src/utils.jl

walk(x, inner, outer) = outer(x)
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))

"""
    postwalk(f, expr)
Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *after* they have been transformed by the walk. See also
`prewalk`.
"""
postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

"""
    prewalk(f, expr)
Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *before* they have been transformed by the walk, and the
walk will be applied to whatever `f` returns.
This makes `prewalk` somewhat prone to infinite loops; you probably want to try
`postwalk` first.
"""
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)

replace(ex, s, s′) = prewalk(x -> x == s ? s′ : x, ex)

const MacroTools_prewalk = prewalk
const MacroTools_postwalk = postwalk

#========== prettify ==========#

using Base.Meta: isexpr

verbosetidy(expr) = MacroTools_postwalk(expr) do ex
        if isexpr(ex, :block)
            args = filter(x -> !(x isa LineNumberNode || x == nothing), ex.args)
            if length(args) == 1 && Meta.isexpr(args[1], :block)
                # disallow block(block(stuff))
                args[1]
            else
                Expr(ex.head, args...)
            end
        elseif isexpr(ex, :macrocall) && length(ex.args) >= 2
            # line number after macro name can't be dropped, but can be nothing:
            Expr(ex.head, ex.args[1], nothing, filter(x -> !(x isa LineNumberNode), ex.args[3:end])...)
        else
            ex
        end
    end

#========== the end ==========#


text:
export @gencxxf

using Cxx

const JlToCxxVarNameMap = Dict{Any,Any}(
    :AbstractCvMat => :(CVCore.handle),
    :Mat => :(CVCore.handle),
    :UMat => :(CVCore.handle),
)

function gencxxf_impl(callexpr::Expr, cxxname::String, varmap::Dict)
    while callexpr.head == :block
        if callexpr.args[1].head == :line
            callexpr = callexpr.args[2]
        else
            callexpr = callexpr.args[1]
        end
    end
    @assert callexpr.head == :call
    name = first(callexpr.args)
    args = callexpr.args[2:end]

    # @show name
    # @show args

    if args[1].head == :tuple
        args = args[1].args
    end

    cxxargs = []
    for arg in args
        if isa(arg, Symbol)
            push!(cxxargs, arg)
        elseif arg.head != :(::)
            push!(cxxargs, arg.args[1])
        else
            cxxarg = copy(arg)
            typname = arg.args[2]
            varname = arg.args[1]
            # Union
            if !isa(typname, Symbol) && typname.head == :curly && typname.args[1] == :Union
                if all(map(x -> haskey(varmap, x), typname.args[2:end]))
                    cxxarg.args[1] = Expr(:call, varmap[typname.args[2]], varname)
                else
                    error("Inconsistent Union type")
                end
            end
            if haskey(varmap, typname)
                cxxarg.args[1] = Expr(:call, varmap[typname], varname)
            end
            push!(cxxargs, cxxarg.args[1])
        end
    end

    icxxstr = """
    $(cxxname)($(join(map(x -> string("\$(", x, ")"), cxxargs), ", ")));
    """
    # @show icxxstr
    icxxcall = Expr(:macrocall, Symbol("@icxx_str"), icxxstr)

    f = esc(quote
            function $name($(args...))
                $icxxcall
            end
        end)
    # @show f
    f
end

macro gencxxf(callexpr, cxxname)
    gencxxf_impl(callexpr, cxxname, JlToCxxVarNameMap)
end


text:
"""
    flux_balance_analysis_vec(args...)::Union{Vector{Float64},Nothing}

A variant of FBA that returns a vector of fluxes in the same order as reactions
of the model, if the solution is found.

Arguments are passed to [`flux_balance_analysis`](@ref).
"""
function flux_balance_analysis_vec(args...; kwargs...)::Union{Vector{Float64},Nothing}
    optmodel = flux_balance_analysis(args...; kwargs...)
    is_solved(optmodel) || return nothing
    value.(optmodel[:x])
end

"""
    flux_balance_analysis_dict(model::MetabolicModel, args...)::Union{Dict{String, Float64},Nothing}

A variant of FBA that returns a dictionary assigning fluxes to reactions, if
the solution is found. Arguments are passed to [`flux_balance_analysis`](@ref).
"""
function flux_balance_analysis_dict(
    model::MetabolicModel,
    args...;
    kwargs...,
)::Union{Dict{String,Float64},Nothing}
    v = flux_balance_analysis_vec(model, args...; kwargs...)
    isnothing(v) && return nothing
    Dict(zip(reactions(model), v))
end

"""
    flux_balance_analysis(
        model::M,
        optimizer;
        modifications = [],
    ) where {M<:MetabolicModel}

Run flux balance analysis (FBA) on the `model` optionally specifying
`modifications` to the problem.  Basically, FBA solves this optimization problem:
```
max cᵀx
s.t. S x = b
     xₗ ≤ x ≤ xᵤ
```
See "Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat
Biotechnol 28, 245–248 (2010). https://doi.org/10.1038/nbt.1614" for more
information.

The `optimizer` must be set to a `JuMP`-compatible optimizer, such as
`GLPK.Optimizer` or `Tulip.Optimizer`

Optionally, you may specify one or more modifications to be applied to the
model before the analysis, such as
[`change_optimizer_attribute`](@ref),[`change_objective`](@ref), and
[`change_sense`](@ref).

Returns an optimized `JuMP` model.

# Example
```
model = load_model(StandardModel, "e_coli_core.json")
biomass = findfirst(model.reactions, "BIOMASS_Ecoli_core_w_GAM")
solution = flux_balance_analysis(model, GLPK.optimizer; modifications=[change_objective(biomass)])
```

"""
function flux_balance_analysis(
    model::M,
    optimizer;
    modifications = [],
) where {M<:MetabolicModel}

    opt_model = make_optimization_model(model, optimizer)

    for mod in modifications
        mod(model, opt_model)
    end

    COBREXA.JuMP.optimize!(opt_model)
    return opt_model
end


text:
using Mimi
using Distributions
using Query
using DataFrames
using IterTools
using DelimitedFiles

using Test

using Mimi: modelinstance, compinstance, get_var_value, OUTER, INNER, ReshapedDistribution

include("test-model/test-model.jl")
using .TestModel
m = create_model()

N = 100

sd = @defsim begin
    # Define random variables. The rv() is required to disambiguate an
    # RV definition name = Dist(args...) from application of a distribution
    # to an model parameter. This makes the (less common) naming of an
    # RV slightly more burdensome, but it's only required when defining
    # correlations or sharing an RV across parameters.
    rv(name1) = Normal(1, 0.2)
    rv(name2) = Uniform(0.75, 1.25)
    rv(name3) = LogNormal(20, 4)

    # assign RVs to model Parameters
    grosseconomy.share = Uniform(0.2, 0.8)
    emissions.sigma[:, Region1] *= name2
    emissions.sigma[2020:5:2050, (Region2, Region3)] *= Uniform(0.8, 1.2)

    grosseconomy.tester = ReshapedDistribution([20, 3], Dirichlet(20*3, 1))

    grosseconomy.depk = [Region1 => Uniform(0.08, 0.14),
            Region2 => Uniform(0.10, 1.50),
            Region3 => Uniform(0.10, 0.20)]

    # indicate which parameters to save for each model run. Specify
    # a parameter name or [later] some slice of its data, similar to the
    # assignment of RVs, above.
    save(grosseconomy.K, grosseconomy.YGROSS, emissions.E, emissions.E_Global)
end

output_dir = joinpath(tempdir(), "mcs")
si = run(sd, m, N; results_output_dir = output_dir, trials_output_filename = joinpath(output_dir, "trialdata.csv"))

d = readdlm(joinpath(output_dir, "trialdata.csv"), ',')
@test size(d)[1] == N+1 # extra row for column names


text:
export isbase58

function isbase58(str::AbstractString)::Bool
    base58Reg = r"^[A-HJ-NP-Za-km-z1-9]*$"

    return contains(str, base58Reg)
end


text:
import Base.vec

struct NominalState{T<:Real}
    p::Vector{T}
    v::Vector{T}
    q::Quaternion{T}
    ab::Vector{T}
    ωb::Vector{T}
    g::Vector{T}
end

function NominalState(x::AbstractVector)
    NominalState(
        x[1:3],   # position
        x[4:6],   # velocity
        Quaternion(x[7], x[8], x[9], x[10]), # attitude
        x[11:13], # accelerometer bias
        x[14:16], # gyrometer bias
        x[17:19]  # gravity vector
    )
end

function vec(nominal::NominalState)::AbstractVector
    [
        nominal.p
        nominal.v
        real(nominal.q)
        Quaternions.imag(nominal.q)
        nominal.ab
        nominal.ωb
        nominal.g
    ]
end

struct ErrorState{T<:Real}
    δp::Vector{T}
    δv::Vector{T}
    δθ::Vector{T}
    δab::Vector{T}
    δωb::Vector{T}
    δg::Vector{T}
end

function ErrorState(x::AbstractVector)
    ErrorState(
        x[1:3],   # position error
        x[4:6],   # velocity error
        x[7:9],   # angles error
        x[10:12], # accelerometer bias error
        x[13:15], # gyrometer bias error
        x[16:18]  # gravity error
    )
end

function vec(error::ErrorState)::AbstractVector
    [
        error.δp
        error.δv
        error.δθ
        error.δab
        error.δωb
        error.δg
    ]
end

struct Measurement{T<:Real}
    am::Vector{T}
    ωm::Vector{T}
end

struct Noise{T<:Real}
    σan::T # velocity noise Vi (m^2/s^2)
    σωn::T # orientation noise Θi (rad^2)
    σaw::T # accelero bias noise Ai (m^2/s^4)
    σωw::T # gyro bias noise Ωi (rad^2/s^2)
end


text:
type IC <: InfoCrit
   likelihood ::Float64
   df         ::Float64
   method     ::ASCIIString
   crterm     ::Vector{Float64}
   value      ::Float64
end

max_stars_repo_name: ilovealldata/SelectMixedModels

max_issues_repo_name: ilovealldata/SelectMixedModels

max_forks_repo_name: ilovealldata/SelectMixedModels


text:
# linear cone
line_mat_lc = LineBasicMaterial(color=color_lc, linewidth=10.0)
points_lc = Vector{Point{3,Float64}}()
for (i, xt) in enumerate(storage_lc.x[1])
    k = xt
    push!(points_lc, Point(k[1], k[2] + 0.04, k[3]))
end

setobject!(vis[:path_lc], Line(points_lc, line_mat_lc))

# nonlinear cone
line_mat_nc = Dojo.LineBasicMaterial(color=color_nc, linewidth=25.0)
points_nc = Vector{Point{3,Float64}}()
for (i, xt) in enumerate(storage_nc.x[1])
    k = xt
    push!(points_nc, Point(k[1], k[2] - 0.04, k[3] + 0.0))
end
setobject!(vis[:path_nc], MeshCat.Line(points_nc, line_mat_nc))

# MuJoCo
line_mat_mj = LineBasicMaterial(color=color_mj, linewidth=25.0)
points_mj = Vector{Point{3,Float64}}()
for (i, xt) in enumerate(storage_mj.x[1])
    k = xt
    push!(points_mj, Point(k[1], k[2], k[3]))
end
setobject!(vis[:path_mj], MeshCat.Line(points_mj, line_mat_mj))

settransform!(vis[:lc], MeshCat.Translation(0,+0.04,0))
settransform!(vis[:nc], MeshCat.Translation(0,-0.04,0))
settransform!(vis[:mj], MeshCat.Translation(0,+0.00,0))


text:
@testset "Testing modeltools.jl" begin

    @testset "Log Utility" begin
    # Test log utility
    u = LogUtility(2.0)
    @test isapprox(u(1.0), 0.0)
    @test isapprox(2.0*log(2.5), u(2.5))
    @test isapprox(u.(ones(5)), zeros(5))
    @test u(2.0) > u(1.0)  # Ensure it is increasing
    @test isapprox(derivative(u, 1.0), 2.0)
    @test isapprox(derivative(u, 2.0), 1.0)

    # Test "extrapolation evaluations"
    @test isapprox(u(0.5e-10), u(1e-10) + derivative(u, 1e-10)*(0.5e-10 - 1e-10))
    @test u(-0.5) < u(-0.1)  # Make sure it doesn't fail with negative values

    @test isapprox(LogUtility().ξ, 1.0)  # Test default constructor

    end

    @testset "CRRA Utility" begin
    # Test CRRA utility
    u = CRRAUtility(2.0)
    @test isapprox(u(1.0), 0.0)
    @test isapprox((2.5^(-1.0) - 1.0) / (-1.0), u(2.5))
    @test isapprox(u.(ones(5)), zeros(5))
    @test u(5.0) > u(3.0)  # Ensure it is increasing
    @test isapprox(derivative(u, 1.0), 1.0)
    @test isapprox(derivative(u, 2.0), 0.25)

    # Test "extrapolation evaluations"
    @test isapprox(u(0.5e-10), u(1e-10) + derivative(u, 1e-10)*(0.5e-10 - 1e-10))
    @test u(-0.5) < u(-0.1)  # Make sure it doesn't fail with negative values

    @test_throws ErrorException CRRAUtility(1.0)  # Test error throwing at γ=1.0
    end

    @testset "CFE Utility" begin
    # Test CFE Utility
    v = CFEUtility(2.0)
    @test isapprox(v(1.0), -2.0/3.0)
    @test isapprox(v(0.5), -v.ξ * 0.5^(1.0 + 1.0/v.ϕ) / (1.0 + 1.0/v.ϕ))
    @test v(0.5) > v(0.85)  # Ensure it is decreasing
    @test isapprox(derivative(v, 0.25), -0.5)
    @test isapprox(derivative(v, 1.0), -1.0)

    # Test "extrapolation evaluations"
    @test isapprox(v(0.5e-10), v(1e-10) + derivative(v, 1e-10)*(0.5e-10 - 1e-10))
    @test v(-0.5) > v(-0.1)  # Make sure it doesn't fail with negative values

    @test_throws ErrorException CRRAUtility(1.0)  # Test error throwing at ϕ=1.0
    end

    @testset "Elliptical Utility" begin
    # Test Elliptical Utility
    u = EllipticalUtility(1.0, 2.0)
    @test isapprox(u(sqrt(0.75)), 0.5)
    @test isapprox(u(sqrt(0.51)), 0.7)
    @test u(0.5) > u(0.95)  # Ensure it is decreasing
    @test isapprox(derivative(u, 0.0), 0.0)
    @test isapprox(derivative(u, sqrt(0.5)), -1.0)

    # Test default values
    @test isapprox(EllipticalUtility().b, 0.5223)
    @test isapprox(EllipticalUtility().μ, 2.2926)

    end

end

module Test_at_def_sim
    using QuantEcon
    using Base.Test

    @def_sim Simulation (T => Float64,) struct Observation{T<:Number}
        c::T
        k::T
        i_z::Int
    end

    @testset "@def_sim" begin
        @test isdefined(:Observation)
        @test isdefined(:Simulation)
        @test method_exists(Simulation, Tuple{NTuple})
        @test method_exists(Base.endof, Tuple{Simulation})
        @test method_exists(Base.length, Tuple{Simulation})
        @test method_exists(Base.start, Tuple{Simulation})
        @test method_exists(Base.next, Tuple{Simulation,Int})
        @test method_exists(Base.done, Tuple{Simulation,Int})
        @test method_exists(Base.getindex, Tuple{Simulation,Int})

        sim = Simulation([0.1, 0.2], [1.1, 1.2], [3, 4])
        @test isa(sim, Simulation{1,Float64})
        @test length(sim) == 2
        @test endof(sim) == 2
        obs12 = [Observation(0.1, 1.1, 3), Observation(0.2, 1.2, 4)]

        for (i, have) in enumerate(sim)
            @test have == obs12[i]
            @test obs12[i] == @inferred sim[i]
        end

        sim10 = @inferred Simulation((2, 1, 2, 1, 2, 1, 2, 1, 2, 1))
        @test isa(sim10, Simulation{10,Float64})
        @test length(sim10) == 32
        @test endof(sim10) == 32

        sim4 = @inferred Simulation(rand(4, 4, 4, 4), rand(4, 4, 4, 4), rand(Int, 4, 4, 4, 4))
        @test isa(sim4, Simulation{4,Float64})
        @test length(sim4) == 4*4*4*4
        @test endof(sim4) == 4*4*4*4

        sim2 = @inferred Simulation(rand(Float16, 4, 4), rand(Float16, 4, 4), rand(Int, 4,4))
        @test isa(sim2, Simulation{2,Float16})
        @test length(sim2) == 4*4
        @test endof(sim2) == 4*4
    end
end


text:
#
# Errors
#
# (many of) the errors thrown

####################################################################################################
#                                   Dimensions
empty_support_error(name::Symbol) = throw(ArgumentError(repr(name) * "'s support is empty"))

non_unique_support_error() = throw(ArgumentError("support must have unique values"))

not_in_dimension_error(v, d::Dimension) =
        throw(ArgumentError("value " * repr(v) * " is not in dimension " * repr(d)))

####################################################################################################
#                                   Factors

non_unique_dims_error() = throw(ArgumentError("Dimension names must be unique"))

not_enough_dims_error() =
        throw(DimensionMismatch("`potential`s number of " *
                "dimensions must match the length of `dimensions`"))

invalid_dim_sizes() = throw(DimensionMismatch("dimension lengths must match shape of `potential`"))

not_in_factor_error(name::Symbol, φ::Factor) =
        throw(ArgumentError(repr(name) * " is not in factor " * repr(φ)))

invalid_dims_error(func, got) =
        throw(TypeError(func, "type of dimensions", Union{Symbol, AbstractVector{Symbol}}, got))

assignment_missing_error(name::Symbol) =
        throw(ArgumentError("Assignment is missing variable " * repr(name)))

nonsingleton_assignment_error(name::Symbol) =
        throw(ArgumentError("Assignment to " * repr(name) * " should be a single value"))


text:
#!/usr/bin/env julia
using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, Pose2D, PoseStamped, Vector3, Twist
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan
@rosimport gazebo_msgs.msg: ModelState, ModelStates
@rosimport gazebo_msgs.srv: SetModelState, GetModelState, GetWorldProperties

rostypegen()
using geometry_msgs.msg
using std_srvs.srv
using nav_msgs.srv
using nav_msgs.srv.GetPlan
using gazebo_msgs.srv
using gazebo_msgs.msg

function callback(twist_msg::Twist, get_state::ServiceProxy{GetModelState}, set_state::ServiceProxy{SetModelState})
    modelName = RobotOS.get_param("robotName")

    # Get the current position of the Gazebo model
    gs = GetModelStateRequest()
    gs.model_name = modelName
    gs_r = get_state(gs)

    if !gs_r.success
        error(string(" calling /gazebo/get_model_state service: ", gs_r.status_message))
    end

    # Define the robot state
    ms = ModelState()
    ms.model_name = modelName
    ms.pose = gs_r.pose
    ms.twist = twist_msg

    # Set the state of the Gazebo model
    ss = SetModelStateRequest()
    ss.model_state = ms
    ss_r = set_state(ss)

    if !ss_r.success
        error(string(" calling /gazebo/set_model_state service: ", ss_r.status_message))
    end

end

function main()
    init_node("rosjl_example")

    # Set up service to set Gazebo model state
    const set_state = ServiceProxy("/gazebo/set_model_state", SetModelState)
    println("Waiting for 'gazebo/set_model_state' service...")
    wait_for_service("gazebo/set_model_state")

    # Set up service to get Gazebo model state
    const get_state = ServiceProxy("/gazebo/get_model_state", GetModelState)
    println("Waiting for 'gazebo/get_model_state' service...")
    wait_for_service("gazebo/get_model_state")

    sub = Subscriber{Twist}("/cmd_vel",callback,(get_state,set_state), queue_size=10)  # pub is the additional callback param is given to forward on
    spin()
end

if ! isinteractive()
    main()
end

#= TODO

# Teliop
function callback(msg::Twist, pub_obj::Publisher{ModelState)
    pt_msg = Pose(msg.linear.x, msg.linear.x, msg.linear.z, msg.angular.x, msg.angular.y, msg.angular.z)

    # Define the robot state
    ms = ModelState()
    ms.model_name = modelName
    ms.pose = vehPose

    publish(pub_obj, pt_msg)
end

function loop(pub_obj)
    loop_rate = Rate(5.0)
    while ! is_shutdown()
        npt = Point(rand(), rand(), 0.0)
        publish(pub_obj, npt)
        rossleep(loop_rate)
    end
end
function main()
    init_node("rosjl_example")
    pub = Publisher{"/gazebo/model_states", ModelState}
    sub = Subscriber{Twist}("/cmd_vel",callback,(pub,), queue_size=10)  # pub is the additional callback param is given to forward on
    loop(pub)
end

if ! isinteractive()
    main()
end

# Teliop
# wait for all services to be in some ready state, seehttp://docs.ros.org/electric/api/srs_user_tests/html/prepare__robot__manip__sim_8py_source.html
# start laucnh file in a paused state and have some sin_init file start it

 sim = get_param("/use_sim_time")
 if sim
      loginfo('Waiting until simulated robot is prepared for the task...')

      const get_world = ServiceProxy("/gazebo/get_world_properties",GetWorldProperties)
      println("Waiting for '/gazebo/get_world_properties' service...")
      wait_for_service("/gazebo/get_world_properties")

      wp_r=GetWorldPropertiesRequest()

      wp=get_world(wp_r)

wait_for_message('/sim_robot_init',EmptyMsg)
wait_for_service("/gazebo/pause_physics")

loginfo("Pausing physics")

=#
# world files with obstacles
# LiDAR data processing algorithm
# figure out elipses


text:
#=
    CanonicalIndexError (since Julia 1.8)
    https://github.com/JuliaLang/julia/pull/43852
=#

function _throw_canonical_error(msg::String, T::Type)
    @static if isdefined(Base, :CanonicalIndexError)
        throw(CanonicalIndexError(msg, T))
    else
        error(msg, T)
    end
end

@static if !isdefined(Base, :CanonicalIndexError)
    const CanonicalIndexError = ErrorException
end


text:
#testing seed is preserved

Fs = 4e4                    # Sampling frequency
c = 343                     # create new acoustic env with default values
Nt = round(Int64,0.1*Fs)    # Number of time samples, 0.25 secs
L  = (4., 4.,  4.)          # Room dimensions
Lx,Ly,Lz = L
xs = (2., 1.5, 1.)      # Source position
xr = (1., 2.,  2.)   # Receiver position

T60 = 0.1           # Reverberation Time
beta = ImageMethodReverb.revTime2beta(L, T60, c)
Tw = 20              # samples of Low pass filter
Fc = 0.9             # cut-off frequency

println("testing β for a given T60 gives correct T60")
@time h, seed  = rim(xs,xr,L,beta,Nt,Fs; N = (4,4,4))
rt, edc = revTime(h, Fs)
@test norm(rt.-T60)/norm(T60)< 0.25 #testing rev time

println("testing seed is preserved")
@time h2,  = rim(xs,xr,L,beta,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(h-h2) < 1e-8

println("testing calls")
# with T60
h2,  = rim(xs,xr,L,T60,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(h-h2) < 1e-8

# with source signal
s = randn(Nt)
p,  = rim(s,xs,xr,L,T60,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(p-conv(h,s)) < 1e-8

# multi mics
xr = [xr; (Lx-0.5,Ly,Lz)./2]
H,  = rim(xs,xr,L,beta,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(H[:,1]-h2) < 1e-8

# multi mics with source
P,  = rim(s,xs,xr,L,beta,Nt,Fs; N = (4,4,4), seed = seed)
@test norm(P[:,1]-p) < 1e-8

## original image source method
println("test with no fractional delay")
xs = (2., 1.5, 1.)      # Source position
xr = (1., 2.,  2.)   # Receiver position
h2, = rim(xs,xr,L,beta,Nt,Fs; Tw = 0, Fc = 0., N = (4,4,4), seed = seed)
@time h2, = rim(xs,xr,L,beta,Nt,Fs; Tw = 0, Fc = 0., N = (4,4,4), seed = seed)

println("test line of sight")
xs = (1., 1., 1.)
for i = 1:10
    d = (0.5 .*randn(3)...,) # dandom dispacement
    h_rand, = rim(xs, xs .+ d, L, 0.001, Nt,Fs; Rd = 1e-1, N = (0,0,0))
    @test argmax(h_rand)*1/Fs - norm(xr .- xs)/343 < 1e-4
end

# testing errors
@test_throws ErrorException rim(xs, (1,1,1.1), (1,1,1), beta, Nt,Fs)
@test_throws ErrorException rim(xs, (1,1, -1), (1,1,1), beta, Nt,Fs)
@test_throws ErrorException rim((1,1,1.1), xs,  (1,1,1), beta, Nt,Fs)
@test_throws ErrorException rim((1,1, -1), xs,  (1,1,1), beta, Nt,Fs)
@test_throws ErrorException rim(xs, xr, (1,1,1), beta, Nt,Fs; N = (-1,4,4))
@test_throws ErrorException rim(xr, xs,  (1,1,1), -0.1, Nt,Fs)


text:
# module adapters_mysql

module Options

arguments = (
    username = "root",
    password = "",
    hostname = "localhost",
    port = 3306,
    db = "mysqltest",
    unix_socket = Sys.islinux() ? "/var/run/mysqld/mysqld.sock" : "/tmp/mysql.sock",
)

end # module Options

# module adapters_mysql


text:
# SSigned
struct SInt128{V} <: Signed
    function SInt128{V}() where V
        !(typeof(V) === Int128) && throw(ArgumentError("SInt128 only supports static Int128 storage, got $(typeof(V))"))
        new{V}()
    end
end

struct SInt64{V} <: Signed
    function SInt64{V}() where V
        !(typeof(V) === Int64) && throw(ArgumentError("SInt64 only supports static Int64 storage, got $(typeof(V))"))
        new{V}()
    end
end

const SInt{V} = SInt64{V}

struct SInt32{V} <: Signed
    function SInt32{V}() where V
        !(typeof(V) === Int32) && throw(ArgumentError("SInt32 only supports static Int32 storage, got $(typeof(V))"))
        new{V}()
    end
end

struct SInt16{V} <: Signed
    function SInt16{V}() where V
        !(typeof(V) === Int16) && throw(ArgumentError("SInt16 only supports static Int16 storage, got $(typeof(V))"))
        new{V}()
    end
end

struct SInt8{V} <: Signed
    function SInt8{V}() where V
        !(typeof(V) === Int8) && throw(ArgumentError("SInt8 only supports static Int8 storage, got $(typeof(V))"))
        new{V}()
    end
end

const SSigned{V} = Union{SInt8{V},SInt16{V},SInt32{V},SInt64{V},SInt128{V}}

SSigned(x::Int8) = SInt8(x)
SSigned(x::Int16) = SInt16(x)
SSigned(x::Int32) = SInt32(x)
SSigned(x::Int64) = SInt64(x)
SSigned(x::Int128) = SInt128(x)

Base.widen(::Type{<:SInt8{X}}) where X = SInt16
Base.widen(::Type{<:SInt16{X}}) where X = SInt32
Base.widen(::Type{<:SInt32{X}}) where X= SInt64
Base.widen(::Type{<:SInt64{X}}) where X = SInt128
Base.widen(::Type{<:SInt128{X}}) where X = SInt128  # Note: there's no `SBigInt` so max out at 128

Base.Signed(x::SSigned) = x

function SSigned(val::Val{V}) where V
    if V isa Int8
        SInt8{V}()
    elseif V isa Int16
        SInt16{V}()
    elseif V isa Int32
        SInt32{V}()
    elseif V isa Int128
        SInt128{V}()
    else
        SInt64{V}()
    end
end

const AbstractInt128 = Union{SInt128,Int128}
const AbstractInt64  = Union{SInt64,Int64}
const AbstractInt32  = Union{SInt32,Int32}
const AbstractInt16  = Union{SInt16,Int16}
const AbstractInt8   = Union{SInt8,Int8}

SSI2BSI = Dict(SInt128 => Int128,
               SInt16 => Int16,
               SInt32 => Int32,
               SInt64 => Int64,
               SInt8 => Int8)

BSI2SSI = Dict(Int128 => SInt128,
               Int16 => SInt16,
               Int32 => SInt32,
               Int64 => SInt64,
               Int8 => SInt8)

for (ST,BT) in SSI2BSI
    defbasics(ST, BT)
    defbool(ST, BT)
end

SInt128One, SInt128Zero = defmath(SInt128, Int128)
SInt128NegOne = - SInt128One

SInt64One, SInt64Zero = defmath(SInt64, Int64)
const SZero = SInt64Zero
const SOne = SInt64One

SInt64NegOne = - SInt64One

SInt32One, SInt32Zero = defmath(SInt32, Int32)
SInt32NegOne = - SInt32One

SInt16One, SInt16Zero = defmath(SInt16, Int16)
SInt16NegOne = - SInt16One

SInt8One, SInt8Zero = defmath(SInt8, Int8)
SInt8NegOne = - SInt8One


text:
"""
    log_network(state, tick_nr)

Log the current state of a simulation.

# Arguments
- `state`: A tuple of the current graph and agent_list
- `tick_nr`: Current tick in the simulation

See also: [`tick!`](@ref), [`simulate`](@ref)
"""
function log_network(
    state::Tuple{AbstractGraph, AbstractArray},
    tick_nr::Int64
)

    graph, agent_list = state
    agent_opinion = [a.opinion for a in agent_list]
    agent_perceiv_publ_opinion = [a.perceiv_publ_opinion for a in agent_list]
    agent_check_regularity = [a.check_regularity for a in agent_list]
    agent_inclin_interact = [a.inclin_interact for a in agent_list]
    agent_desired_input_count = [a.desired_input_count for a in agent_list]
    agent_inactive_ticks = [a.inactive_ticks for a in agent_list]
    agent_active_state = [a.active for a in agent_list]
    agent_indegree = indegree(graph)
    agent_outdegree = outdegree(graph)

    return DataFrame(
        TickNr=tick_nr,
        AgentID=1:length(agent_list),
        Opinion=agent_opinion,
        PerceivPublOpinion=agent_perceiv_publ_opinion,
        CheckRegularity=agent_check_regularity,
        InclinInteract=agent_inclin_interact,
        DesiredInputCount=agent_desired_input_count,
        InactiveTicks=agent_inactive_ticks,
        Indegree=agent_indegree,
        Outdegree=agent_outdegree,
        ActiveState=agent_active_state
    )

end

# suppress output of include()
;


text:
using MinimalRLCore
using Random
using Test
using StableRNGs

mutable struct TestEnvEpis <: MinimalRLCore.AbstractEnvironment
    size::Int
    state::Int
end

MinimalRLCore.get_state(env::TestEnvEpis) = env.state
MinimalRLCore.get_reward(env::TestEnvEpis) = env.state == env.size ? 1.0 : 0.0
MinimalRLCore.get_actions(env::TestEnvEpis) = (1, 2)
MinimalRLCore.is_terminal(env::TestEnvEpis) = env.state == env.size

MinimalRLCore.reset!(env::TestEnvEpis, rng::Random.AbstractRNG=Random.GLOBAL_RNG) =
    env.state = rand(rng, 1:env.size)

MinimalRLCore.environment_step!(env::TestEnvEpis, action, rng::AbstractRNG=Random.GLOBAL_RNG) =
    env.state += action == 1 ? -1 : 1

mutable struct TestAgentEps <: AbstractAgent
    action::Int
end
# TestAgentEps(action) = TestAgentEps(action)

function MinimalRLCore.start!(agent::TestAgentEps, s, args...)
    agent.action
end

function MinimalRLCore.step!(agent::TestAgentEps, s, r, t, args...)
    agent.action
end

function test_episode()

    res = [(4, 2, 5, 0.0, false),
           (5, 2, 6, 0.0, false),
           (6, 2, 7, 0.0, false),
           (7, 2, 8, 0.0, false),
           (8, 2, 9, 0.0, false),
           (9, 2, 10, 1.0, true)]

    env = TestEnvEpis(10, 1)
    agent = TestAgentEps(2)
    @testset "Test Episode" begin
        ret = []
        rng = StableRNGs.StableRNG(1)
        total_rew, steps = MinimalRLCore.run_episode!(env, agent, 10000, rng) do (sars)
           push!(ret, sars)
        end
        println(ret)
        @test all(ret .== res)
        @test total_rew == 1.0
        @test steps == 6
    end

    @testset "Test Episode Max Steps" begin
        ret = []
        rng = StableRNGs.StableRNG(1)
        total_rew, steps = MinimalRLCore.run_episode!(env, agent, 2, rng) do (sars)
           push!(ret, sars)
        end
        @test all(ret .== res[1:2])
        @test total_rew == 0.0
        @test steps == 2
    end

end


text:
using LightGraphs
using SimpleWeightedGraphs
using LinearAlgebra
using GraphPlot
using Colors
using Serialization
using ..ChemistryFeaturization.Utils.GraphBuilding

# TO CONSIDER: store ref to featurization rather than the thing itself? Does this matter for any performance we care about?
"""
    AtomGraph

A type representing an atomic structure as a graph (`gr`).

# Fields
- `graph::SimpleWeightedGraph{<:Integer,<:Real}`: the graph representing the structure. See
  [`build_graph`](@ref) for more on generating the weights.
- `elements::Vector{String}`: list of elemental symbols corresponding to each node of the
  graph
- `laplacian::Matrix{<:Real}`: Normalized graph Laplacian matrix, stored to speed up
  convolution operations by avoiding recomputing it every pass.
- `encoded_atom_features::Union{Matrix{<:Real},Nothing}`: Feature matrix of size (# features, # nodes). AtomGraph can
  be initialized without defining this field, but if it is defined, the subsequent field must be also.
- `featurization::Union{AbstractFeaturization,Nothing}`: Featurization scheme specification to maintain
  "decodability" of features.
- `id::String`: Optional, an identifier, e.g. to correspond with tags/labels of an imported
  dataset.
"""
mutable struct AtomGraph <: AbstractAtoms
    graph::SimpleWeightedGraph{<:Integer,<:Real}
    elements::Vector{String}
    laplacian::Matrix{<:Real} # wanted to use LightGraphs.LinAlg.NormalizedGraphLaplacian but seems this doesn't support weighted graphs?
    encoded_features::Union{Matrix{<:Real},Nothing} # if we add edge features this type will have to relax
    featurization::Union{AbstractFeaturization,Nothing}
    id::String # or maybe we let it be a number too?
end

# first, the basic constructor
"""
    AtomGraph(gr, elements, features, featurization, id="")
    AtomGraph(gr, elements, id="")
    AtomGraph(adj, elements, features, featurization, id="")
    AtomGraph(adj, elements, id="")

Construct an AtomGraph object, either directly from a SimpleWeightedGraph `gr` or from an
adjacency matrix `adj`, along with, at minimum, the list of elemental symbols `elements`
representing each node. Note that the object can be initialized without features, but if
features are provided, so too must be the featurization scheme, in order to maintain
"decodability" of features.
"""
function AtomGraph(
    graph::SimpleWeightedGraph{A,B},
    elements::Vector{String},
    features::Matrix{<:Real},
    featurization::AbstractFeaturization,
    id = "",
) where {B<:Real,A<:Integer}
    # check that elements is the right length
    num_atoms = size(graph)[1]
    @assert length(elements) == num_atoms "Element list length doesn't match graph size!"

    # check that features is the right dimensions (# features x # nodes)
    expected_feature_length = sum(f.num_bins for f in featurization)
    @assert size(features) == (expected_feature_length, num_atoms) "Feature matrix is of wrong dimension! It should be of size (# features, # nodes)"

    # if all these are good, calculate laplacian and build the thing
    laplacian = normalized_laplacian(graph)
    AtomGraph(graph, elements, laplacian, features, featurization, id)
end

# one without features or featurization initialized yet
function AtomGraph(
    graph::SimpleWeightedGraph{A,B},
    elements::Vector{String},
    id = "",
) where {B<:Real,A<:Integer}
    # check that elements is the right length
    num_atoms = size(graph)[1]
    @assert length(elements) == num_atoms "Element list length doesn't match graph size!"

    laplacian = B.(normalized_laplacian(graph))
    AtomGraph(graph, elements, laplacian, nothing, nothing, id)
end

# initialize directly from adjacency matrix
AtomGraph(
    adj::Array{R1},
    elements::Vector{String},
    features::Matrix{R2},
    featurization::AbstractFeaturization,
    id = "",
) where {R2<:Real,R1<:Real} =
    AtomGraph(SimpleWeightedGraph{R}(adj), elements, features, featurization, id)

AtomGraph(adj::Array{R}, elements::Vector{String}, id = "") where {R<:Real} =
    AtomGraph(SimpleWeightedGraph(adj), elements, id)

function AtomGraph(
    input_file_path::String,
    id::String = "",
    output_file_path::Union{String,Nothing} = nothing,
    featurization::Union{AbstractFeaturization,Nothing} = nothing;
    overwrite_file::Bool = false,
    use_voronoi::Bool = false,
    cutoff_radius::Real = 8.0,
    max_num_nbr::Integer = 12,
    dist_decay_func::Function = inverse_square,
    normalize_weights::Bool = true,
)

    local ag
    local to_build_graph = true
    to_serialize = !isnothing(output_file_path)

    if to_serialize
        if isfile(output_file_path)
            if overwrite_file
                @info "Output file already exists and `overwrite_file` is set to false; returning deserialized AtomGraph at $output_file_path. If you wanted to rebuild the graph, set `overwrite=true`."
                to_build_graph = false
                ag = deserialize(output_file_path)
            end
        end
    end

    if to_build_graph
        if splitext(input_file_path)[end] == ".jls"
            ag = deserialize(input_file_path)
            ag.id = id
        else
            try
                ag = AtomGraph(
                    build_graph(
                        input_file_path,
                        use_voronoi = use_voronoi,
                        cutoff_radius = cutoff_radius,
                        max_num_nbr = max_num_nbr,
                        dist_decay_func = dist_decay_func,
                        normalize_weights = normalize_weights,
                    )...,
                    id,
                )
            catch
                @warn "Unable to build graph for $input_file_path"
                return missing
            end
        end
    end

    if !isnothing(featurization)
        featurize!(ag, featurization)
    end

    if to_serialize
        serialize(output_file_path, ag)
    end

    return ag
end

# pretty printing, short version
function Base.show(io::IO, ag::AtomGraph)
    st = "AtomGraph $(ag.id) with $(nv(ag.graph)) nodes, $(ne(ag.graph)) edges"
    if !isnothing(ag.featurization)
        st = string(st, ", feature vector length $(size(ag.encoded_features)[1])")
    end
    print(io, st)
end

# pretty printing, long version
function Base.show(io::IO, ::MIME"text/plain", ag::AtomGraph)
    st = "AtomGraph $(ag.id) with $(nv(ag.graph)) nodes, $(ne(ag.graph)) edges\n   atoms: $(ag.elements)\n   feature vector length: "
    if isnothing(ag.featurization)
        st = string(st, "uninitialized\n   featurization: uninitialized")
    else
        st = string(
            st,
            "$(size(ag.encoded_features)[1])\n   featurization: ",
            ag.featurization,
        )
    end
    print(io, st)

end

"""
    normalized_laplacian(graph)

Compute the normalized graph Laplacian matrix of the input graph, defined as

``I - D^{-1/2} A D^{-1/2}``

where ``A`` is the adjacency matrix and ``D`` is the degree matrix.
"""
function normalized_laplacian(g::G) where {G<:LightGraphs.AbstractGraph}
    a = adjacency_matrix(g)
    d = vec(sum(a, dims = 1))
    inv_sqrt_d = diagm(0 => d .^ (-0.5f0))
    laplacian = Float32.(I - inv_sqrt_d * a * inv_sqrt_d)
    !any(isnan, laplacian) || throw(
        ArgumentError(
            "NaN values in graph Laplacian! This is most likely due to atomic separations larger than the specified cutoff distance leading to block zeros in the adjacency matrix...try increasing the cutoff distance or inspecting your structure to ensure the file is correct.",
        ),
    )
    return laplacian
end

normalized_laplacian(ag::AtomGraph) = ag.laplacian

# now visualization stuff...

"Get a list of colors to use for graph visualization."
function graph_colors(atno_list, seed_color = colorant"cyan4")
    atom_types = unique(atno_list)
    atom_type_inds = Dict(atom_types[i] => i for i = 1:length(atom_types))
    color_inds = [atom_type_inds[i] for i in atno_list]
    colors = distinguishable_colors(length(atom_types), seed_color)
    return colors[color_inds]
end

"Helper function for sorting because edge ordering isn't preserved when converting to SimpleGraph."
function lt_edge(
    e1::SimpleWeightedGraphs.SimpleWeightedEdge{<:Integer,<:Real},
    e2::SimpleWeightedGraphs.SimpleWeightedEdge{<:Integer,<:Real},
)
    if e1.src < e2.src
        return true
    elseif e1.dst < e2.dst
        return true
    else
        return false
    end
end

"Compute edge widths (proportional to weights on graph) for graph visualization."
function graph_edgewidths(ag::AtomGraph)
    edgewidths = []
    edges_sorted = sort([e for e in edges(ag.graph)], lt = lt_edge)
    for e in edges_sorted
        append!(edgewidths, e.weight)
    end
    return edgewidths
end

"Visualize a given graph."
function visualize(ag::AtomGraph)
    # gplot doesn't work on weighted graphs
    sg = SimpleGraph(adjacency_matrix(ag))
    plt = gplot(
        sg,
        nodefillc = graph_colors(ag.elements),
        nodelabel = ag.elements,
        edgelinewidth = graph_edgewidths(ag),
    )
    display(plt)
end


text:
module SeqMaker

using JSON
using OpenGene
using OpenGene.Reference

const SEQMAKER_VERSION = "V1"

# package code goes here

include("compat.jl")
include("panel.jl")
include("sample.jl")
include("snv.jl")
include("cnv.jl")
include("fusion.jl")
include("sequencing.jl")
include("maker.jl")

# test
export ngs,
    alk,
    wgs

function ngs(outdir="output", panel_file="", profile_file=""; depth=0)
    if panel_file==""
        panel_file = joinpath(dirname(@__FILE__), "../data/panels/lung_cancer_hg19.bed")
    end
    if profile_file==""
        profile_file = joinpath(dirname(@__FILE__), "../data/profiles/example.json")
    end
    make_seq(panel_file, profile_file, outdir, depth)
end

function alk(profile_file=""; depth=50)
    if profile_file==""
        profile_file = joinpath(dirname(@__FILE__), "../data/profiles/fusion.json")
    end
    return ngs("alk", joinpath(dirname(@__FILE__), "../data/panels/alk.bed"), profile_file, depth=depth)
end

function wgs(outdir="wgs", profile_file=""; depth=30)
    return ngs(outdir, joinpath(dirname(@__FILE__), "../data/panels/hg19_whole_genome.bed"), profile_file, depth=depth)
end

end # module


text:
function _initialize_destination(destination, overwrite)
    if overwrite
        rm(destination; force = true, recursive = true)
    end
    if !isdir(destination)
        mkpath(destination)
    end
    return destination
end

function _append_to_path(path, thing_to_append)
  thing_to_append = replace(thing_to_append, r"[^0-9a-zA-Z]" => "")
  return joinpath(path, thing_to_append)
end

"""
    download_series(series_id::AbstractString, destination::AbstractString = "./", overwrite::Boolean = true)

Downloads images belonging to series with `series_id` and extracts them to `destination` folder.
If the destination folder already exists, then it will be overwritten by default unless `overwrite = false`.
"""
function download_series(series_id::AbstractString, destination = "./", overwrite = true)
    _initialize_destination(destination, overwrite)
    zip_file = joinpath(destination, "downloaded.zip")
    tcia_images(series = series_id, file = zip_file)
    unzip_command = `unzip -o $zip_file -d $destination`
    run(unzip_command)
    rm(zip_file)
    return destination
end

"""
    download_series(df::DataFrame, destination::AbstractString = "./"; append_desc::Boolean = true, overwrite::Boolean = true)

Downloads all images from the series in the dataframe `df` and then extracts them to `destination` folder.
The `df` can be obtained through the `tcia_series()` function.
By default, the series description will be appended to the path unless `append_desc = false`.
If the destination folder already exists, then it will be overwritten by default unless `overwrite = false`.
"""
function download_series(series_df::DataFrames.DataFrame, destination = "./"; append_desc = true, overwrite = true)
    return [download_series(row, destination; append_desc=append_desc, overwrite=overwrite) for row in eachrow(series_df)]
end

function download_series(series::DataFrames.DataFrameRow, destination = "./"; append_desc = true, overwrite = true)
    series_id = series.SeriesInstanceUID
    if append_desc
        destination = _append_to_path(destination, series.SeriesDescription)
    end
    return download_series(series_id, destination, overwrite)
end

"""
    download_series(arr::Array, destination::AbstractString = "./"; append_desc::Boolean = true, overwrite::Boolean = true)

Downloads all images from the series in the array `arr` and then extracts them to `destination` folder.
The `arr` can be obtained through the `tcia_series(..., format = "json")` command.
By default, the series description will be appended to the path unless `append_desc = false`.
If the destination folder already exists, then it will be overwritten by default unless `overwrite = false`.
"""
function download_series(series_array::Array, destination = "./"; append_desc = true, overwrite = true)
    return [download_series(series, destination; append_desc=append_desc, overwrite=overwrite) for series in series_array]
end

function download_series(series::Dict, destination = "./"; append_desc = true, overwrite = true)
    series_id = series["SeriesInstanceUID"]
    if append_desc
        destination = _append_to_path(destination, series["SeriesDescription"])
    end
    return download_series(series_id, destination, overwrite)
end


text:
# using these functions
using DelimitedFiles, Polynomials

# function definitions

# get user input for left baseline points
function get_baseline_points_l()
	println("Enter desired points for left baseline correction, separated by spaces:")
	string = readline()
	splitstring = split(string, " ")
	return [parse(Int, i) for i in splitstring]
end

# get user input for right baseline points
function get_baseline_points_r()
	println("Enter desired points for right baseline correction, separated by spaces:")
	string = readline()
	splitstring = split(string, " ")
	return [parse(Int, i) for i in splitstring]
end

# get order of polynomial for fitting
function get_polynomial_order()
	println("Enter polynomial order:")
	return parse(Int, readline())
end

# get rid of first 250 ps and return array of remaining spectra
function droptimes()
	println("getting input file")
	arr = readdlm("all_procesd_spectra_.csv", ',')
	pixels = @view arr[:, 1]
	longdelays = @view arr[:, 42:end]
	smallarr = hcat(pixels, longdelays)
	return smallarr
end

# stitch all times back together
function stitch(corrected)
	println("compiling output")
	original = readdlm("all_procesd_spectra_.csv", ',')
	shortdelays = @view original[2:end, 1:41]
	delays = @view original[1, :]
	arrout = vcat(transpose(delays), hcat(shortdelays, corrected))
	return arrout
end

# get left y-values for polynomial fitting
function ysforpolyfitl(arrin)
	spectral = @view arrin[2:129, 2:end]
	indexl = @view arrin[2:129, 1]
	yindiciesl = [findfirst(isequal(n), indexl) for n in fitpointsl]

	println("generate y-values for polynomial fit (left)")
	ys = [spectral[k, 1] for k in yindiciesl]
	for i in 2:length(spectral[1, : ])
		temp =  [spectral[j, i] for j in yindiciesl]
		ys = hcat(ys, temp)
	end
	return ys
end

# get right y-values for polynomial fitting
function ysforpolyfitr(arrin)
	spectrar = @view arrin[130:end, 2:end]
	indexr = @view arrin[130:end, 1]
	yindiciesr = [findfirst(isequal(n), indexr) for n in fitpointsr]

	println("generate y-values for polynomial fit (right)")
	ys = [spectrar[k, 1] for k in yindiciesr]
	for i in 2:length(spectrar[1, : ])
		temp =  [spectrar[j, i] for j in yindiciesr]
		ys = hcat(ys, temp)
	end
	return ys
end

# generate array of baseline subtractions for left detector
function genbasesubsl(arrin, polys)
	indexl = @view arrin[2:129, 1]
	println("generating baseline subtractions (left)")
	subslist = [polyval(polys[i], indexl) for i in 1:length(polys)]
	subs = subslist[1]
	for i in 2:length(subslist)
		subs = hcat(subs, subslist[i])
	end
	return subs
end

# generate array of baseline subtractions for right detector
function genbasesubsr(arrin, polys)
	indexr = @view arrin[130:end, 1]
	println("generating baseline subtractions (right)")
	subslist = [polyval(polys[i], indexr) for i in 1:length(polys)]
	subs = subslist[1]
	for i in 2:length(subslist)
		subs = hcat(subs, subslist[i])
	end
	return subs
end

# get array of left detector polynomials
function polynomialsl(yvalues)
	println("fitting polynomials (left)")
	polys = [polyfit(fitpointsl, yvalues[:, i], order) for i in 1:length(yvalues[1, : ])]
	return polys
end

# get array of right detector polynomials
function polynomialsr(yvalues)
	println("fitting polynomials (right)")
	polys = [polyfit(fitpointsr, yvalues[:, i], order) for i in 1:length(yvalues[1, : ])]
	return polys
end

# get array of corrected baseline corrected spectra
function subtractbase(arrin, subs)
	spectra = @view arrin[2:end, 2:end]
	println("subtracting baselines")
	corrected = spectra .- subs
	return corrected
end

# code below here
const fitpointsl = get_baseline_points_l()
const fitpointsr = get_baseline_points_r()
const order = get_polynomial_order()

const arrin = droptimes()
const subs = vcat(genbasesubsl(arrin, polynomialsl(ysforpolyfitl(arrin))), genbasesubsr(arrin, polynomialsr(ysforpolyfitr(arrin))))
const corrected = subtractbase(arrin, subs)

const outarr = stitch(corrected)

writedlm("testset.csv", outarr, ',')
println("complete")
exit()

max_stars_repo_name: mtmcrobie/ultra-proc

max_issues_repo_name: mtmcrobie/ultra-proc

max_forks_repo_name: mtmcrobie/ultra-proc


text:
using Chain
using Test

@testset "1" begin
    x = [1, 2, 3]
    # one symbol
    y = @chain x begin
        sum
    end
    @test y == sum(x)

    # two expressions
    z = @chain x begin
        *(3)
        sum
    end
    @test z == sum(x .* 3)

    # interleaved expressions
    called = false
    zz = @chain x begin
        .*(3)
        @aside @assert sum(_) / length(_) == 6 # this doesn't change anything
        @aside called = true # this also doesn't do the _ insertion and doesn't change anything
        sum
    end
    @test zz == z
    @test called

    zzz = @chain x begin
        _ .* 3
        sum
    end
    @test zzz == z
end

@testset "2" begin
    x = 1:4
    y = @chain x begin
        filter(isodd, _)
        map(-, _)
        sum
        _ ^ 2
    end
    @test y == 16
end

@testset "nested begin" begin
    x = 1:5
    y = @chain x begin
        begin
            z = sum(_) + 3
            z - 7
        end
        sqrt
    end
    @test y == sqrt(sum(x) + 3 - 7)
end

@testset "invalid invocations" begin
    # just one argument
    @test_throws LoadError eval(quote
        @chain [1, 2, 3]
    end)

    # no begin block
    @test_throws LoadError eval(quote
        @chain [1, 2, 3] sum
    end)

    # empty
    @test_throws LoadError eval(quote
        @chain [1, 2, 3] begin
        end
    end)

    # let block
    @test_throws LoadError eval(quote
        @chain [1, 2, 3] let
            sum
        end
    end)

    # variable defined in chain block doesn't leak out
    z = @chain [1, 2, 3] begin
        @aside inside_var = 5
        @aside @test inside_var == 5
        sum(_) + inside_var
    end
    @test z == 11
    @test_throws UndefVarError inside_var
end

@testset "nested chains" begin
    x = 1:5
    local z
    y = @chain x begin
        _ * 2
        @aside @chain _ begin
            sum(_)
            _ * 2
            @aside z = _
        end
        sum
    end
    @test y == sum(x * 2)
    @test z != x * 2
end

@testset "broadcast macro symbol" begin
    x = 1:5
    y = @chain x begin
        @. sin
        sum
    end
    @test y == sum(sin.(x))

    ## leave non-symbol invocations intact
    yy = @chain x begin
        @. sin(_)
        sum
    end
    @test yy == sum(sin.(x))
end

macro sin(exp)
    :(sin($(esc(exp))))
end

macro broadcastminus(exp1, exp2)
    :(broadcast(-, $(esc(exp1)), $(esc(exp2))))
end

@testset "splicing into macro calls" begin

    x = 1
    y = @chain x begin
        @sin
    end
    @test y == sin(x)

    y = @chain x begin
        @sin()
    end
    @test y == sin(x)

    xx = [1, 2, 3, 4]
    yy = @chain xx begin
        @broadcastminus(2.5)
    end
    @test yy == broadcast(-, xx, 2.5)

    xxx = [1, 2, 3, 4]
    yyy = @chain xxx begin
        @broadcastminus(2.5, _)
    end
    @test yyy == broadcast(-, 2.5, xxx)
end

@testset "single arg version" begin
    x = [1, 2, 3]

    xx = @chain begin
        x
    end
    @test xx == x

    # this has a different internal structure (one LineNumberNode missing I think)
    @test x == @chain begin
        x
    end

    @test sum(x) == @chain begin
        x
        sum
    end

    y = @chain begin
        x
        sum
    end
    @test y == sum(x)

    z = @chain begin
        x
        @. sqrt
        sum(_)
    end
    @test z == sum(sqrt.(x))

    @test sum == @chain begin
        sum
    end
end

@testset "invalid single arg versions" begin
    # empty
    @test_throws LoadError eval(quote
        @chain begin
        end
    end)

    # rvalue _ errors
    @test_throws ErrorException eval(quote
        @chain begin
            _
        end
    end)

    @test_throws ErrorException eval(quote
        @chain begin
            sum(_)
        end
    end)
end

@testset "handling keyword argments" begin
    f(a; kwarg) = (a, kwarg)
    @test (:a, :kwarg) == @chain begin
        :a
        f(kwarg = :kwarg)
    end
    @test (:a, :kwarg) == @chain begin
        :a
        f(; kwarg = :kwarg)
    end
end

# issue 13
@testset "no argument call" begin
    x = 1
    y = @chain x begin
        sin()
    end
    @test y == sin(x)
end

# issue 13
@testset "broadcasting calls" begin

    xs = [1, 2, 3]
    ys = @chain xs begin
        sin.()
    end
    @test ys == sin.(xs)

    add(x, y) = x + y

    zs = [4, 5, 6]
    sums = @chain xs begin
        add.(zs)
    end
    @test sums == add.(xs, zs)
end