text:
"""
    struct NystromDOF

Structure containing information of a degree-of-freedom in Nyström methods.
"""
struct NystromDOF{N,T,M,NM}
    # index::Int64         # QNode global index
    # element_index::Int64 # Index of the element to which QNode belongs
    coords::SVector{N,T}           # *Lifted* quadrature nodes
    weight::T                      # *Lifted* quadrature weigth
    jacobian::SMatrix{N,M,T,NM}    # Jacobian matrix at qnode
    locidx::Int                    # Index within element
    globidx::Int                   # Index in global system
end

coords(dof::NystromDOF)   = dof.coords
center(dof::NystromDOF)   = dof.coords
weight(dof::NystromDOF)   = dof.weight
jacobian(dof::NystromDOF) = dof.jacobian
normal(dof::NystromDOF)   = normal(jacobian(dof))

function NystromDOF(q::SVector{N,T},w::T,jac::SMatrix{N,_,T}) where {N,_,T}
    NystromDOF{N,T}(q,w,jac)
end

function Base.show(io::IO,dof::NystromDOF)
    print(io,
    "Nystrom DOF:
        coords  = $(dof.coords)
        locidx  = $(dof.locidx)
        globidx = $(dof.globidx)")
end

Interpolation.HyperRectangle(qnodes::Vector{<:NystromDOF}) = HyperRectangle(coords(q) for q in qnodes)

"""
    struct NystromMesh{N,T} <: AbstractMesh{N,T}

A mesh data structure for solving boundary integral equation using Nyström
methods.

A `NystromMesh` can be constructed from a `mesh::AbstractMesh` and a dictionary
`etype2qrule` mapping element types in `mesh` (the keys) to an appropriate
quadrature rule for integration over elements of that type (the value).

The degrees of freedom in a `NystromMesh` are associated to nodal values at the
quadrature nodes, and are reprensented using [`NystromDOF`](@ref).
"""
struct NystromMesh{N,T,M,NM} <: AbstractMesh{N,T}
    elements::Dict{DataType,Any}                                # Geometrical elements
    etype2qrule::Dict{DataType,AbstractQuadratureRule}          # Element type --> quadrature
    dofs::Vector{NystromDOF{N,T,M,NM}}                          # Degrees of freedom
    elt2dof::Dict{DataType,Matrix{Int}}                         # Elements --> dofs
    ent2elt::Dict{AbstractEntity,Dict{DataType,Vector{Int}}}    # Entitiy --> element
end

function Base.show(io::IO,msh::NystromMesh)
    print(io," NystromMesh with $(length(dofs(msh))) DOFs")
end

# getters
elements(m::NystromMesh)             = m.elements
elements(m::NystromMesh,E::DataType) = m.elements[E]
etype2qrule(m)   = m.etype2qrule
etype2qrule(m,E) = etype2qrule(m)[E]
dofs(m::NystromMesh)  = m.dofs
elt2dof(m::NystromMesh) = m.elt2dof
elt2dof(m::NystromMesh,E::DataType) = m.elt2dof[E]
ent2elt(m::NystromMesh) = m.ent2elt
ent2elt(m::NystromMesh,ent::AbstractEntity) = m.ent2elt[ent]

# generators for iterating over fields of dofs
qcoords(m::NystromMesh)  = (coords(q) for q in dofs(m))
qweights(m::NystromMesh) = (weight(q) for q in dofs(m))
qjacobians(m::NystromMesh) = (jacobian(q) for q in dofs(m))
qnormals(m::NystromMesh) = (normal(q) for q in dofs(m))

# entities and domain
Geometry.entities(mesh::NystromMesh) = collect(keys(mesh.ent2elt))
domain(mesh::NystromMesh)   = Domain(entities(mesh))

Base.keys(m::NystromMesh) = keys(elements(m))

# other maps constructed from ent2elt and elt2dof
function dom2elt(mesh::NystromMesh,domain::Domain)
    dict = Dict{DataType,Vector{Int}}()
    for ent in entities(domain)
        other = ent2elt(mesh,ent)
        mergewith!(append!,dict,other)
    end
    return dict
end

function dom2dof(mesh::NystromMesh,domain::Domain)
    idxs = Int[]
    for ent in entities(domain)
        dict = ent2elt(mesh,ent)
        for (E,tags) in dict
            append!(idxs,view(elt2dof(mesh,E),:,tags))
        end
    end
    return idxs
end
dom2dof(mesh,ent::AbstractEntity) = dom2dof(mesh,Domain(ent))

ent2dof(mesh,ent::AbstractEntity) = dom2dof(mesh,ent)

function NystromMesh(msh::AbstractMesh{N,T},Ω::Domain,e2qrule) where {N,T}
    # initialize mesh with empty fields
    M       = geometric_dimension(Ω) |> Int
    NM      = N*M
    nys_msh = NystromMesh{N,T,M,NM}(
        Dict{DataType,Any}(),                               # elemements
        e2qrule,                                            # etype2qrule
        NystromDOF{N,T,M,NM}[],                             # dofs
        Dict{DataType,Matrix{Int}}(),                       # elt2dof
        Dict{AbstractEntity,Dict{DataType,Vector{Int}}}()   # ent2elt
    )
    # loop over entities, generate quadrature, and fill in various fields
    for ent in entities(Ω)
        ent_msh = view(msh,ent)
        dict    = Dict{DataType,Vector{Int}}() # store ent2elt for current ent
        for E in keys(ent_msh) # element types
            @assert haskey(e2qrule,E) "no quadrature rule found for element of type $E"
            qrule  = e2qrule[E]
            els    = get!(nys_msh.elements,E,Vector{E}())
            istart = length(els)+1
            iter   = ElementIterator(ent_msh,E)
            _build_nystrom_mesh!(nys_msh,iter,qrule)
            iend    = length(els)
            dict[E] = collect(istart:iend)
        end
        nys_msh.ent2elt[ent] = dict # new entry
    end
    return nys_msh
end
function NystromMesh(msh::GenericMesh, Ω::Domain;order)
    etypes = keys(view(msh,Ω))
    e2qrule = Dict(E=>qrule_for_reference_shape(domain(E),order) for E in etypes)
    NystromMesh(msh,Ω,e2qrule)
end
function NystromMesh(submesh::SubMesh,args...;kwargs...)
    NystromMesh(parent(submesh),domain(submesh),args...;kwargs...)
end

@noinline function _build_nystrom_mesh!(msh,iter::ElementIterator,qrule::AbstractQuadratureRule)
    E               = eltype(iter)
    els::Vector{E}  = msh.elements[E]
    el2dofs         = get!(msh.elt2dof,E,Matrix{Int}(undef,0,0))
    el2dofs         = vec(el2dofs)
    x̂,ŵ             = qrule() #nodes and weights on reference element
    num_nodes       = length(ŵ)
    for el in iter
        # add element
        push!(els,el)
        # and all qnodes for that element
        for i in 1:num_nodes
            x       = el(x̂[i])
            jac     = jacobian(el,x̂[i])
            μ       = integration_measure(jac)
            w       = μ*ŵ[i]
            # global index of the dof that is going to be created
            globidx = length(dofs(msh)) + 1
            locidx  = i
            dof     = NystromDOF(x,w,jac,locidx,globidx)
            push!(dofs(msh),dof)
            push!(el2dofs,globidx)
        end
    end
    el2dofs        = reshape(el2dofs,num_nodes,:)
    msh.elt2dof[E] = el2dofs
    return msh
end

# ElementIterator for NystromMesh
function Base.size(iter::ElementIterator{<:AbstractElement,<:NystromMesh})
    E = eltype(iter)
    M = mesh(iter)
    return size(elements(M,E))
end

function Base.getindex(iter::ElementIterator{<:AbstractElement,<:NystromMesh},i::Int)
    E = eltype(iter)
    M = mesh(iter)
    return elements(M,E)[i]
end

function Base.iterate(iter::ElementIterator{<:AbstractElement,<:NystromMesh}, state=1)
    state > length(iter) && (return nothing)
    iter[state], state + 1
end

"""
    near_interaction_list(pts,Y;atol)

Given target points `pts` and a `Y::NystromMesh`, return a dictionary with keys
given by element types of the source mesh `Y`. Each key has a value given by a
vector whose `i`-th entry encodes information about the points in `X` which are
close to a given element of the key type.
"""
function near_interaction_list(pts,Y::NystromMesh;atol)
    dict = Dict{DataType,Vector{Vector{Tuple{Int,Int}}}}()
    for E in keys(Y)
        ilist = _etype_near_interaction_list(pts,Y,E,atol)
        push!(dict,E=>ilist)
    end
    return dict
end

function _etype_near_interaction_list(pts,Y,E,atol)
    ilist    = Vector{Vector{Tuple{Int,Int}}}()
    e2n      = elt2dof(Y,E)
    npts,nel = size(e2n)
    for n in 1:nel
        ynodes = @views dofs(Y)[e2n[:,n]]
        inear  = _near_interaction_list(pts,ynodes,atol)
        push!(ilist,inear)
    end
    ilist
end

function _near_interaction_list(pts,ynodes,atol)
    ilist    = Tuple{Int,Int}[]
    for (i,x) in enumerate(pts)
        d,j = findmin([norm(coords(x)-coords(y)) for y in ynodes])
        if d ≤ atol
            push!(ilist,(i,j))
        end
    end
    return ilist
end

# function Base.append!(msh1::NystromMesh,msh2::NystromMesh)
#     Ω = domain(msh2)
#     @assert ambient_dimension(msh1) == ambient_dimension(msh2)
#     @assert eltype(msh1) == eltype(msh2)
#     # extract relevant quadrature rules
#     etype2dof = Dict{DataType,Vector{Int}}()
#     for (E,q) in etypes(msh2,Ω)
#         if haskey(msh1.elements,E)
#             etype2dof[E]        = msh1.elt2dof[E] |> vec
#             @assert msh1.etype2qrule[E] == q
#         else
#             etype2dof[E]        = Int[]
#             msh1.elements[E]    = E[]
#             msh1.etype2qrule[E] = q
#         end
#     end
#     # loop over entities and extract the information
#     for ent in entities(msh2)
#         haskey(msh1.ent2elt,ent) && continue # entity already meshed
#         msh1.ent2elt[ent] = Dict{DataType,Vector{Int}}()
#         for (E,v) in ent2elt(msh2,ent)
#             # add dofs information
#             msh2_dofs  = msh2.elt2dof[E][:,v]
#             idx_start = length(msh1.qnodes) + 1
#             append!(msh1.qnodes,msh2.qnodes[msh2_dofs])
#             append!(msh1.qweights,msh2.qweights[msh2_dofs])
#             if !isempty(qnormals(msh2)) # otherwise assume it did not have normal stored
#                 append!(msh1.qnormals,msh2.qnormals[msh2_dofs])
#             end
#             idx_end = length(msh1.qnodes) # end index of dofs added
#             append!(etype2dof[E],collect(idx_start:idx_end))
#             # append elements
#             idx_start = length(msh1.elements[E]) + 1
#             append!(msh1.elements[E],msh2.elements[E][v])
#             idx_end   = length(msh1.elements[E])
#             # and mapping ent -> elt
#             msh1.ent2elt[ent][E] = collect(idx_start:idx_end)
#         end
#     end
#     # finally do some rehshaping
#     for (E,v) in etype2dof
#         dof_per_el = size(msh2.elt2dof[E],1)
#         msh1.elt2dof[E] = reshape(v,dof_per_el,:)
#     end
#     return msh1
# end

function diagonal_ncross_and_jacobian_matrices(nmesh)
    qnodes = dofs(nmesh)
    nmatrix = Diagonal([cross_product_matrix(normal(q)) for q in qnodes])
    jmatrix = Diagonal([jacobian(q) for q in qnodes])
    dual_jmatrix = Diagonal([SMatrix{2,3,Float64,6}(pinv(jacobian(q))) for q in qnodes])
    return nmatrix, jmatrix, dual_jmatrix
end
diagonal_ncross_matrix(nmesh) = diagonal_ncross_and_jacobian_matrices(nmesh)[1]
diagonal_jacobian_matrix(nmesh) = diagonal_ncross_and_jacobian_matrices(nmesh)[2]
diagonal_dualjacobian_matrix(nmesh) = diagonal_ncross_and_jacobian_matrices(nmesh)[3]

# Plot recipes

struct PlotQuadrature end

@recipe function f(mesh::NystromMesh)
    # plot elements
    label --> ""
    grid  --> false
    aspect_ratio --> :equal
    for (E,els) in mesh.elements
        for el in els
            @series begin
                el
            end
        end
    end
    # # plot quadrature
    # PlotQuadrature(),mesh
end

@recipe function f(::PlotQuadrature,mesh::NystromMesh)
    marker --> :circle
    pts = qcoords(mesh) |> collect
    WavePropBase.IO.PlotPoints(), pts
end

max_stars_repo_name: WaveProp/Nystrom

max_issues_repo_name: WaveProp/Nystrom

max_forks_repo_name: WaveProp/Nystrom


text:
##############################################################################
#
# Metadata field getters and setters
#
##############################################################################

import Languages.name

name(d::AbstractDocument) = d.metadata.name
language(d::AbstractDocument) = d.metadata.language
author(d::AbstractDocument) = d.metadata.author
timestamp(d::AbstractDocument) = d.metadata.timestamp

function name!(d::AbstractDocument, nv::AbstractString)
    d.metadata.name = nv
end

function language!{T <: Language}(d::AbstractDocument, nv::Type{T})
    d.metadata.language = nv
end

function author!(d::AbstractDocument, nv::AbstractString)
    d.metadata.author = nv
end

function timestamp!(d::AbstractDocument, nv::AbstractString)
    d.metadata.timestamp = nv
end

##############################################################################
#
# Vectorized getters for an entire Corpus
#
# TODO: Add vectorized setters
#
##############################################################################

names(c::Corpus) = map(d -> name(d), documents(c))
languages(c::Corpus) = map(d -> language(d), documents(c))
authors(c::Corpus) = map(d -> author(d), documents(c))
timestamps(c::Corpus) = map(d -> timestamp(d), documents(c))

ext: jl


text:
# Autogenerated wrapper script for SDL2_mixer_jll for x86_64-apple-darwin
export libsdl2_mixer

using SDL2_jll
using libvorbis_jll
using FLAC_jll
using mpg123_jll
JLLWrappers.@generate_wrapper_header("SDL2_mixer")
JLLWrappers.@declare_library_product(libsdl2_mixer, "@rpath/libSDL2_mixer-2.0.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(SDL2_jll, libvorbis_jll, FLAC_jll, mpg123_jll)
    JLLWrappers.@init_library_product(
        libsdl2_mixer,
        "lib/libSDL2_mixer.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Statistics
using Plots
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using JLD2
using Random
using Distributions
using FITSIO
using LineSearches
using SparseArrays

push!(LOAD_PATH, pwd()*"/../../../main")
using DHC_2DUtils
push!(LOAD_PATH, pwd()*"/../../../scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using ReconFuncs

numfile = Base.parse(Int, ENV["SLURM_ARRAY_TASK_ID"])
println(numfile, ARGS[1], ARGS[2])

if ARGS[1]=="log"
    logbool=true
else
    if ARGS[1]!="reg" error("Invalid log arg") end
    logbool=false
end

if ARGS[2]=="apd"
    apdbool = true
else
    if ARGS[2]!="nonapd" error("Invalid apd arg") end
    apdbool=false
end

if ARGS[3]=="iso"
    isobool = true
else
    if ARGS[3]!="noiso" error("Invalid iso arg") end
    isobool=false
end

direc = ARGS[4]
datfile = direc * "Data_" * string(numfile) * ".jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]
Nx=size(true_img)[1]
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
Nf = size(filter_hash["filt_index"])[1]
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>apdbool, :iso=>isobool)
if dhc_args[:iso]
    coeff_mask = falses(2+S1iso+S2iso)
    coeff_mask[S1iso+3:end] .= true
else #Not iso
    coeff_mask = falses(2+Nf+Nf^2)
    coeff_mask[Nf+3:end] .= triu(trues(Nf, Nf))[:]
end
recon_settings = Dict([("log", logbool), ("covar_type", "sfd_dbn"), ("target_type", "sfd_dbn"), ("Invcov_matrix", ARGS[6]), ("epsvalue", 1e-10)])
println("Apdbool", apdbool)
println("Isobool", isobool)
println(recon_settings)

s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings, safety=nothing)
invvar1 = diag(s2icov)

#Check wrt calculation here
dbnimg = readsfd(Nx, logbool=logbool)
dbn = get_dbn_coeffs(dbnimg, filter_hash, dhc_args, coeff_mask = nothing)
invvar2 = std(dbn, dims=1)[:]
println(size(invvar2))
invvar2 = (invvar2.^(-2))[coeff_mask]
println("InvVar1", invvar1[:10])
println("InvVar2", invvar2[:10])
println("Mean Abs Error = ", mean(abs.(invvar1 .- invvar2)))
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
module TestReflection
using Test
using Configurations

@option struct OptionA{T}
    x::T = 2
end

@option struct OptionB
    x::Int = 2
end

@testset "field_default" begin
    @testset "error case" begin
        @test_throws ErrorException field_default(ComplexF32, :im)
        @test_throws ErrorException field_default(OptionA, :x)
        @test_throws ErrorException field_default(OptionB, :y)
    end
end

@testset "type_alias" begin
    @test_throws ErrorException type_alias(Complex)
end

end


text:
import ..Common.test_constructor

function run_test_constructor()
    test_constructor(
        Dict(
            :atoms => [Atom(index=10)]
        ),
        Molecule([Atom(index=10)])
    )
end


text:
int f(int x, int x) {
   return x;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
@info "Loading packages..."
@time begin
    using VisualStringDistances, UnbalancedOptimalTransport
    using AbstractPlotting: px
    using UnbalancedOptimalTransport: KL, Balanced
    using Makie, MakieLayout
    using GeometryBasics: Point2f0
    using PlotUtils: RGBA, RGB
end

function hide_decorations!(ax)
    ax.xticksvisible=false
    ax.yticksvisible=false
    ax.xticklabelsvisible=false
    ax.yticklabelsvisible=false
    ax.bottomspinevisible = false
    ax.leftspinevisible = false
    ax.topspinevisible = false
    ax.rightspinevisible = false
    ax.xgridvisible=false
    ax.ygridvisible=false
end

function imgrot(v)
    Point2f0(v[2], 1-v[1])
end

# transparent to black for 0 to 1, then becomes redder from 1 to 2
const CMAP = cgrad([RGBA{Float64}(0.0,0.0,0.0,0.0), RGBA{Float64}(0.0,0.0,0.0,1.0), RGBA{Float64}(1.0,0.0,0.0,1.0)], [0,1,1.5])

density_to_color(d) = get(CMAP, d/2)

function animate_words(string1, string2;
    normalize_density = false,
    D = KL(1.0),
    random_colors = false,
    kwargs...
)
    w1 = word_measure(string1)
    w2 = word_measure(string2)

    if normalize_density
        w1 = DiscreteMeasure(w1.density / sum(w1.density), w1.set)
        w2 = DiscreteMeasure(w2.density / sum(w2.density), w2.set)
    end

    if random_colors
        # doesn't matter how we chose `π`
        π = ones(length(w1.set), length(w2.set))
    else
        π = optimal_coupling!(D, w1, w2)
    end

    scene, layout = layoutscene()
    layout[1,1] = ax = LAxis(scene)
    display(scene)
    animate_coupling!(scene, ax, π, w1.set, w2.set; random_colors = random_colors, kwargs...)
    return
end

function animate_coupling!(scene, ax, π, coords1, coords2;
    duration = 2,
    total_frames = 50*duration,
    pause_frames = total_frames ÷ 4,
    move_frames = total_frames - 2*pause_frames,
    markersize = 20px,
    random_colors = false,
    save_path = nothing,
    α = 0.7,
    compression = 20
)

    x_min = min(minimum([x[1] for x in imgrot.(coords1)]), minimum([x[1] for x in imgrot.(coords2)]))
    x_max = max(maximum([x[1] for x in imgrot.(coords1)]), maximum([x[1] for x in imgrot.(coords2)]))
    y_min = min(minimum([x[2] for x in imgrot.(coords1)]), minimum([x[2] for x in imgrot.(coords2)]))
    y_max = max(maximum([x[2] for x in imgrot.(coords1)]), maximum([x[2] for x in imgrot.(coords2)]))

    hide_decorations!(ax)

    s1, s2 = size(π)
    if random_colors
        cvals = collect(Iterators.Flatten(Iterators.repeated([ RGBA(rand(RGB), α) for _ = 1:s1], s2)))
    else
        scale = sqrt(prod(size(π)))
        cvals = density_to_color.(scale * vec(π) ./ sum(π))
    end

    t = 0
    locations = Node([ imgrot( (1-t)*coords1[i] + t*coords2[j]) for j = 1:size(π,2) for i = 1:size(π,1)])

    anim_plt = scatter!(ax, locations, color = cvals, markersize=markersize,
    transparency=true)
    x_pad = (x_max - x_min)*.05
    y_pad = (y_max - y_min)*.05
    limits!(ax, x_min - x_pad, x_max + x_pad, y_min - y_pad, y_max + y_pad)

    move_ts = range(0, 1; length=move_frames)
    do_frame = function(j)
        if (j <= pause_frames) || (j > total_frames - pause_frames)
            save_path === nothing && sleep(duration/total_frames)
        else
            t = move_ts[j - pause_frames]
            locations[] = [ imgrot((1-t)*coords1[i] + t*coords2[j]) for j = 1:size(π,2) for i = 1:size(π,1)]
            save_path === nothing && sleep(duration/total_frames)
        end
    end

    if save_path === nothing
        map(do_frame, 1:total_frames)
    else
        record(do_frame, scene, save_path, 1:total_frames, framerate = round(Int, total_frames / duration), compression=compression)
    end
    return
end

@eval AbstractPlotting begin
   function save(path::String, io::VideoStream;
              framerate::Int = 24, compression = 20)
    close(io.process)
    wait(io.process)
    p, typ = splitext(path)
    if typ == ".mkv"
        cp(io.path, path, force=true)
    elseif typ == ".mp4"
        ffmpeg_exe(`-loglevel quiet -i $(io.path) -crf $compression -c:v libx264 -preset slow -r $framerate -pix_fmt yuv420p -c:a libvo_aacenc -b:a 128k -y $path`)
    elseif typ == ".webm"
        ffmpeg_exe(`-loglevel quiet -i $(io.path) -crf $compression -c:v libvpx-vp9 -threads 16 -b:v 2000k -c:a libvorbis -threads 16 -r $framerate -vf scale=iw:ih -y $path`)
    elseif typ == ".gif"
        filters = "fps=$framerate,scale=iw:ih:flags=lanczos"
        palette_path = dirname(io.path)
        pname = joinpath(palette_path, "palette.bmp")
        isfile(pname) && rm(pname, force = true)
        ffmpeg_exe(`-loglevel quiet -i $(io.path) -vf "$filters,palettegen" -y $pname`)
        ffmpeg_exe(`-loglevel quiet -i $(io.path) -r $framerate -f image2 $(palette_path)/image_%06d.png`)
        ffmpeg_exe(`-loglevel quiet -framerate $framerate -i $(palette_path)/image_%06d.png -i $pname -lavfi "$filters [x]; [x][1:v] paletteuse" -y $path`)
        rm(pname, force = true)
    else
        rm(io.path)
        error("Video type $typ not known")
    end
    rm(io.path)
    return path
end
end

@info "Generating animations..."

@time begin
    animate_words("hello", "heIIo"; D = KL(1.0), save_path=abspath(joinpath(@__DIR__, "..", "..", "docs", "assets", "hello_heIIo.gif")))
    animate_words("hello", "heIIo"; D = Balanced(), normalize_density=true, save_path=abspath(joinpath(@__DIR__, "..", "..", "docs", "assets", "hello_heIIo_balanced.gif")))
end

@info "Done!"


text:
##########
##########
# A check of the distribution of shapes of biomass and richness
# in relevant paper
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-05-17
##########
##########

# set up =======================================================================
using DrWatson
@quickactivate "trait-based-rewiring"

# set packages
using DataFrames, Plots, StatsPlots

include(srcdir("01_eqn_component_functions.jl"))

# set all positional values
organism_id_pos = 1
organism_type_pos = 2
body_size_pos = 3
lower_prey_limit_pos = 4
upper_prey_limit_pos = 5
habitat_midpoint_pos = 6
habitat_lower_limit_pos = 7
habitat_upper_limit_pos = 8
nutritional_value_pos = 9
biomass_pos = 10
network_id_pos = 11
alive_pos = 12

# perform loop to grab all values ==============================================

# draw 100 networks of the biomass/richness shapes and look at how they group
richness_biomass_comp = Array{Any}(undef, (450000, 5))
j = 1
for i in 1:1000
    for richness in ["top-heavy", "bottom-heavy", "uniform"]
        for biomass in ["pyramid", "uniform", "inverted"]
            temp = Array{Any}(undef, (50, 5))
            species_list = network_draw(50, richness)[1]
            species_list = initialize_biomass(species_list, biomass)
            for k in 1:size(species_list, 1)
                temp[k, 1] = richness
                temp[k, 2] = biomass
                temp[k, 3] = species_list[k, body_size_pos]
                temp[k, 4] = species_list[k, habitat_midpoint_pos]
                temp[k, 5] = species_list[k, organism_type_pos]
            end
            richness_biomass_comp[(j:j+49), :] = temp
            j = j + 50
        end
    end
end

# turn array into dataframe
richness_biomass_comp = DataFrame(
    richness_biomass_comp,
    [:richness_shape, :biomass_shape, :body_size, :biomass, :org_type]
)

# make plots of the 9 combos ===================================================

# top-heavy & pyramid
topheavy_pryamid_data = DataFrame(filter(
    (row -> row.richness_shape == "top-heavy" &&
    row.biomass_shape == "pyramid"),
    richness_biomass_comp
))
topheavy_pryamid_plot = @df topheavy_pryamid_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = top-heavy, biomass = pyramid"
)
# top-heavy & uniform
topheavy_uniform_data = filter(
    (row -> row.richness_shape == "top-heavy" &&
    row.biomass_shape == "uniform"),
    richness_biomass_comp
)
topheavy_uniform_plot = @df topheavy_uniform_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = top-heavy, biomass = uniform"
)
# top-heavy and inverted
topheavy_inverted_data = filter(
    (row -> row.richness_shape == "top-heavy" &&
    row.biomass_shape == "inverted"),
    richness_biomass_comp
)
topheavy_inverted_plot = @df topheavy_inverted_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = top-heavy, biomass = inverted"
)
# uniform and pyramid
uniform_pyramid_data = filter(
    (row -> row.richness_shape == "uniform" &&
    row.biomass_shape == "pyramid"),
    richness_biomass_comp
)
uniform_pyramid_plot = @df uniform_pyramid_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = uniform, biomass = pyramid"
)
# uniform and uniform
uniform_uniform_data = filter(
    (row -> row.richness_shape == "uniform" &&
    row.biomass_shape == "uniform"),
    richness_biomass_comp
)
uniform_uniform_plot = @df uniform_uniform_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = uniform, biomass = uniform"
)
# uniform and inverted
uniform_inverted_data = filter(
    (row -> row.richness_shape == "uniform" &&
    row.biomass_shape == "inverted"),
    richness_biomass_comp
)
uniform_inverted_plot = @df uniform_inverted_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = uniform, biomass = inverted"
)
# bottom-heavy and pyramid
bottomheavy_pryamid_data = filter(
    (row -> row.richness_shape == "bottom-heavy" &&
    row.biomass_shape == "pyramid"),
    richness_biomass_comp
)
bottomheavy_pryamid_plot = @df bottomheavy_pryamid_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = bottom-heavy, biomass = pyramid"
)
# bottom-heavy and uniform
bottomheavy_uniform_data = filter(
    (row -> row.richness_shape == "bottom-heavy" &&
    row.biomass_shape == "uniform"),
    richness_biomass_comp
)
bottomheavy_uniform_plot = @df bottomheavy_uniform_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    legend = false
    #title = "Richness = bottom-heavy, biomass = uniform"
)
# bottom-heavy and inverted
bottomheavy_inverted_data = filter(
    (row -> row.richness_shape == "bottom-heavy" &&
    row.biomass_shape == "inverted"),
    richness_biomass_comp
)
bottomheavy_inverted_plot = @df bottomheavy_inverted_data plot(
    :body_size,
    :biomass,
    group = :org_type,
    seriestype = :scatter,
    xlabel = "body size",
    ylabel = "biomass",
    #title = "Richness = bottom-heavy, biomass = inverted",
    legend = false
)
# put all plots together
all_shapes_plots = Plots.plot(
    topheavy_pryamid_plot,
    topheavy_uniform_plot,
    topheavy_inverted_plot,
    uniform_pyramid_plot,
    uniform_uniform_plot,
    uniform_inverted_plot,
    bottomheavy_pryamid_plot,
    bottomheavy_uniform_plot,
    bottomheavy_inverted_plot,
    layout = (3,3),
    #legend = false,
    #layout = Plots.grid(3,3, heights = [0.3, 0.3, 0.3]),
    #xlabel = "body size",
    #ylabel = "biomass"
    title = ["1", "2", "3", "4", "5", "6", "7", "8", "9"]
)
Plots.savefig(all_shapes_plots, "./figs/all_shapes_plots.png")

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
abstract InterpType
type InterpNearest <: InterpType; end
type InterpLinear <: InterpType; end
type InterpQuadratic <: InterpType; end
#type InterpCubic <: InterpType; end   # TODO


text:
module ANOVA
using GLM: ccdf, cholesky!, CoefTable, deviance, dof_residual, FDist, fit,
           ftest, LinearModel, nobs

import Base: show, summary

"""
    Anova(mod, anovatype = 3)

Compute an ANOVA table for the given linear model.
"""
struct Anova
  Source::Vector{String}
  DF::Vector{Float64}
  SS::Vector{Float64}
  MSS::Vector{Float64}
  F::Vector{Float64}
  p::Vector{Float64}
  function Anova(mod; anovatype = 3)
    mf = mod.mf
    mm = mod.mm
    mod = mod.model

    eff = effects(mod) # get effects

    response = mf.terms.eterms[1]
    terms = mf.terms.terms

    ## calculate variables for residuals (this is the full modell) ##
    DFres = dof_residual(mod)           # degrees of freedom
    SSres = deviance(mod)               # calculate sum of squares of weighted residuals
    MSSres = SSres / DFres              # calculate mean sum of squares for residuals

    # check if there is an intercept
    v = mf.terms.intercept

    # get assigned numbers for factors
    factors_assig = unique(mm.assign)

    ## calculate some variables for computation ##
    n = nobs(mod)                                # sample size
    Nfactors = count(x -> x > 0, factors_assig)  # number of factors

    ## create arrays for ANOVA factors ##
    DF = zeros(Nfactors)   # will contain degrees of freedoms for factor
    SS = zeros(Nfactors)   # will contain sum of squares for factor
    MSS = zeros(Nfactors)  # will contain mean sum of squares for factor
    RSS = zeros(Nfactors)  #
    F = zeros(Nfactors)    # will contain F values for factors
    p = zeros(Nfactors)    # will contain pvalue for factors

    for i in 1:Nfactors
      mask = mm.assign .== factors_assig[i+v]
      if anovatype == 3 # if this is a type III ANOVA, we use a similar procedure to R's drop1() function
        fullmod = mod
        newmod = droptermbymask(fullmod, mask)
        RSS[i] = deviance(newmod)
        SS[i] = RSS[i] - deviance(fullmod)
        DF[i] = count(mask)
        MSS[i] = SS[i]/DF[i]
        F[i] = ftest(fullmod,newmod).fstat[1]
      elseif anovatype == 2 # if this is a type II ANOVA
        full_model_mask = falses(length(mm.assign))
        for j in 1:Nfactors
          term = terms[j]
          if (isa(term, Expr) && isa(terms[i], Symbol)) # this is not an interaction term
            if terms[i] ∈ term.args
              full_model_mask = merge_bool_array(full_model_mask, mm.assign .== factors_assig[j+v])
            end
          elseif (isa(term, Expr) && isa(terms[i], Expr)) # this is an interaction term
            if (all(x -> x ∈ term.args, terms[i].args) && term ≠ terms[i])
              full_model_mask = merge_bool_array(full_model_mask, mm.assign .== factors_assig[j+v])
            end
          end
        end
        fullmod = droptermbymask(mod, full_model_mask) # we have to drop all interaction terms containing the current term in the full modell
        newmod = droptermbymask(mod, merge_bool_array(full_model_mask, mask))
        RSS[i] = deviance(newmod)
        SS[i] = RSS[i] - deviance(fullmod)
        DF[i] = count(mask)
        MSS[i] = SS[i] / DF[i]
        F[i] = MSS[i] / MSSres
      else # type I ANOVA
        SS[i] = sum(abs2, eff[mask])
        RSS[i] = 0
        DF[i] = count(mask)
        MSS[i] = SS[i] / DF[i]
        F[i] = MSS[i] / MSSres
      end
      p[i] = ccdf(FDist(DF[i], DFres), F[i])
    end
    new(vcat(string.(terms), "Residuals"),
        vcat(DF, DFres),
        vcat(SS, SSres),
        vcat(MSS, MSSres),
        vcat(F, 0),
        vcat(p, 0))
  end
  Anova(Source::AbstractVector{<:AbstractString},
        DF::AbstractVector{<:Real},
        SS::AbstractVector{<:Real},
        F::AbstractVector{<:Real},
        p::AbstractVector{<:Real}) = new(Source, DF, SS, SS ./ DF, F, p)
end

## helper functions ##

merge_bool_array(a, b) = ifelse.(b, b, a)

# calculate effects for type I ANOVA
effects(mod) = (mod.pp.X / cholesky!(mod.pp).U)' * mod.rr.y

# this will drop variables from the matrix using the mask argument then fit a new linear model
droptermbymask(mod, mask) = fit(LinearModel, mod.pp.X[:,.!mask], mod.rr.y)

summary(io::IO, obj::Anova) =
  println(io, "ANOVA table for linear model")

show(io::IO, obj::Anova) =
  println(io,
          summary(obj),
          CoefTable(hcat(obj.DF, obj.SS, obj.MSS, obj.F, obj.p),
                    ["DF", "SS", "MSS", "F", "p"],
                    obj.Source,
                    5))

export Anova

end # module


text:
cla()
x = range(0, stop=π, length=250)
for α ∈ [1, 4, 7, 10]
    y = @. exp(-√α*x)*sin(2*x/α)*α^2
    plot!(x, y, ls="-", lw=0.05)
end
mx = x[1:10:end]
scatter!(mx, sin.(mx), marker="star", col="salmon", msize=0.3)
xlim(0, π)


text:
gpu(A::Operator, GPUType) = Operator(GPUType(data(A)),dims(A))

gpu(A::Ket, GPUType) = Ket(GPUType(data(A)),dims(A))

gpu(L::Liouvillian, GPUType) = Liouvillian(L.dims, GPUType(SparseMatrixCSC{ComplexF32}(L.L₀)), GPUType.(SparseMatrixCSC{ComplexF32}.(L.Lₙ)), L.fₙ, L.pₙ)


text:
# workaround is to load ImageMagick.jl first, front and center
# see https://github.com/JuliaIO/ImageMagick.jl/issues/142
# and https://github.com/JuliaIO/ImageMagick.jl/issues/130
using ImageMagick

using Pkg
# for now add GraffSDK all the time
Pkg.add(PackageSpec(url = "https://github.com/GearsAD/GraffSDK.jl.git"))

using GraffSDK
using Arena
using Test

include("testBasicAnimations.jl")


text:
# Autogenerated wrapper script for kubectl_jll for x86_64-linux-gnu
export kubectl

JLLWrappers.@generate_wrapper_header("kubectl")
JLLWrappers.@declare_executable_product(kubectl)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        kubectl,
        "bin/kubectl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# If adding models from MLJModels for testing purposes, then do the
# following in the interface file (eg, DecisionTree.jl):

# - change `import ..DecisionTree` to `import DecisionTree`
# - remove wrapping as module

module Models

using MLJBase
import MLJModelInterface: @mlj_model, metadata_model, metadata_pkg
import MLJModelInterface

include("models/Constant.jl")
include("models/DecisionTree.jl")
include("models/NearestNeighbors.jl")
include("models/MultivariateStats.jl")
include("models/Transformers.jl")
include("models/foobarmodel.jl")
include("models/simple_composite_model.jl")
include("models/ensembles.jl")

end


text:
module MDTools

using Zygote
using Unitful
using UnitfulAtomic
using Distances
using LoopVectorization

import Base.*
import Base.+

include("integrators.jl")
include("atoms.jl")
include("potentials.jl")
include("system.jl")

export AbstractMDIntegrator,
       add!,
       CubicBox,
       energy,
       force,
       LeapFrog,
       propagate!

function force!(f, potential, coordinates)
    f .= - gradient(potential,coordinates)[1]
end

force(potential, coordinates) = - gradient(potential, coordinates)[1]

function energy(coordinates, velocity, potential, mass)
    Ek = 0.5*sum( mass.*velocity.^2 )
    Ep = potential(coordinates)
    return Ek, Ep, Ek+Ep
end

greet() = print("Hello World!")

end # module


text:
## NEURAL NETWORK CLASSIFIER

seed!(1234)
N = 300
X = MLJBase.table(rand(Float32, N, 4));
ycont = 2*X.x1 - X.x3 + 0.1*rand(N)
m, M = minimum(ycont), maximum(ycont)
_, a, b, _ = range(m, stop=M, length=4) |> collect
y = map(ycont) do η
    if η < 0.9*a
        'a'
    elseif η < 1.1*b
        'b'
    else
        'c'
    end
end |> categorical;

# TODO: replace Short2 -> Short when
# https://github.com/FluxML/Flux.jl/issues/1372 is resolved:
builder = Short2()
optimiser = Flux.Optimise.ADAM(0.03)

losses = []

@testset_accelerated "NeuralNetworkClassifier" accel begin
    Random.seed!(123)
    basictest(MLJFlux.NeuralNetworkClassifier,
              X,
              y,
              builder,
              optimiser,
              0.85,
              accel)

    train, test = MLJBase.partition(1:N, 0.7)

    # baseline loss (predict constant probability distribution):
    dict = StatsBase.countmap(y[train])
    prob_given_class = Dict{CategoricalArrays.CategoricalValue,Float64}()
    for (k, v) in dict
        prob_given_class[k] = dict[k]/length(train)
    end
    dist = MLJBase.UnivariateFinite(prob_given_class)
    loss_baseline =
        MLJBase.cross_entropy(fill(dist, length(test)), y[test]) |> mean

    # check flux model is an improvement on predicting constant
    # distribution:
    stable_rng = StableRNGs.StableRNG(123)
    model = MLJFlux.NeuralNetworkClassifier(epochs=50,
                                            builder=builder,
                                            optimiser=optimiser,
                                            acceleration=accel,
                                            batch_size=10,
                                            rng=stable_rng)
    @time mach = fit!(machine(model, X, y), rows=train, verbosity=0)
    first_last_training_loss = MLJBase.report(mach)[1][[1, end]]
    push!(losses, first_last_training_loss[2])
    yhat = MLJBase.predict(mach, rows=test);
    @test mean(MLJBase.cross_entropy(yhat, y[test])) < 0.95*loss_baseline

    optimisertest(MLJFlux.NeuralNetworkClassifier,
                  X,
                  y,
                  builder,
                  optimiser,
                  accel)

end

# check different resources (CPU1, CUDALibs, etc)) give about the same loss:
reference = losses[1]
@test all(x->abs(x - reference)/reference < 1e-5, losses[2:end])

true


text:
#=

Implements the empirical cumulative distribution function given an array
of observations.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-08-25

=#

"""
One-dimensional empirical distribution function given a vector of
observations.

##### Fields

- `observations::Vector`: The vector of observations
"""
mutable struct ECDF
    observations::Vector
end

"""
Evaluate the empirical cdf at one or more points

##### Arguments

- `x::Union{Real, Array}`: The point(s) at which to evaluate the ECDF
"""
(e::ECDF)(x::Real) = mean(e.observations .<= x)
(e::ECDF)(x::AbstractArray) = e.(x)

@deprecate ecdf(e::ECDF, x) e(x)


text:
struct Credential
    access_key:: String
    secret_key:: String
end
export Credential

module Methods
    @enum Method GET POST DELETE
    export Method
end
export Methods


text:
#####
##### Best n-card hand rank, from combinations.
#####

function min_hand_rank_full(eval1, cards)
    hand_rank_2 = evaluate5(cards)
    if first(eval1) <= hand_rank_2
        return eval1
    else
        hand_type_2 = hand_type_binary_search(hand_rank_2)
        return (hand_rank_2, hand_type_2, cards)
    end
end

function min_hand_rank_compact(hand_rank_1, cards)
    hand_rank_2 = evaluate5(cards)
    if hand_rank_1 < hand_rank_2
        return hand_rank_1
    else
        return hand_rank_2
    end
end

function best_hand_rank_from_7_cards_full(c)
    cards=(c[1],c[2],c[3],c[4],c[5]) # first combination

    hand_rank = evaluate5(cards)
    hand_type = hand_type_binary_search(hand_rank)
    hr = (hand_rank, hand_type, cards)

    cards=(c[1],c[2],c[3],c[4],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[3],c[4],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[3],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[3],c[5],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[3],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[4],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[4],c[5],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[4],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[5],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[3],c[4],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[3],c[4],c[5],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[3],c[4],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[3],c[5],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[4],c[5],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[2],c[3],c[4],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[2],c[3],c[4],c[5],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[2],c[3],c[4],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[2],c[3],c[5],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[2],c[4],c[5],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    cards=(c[3],c[4],c[5],c[6],c[7]); hr = min_hand_rank_full(hr, cards)
    return hr
end

function best_hand_rank_from_6_cards_full(c)
    cards=(c[1],c[2],c[3],c[4],c[5]); # first combination

    hand_rank = evaluate5(cards);
    hand_type = hand_type_binary_search(hand_rank);
    hr = (hand_rank, hand_type, cards)

    cards=(c[1],c[2],c[3],c[4],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[3],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[2],c[4],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[1],c[3],c[4],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    cards=(c[2],c[3],c[4],c[5],c[6]); hr = min_hand_rank_full(hr, cards)
    return hr
end

function best_hand_rank_from_7_cards_compact(c)
    cards=(c[1],c[2],c[3],c[4],c[5]); hr = evaluate5(cards)
    cards=(c[1],c[2],c[3],c[4],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[3],c[4],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[3],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[3],c[5],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[3],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[4],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[4],c[5],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[4],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[2],c[5],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[3],c[4],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[3],c[4],c[5],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[3],c[4],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[3],c[5],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[1],c[4],c[5],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[2],c[3],c[4],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[2],c[3],c[4],c[5],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[2],c[3],c[4],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[2],c[3],c[5],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[2],c[4],c[5],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    cards=(c[3],c[4],c[5],c[6],c[7]); hr = min_hand_rank_compact(hr, cards)
    return hr
end

function best_hand_rank_from_6_cards_compact(c)
    cards = (c[1],c[2],c[3],c[4],c[5]); hr = evaluate5(cards)
    cards = (c[1],c[2],c[3],c[4],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards = (c[1],c[2],c[3],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards = (c[1],c[2],c[4],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards = (c[1],c[3],c[4],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    cards = (c[2],c[3],c[4],c[5],c[6]); hr = min_hand_rank_compact(hr, cards)
    return hr
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 193ee450-7b32-4394-9401-7b91bce21f31
using LaTeXStrings, DifferentialEquations, LinearAlgebra, Plots, SparseArrays

# ╔═╡ 385b2bda-24d8-464c-9257-c4393a0b1c29
using ForwardDiff

# ╔═╡ 13137d52-2ede-4a63-8eac-43b13def7466
md"# Self Activating Toggle Switch"

# ╔═╡ cf20c08b-5226-4c6c-8741-6857ede519a2
md"_Add required function to implement model_"

# ╔═╡ 1dc242f8-b0a1-11eb-2585-a121968b7921
begin
	mma(x) = 1/(1+x)
	mmr(x) = 1/(1+1/x)
	hfr(n) = x -> 1/(1 + x^n)
	hfa(n) = x -> 1/(1 + (1/x)^n)
	hs(λ, n) = x-> hfr(n)(x) + λ * hfa(n)(x)
	hs(λ, n, x₀) = x-> hs(λ, n)(x/x₀)
end

# ╔═╡ d7ef6de7-145c-4478-aa59-011562a5b8de
md"_Implement model piecewise (separate production and degradation terms) for simulations and computation_"

# ╔═╡ 812b4a91-8c80-4499-abd4-e65a464db31b
begin
	n = 4 # hill coefficient
	λ₁ = 3 # for bistable
	λ₂ = 10 # for tristable
	λₙ = 0.1 # repressor
	x₁ = 80 # cutoff count bistable
	x₂ = 25 # cutoff count tristable
	xₙ = 20 # cuttoff count repress
	p = 5 # unregulated production rate
	d = 0.1 # death rate
	proda(a,b) = p*hs(λ₂, n, x₂)(a)*hs(λₙ, n, xₙ)(b)
	dega(a, b) = d*a
	dadt(a, b) = proda(a, b) - dega(a, b)
	prodb(a, b) = proda(b, a)
	degb(a, b)= d*b
	dbdt(a, b) = prodb(a,b) - degb(a, b)
	dfx = dadt
	dfy = dbdt
end

# ╔═╡ 4d905902-e730-4b79-a516-94ec9e236d59
begin
	fx(x::Vector) = dfx(x[1], x[2])
    fy(x::Vector) = dfy(x[1], x[2])
    f(x::Vector) = [fx(x), fy(x)]
end

# ╔═╡ 8416fdf7-1d4a-4af7-ba12-30f609c5dcdf
md"_Run Helmholtz decomposition using ForwardDiff for the divergence_"

# ╔═╡ 3ffd11ce-5388-4dfa-81ec-535c633ff378
begin
	# savefig(p_pot, "ch4_comp_pot.pdf")
	# savefig(p_rhs, "ch4_comp_rhs.pdf")
end

# ╔═╡ c105a211-66fb-495d-83bc-8ecf5679c37a
md"_Compute Vector Field_"

# ╔═╡ 006d5d3c-1b9e-442e-8fca-d7aa3a3abd35
begin
	n_min = 0
	n_max = 550
end

# ╔═╡ da1ae333-aafb-4ca7-a332-5d17a287c3d9
begin
	dq = (n_max - n_min) / 10.
	xs = n_min:dq:n_max+1
	ys = n_min:dq:n_max+1

	dxs = n_min:.5:n_max+1
	dys = n_min:.5:n_max+1

	df(x, y) = normalize([dfx(x, y), dfy(x,y)]) .* 30

	xxs = [x for x in xs for y in ys]
	yys = [y for x in xs for y in ys]

	quiver(xxs, yys, quiver=df)
	contour!(dxs, dys, dfx, levels=[0], color=:orange, label="X Nullcline", colorbar = false)
	contour!(dxs, dys, dfy, levels=[0], color=:green, label = "Y Nullcline", colorbar = false, xlabel=L"X", ylabel=L"Y", legend=:topright,)
end

# ╔═╡ ee260481-d1ae-4a38-b922-6e10bbfe945e
begin
	random_count() = floor(Int, n_min + rand() * (n_max - n_min))
	rand_u0() = [random_count(), random_count()]
	u0 = rand_u0()
	u0 = [5, 500]
	t = (0.0, 30000.0)
end

# ╔═╡ c077f2a7-6dd8-41a6-88ce-56ee38707081
md"_Stochastic DiffEq of Langevin equation for this system_

Note: `abs` is not supposed to be there but the simulation fails without it."

# ╔═╡ 32a77522-b389-4435-b92f-d05449f8b8d1
begin
	function f(du, u, p, t)
		du[1] = dfx(u...)
		du[2] = dfy(u...)
		du[1], du[2]
	end

	g(du, u, p, t) = begin
	  du[1,1] = sqrt(abs(proda(u[1], u[2])))
	  du[1,2] = -sqrt(abs(degb(u[1], u[2])))
	  du[1,3] = 0
	  du[1,4] = 0
	  du[2,1] = 0
	  du[2,2] = 0
	  du[2,3] = sqrt(abs(prodb(u[1], u[2])))
	  du[2,4] = -sqrt(abs(degb(u[1], u[2])))
	end
end

# ╔═╡ 1cd42539-b8db-4e2b-91a3-2038f299cf43
begin

	minval = 1
	maxval =  505
	nsample = 256
	rhs_xs = range(minval, maxval, length=nsample)
	rhs_ys = range(minval, maxval, length=nsample)

	divergence(F::Function, pt) = sum(
    	diag(
        	ForwardDiff.jacobian(F, pt)
        	)
    	)

	rhs_new = [divergence(f, [x, y]) for x=rhs_xs, y=rhs_ys]
end

# ╔═╡ bfb2c556-d3ec-449e-9f09-30d253bf3938
begin
	p_rhs = heatmap(rhs_xs, rhs_ys, rhs_new, xlabel=L"X", ylabel=L"Y", colorbar=:none)
end

# ╔═╡ 5989e134-f90b-43f0-bb2f-b2ffc5c74cb8
begin
	sde_prob_ssa= SDEProblem(f, g, u0, t, noise_rate_prototype=zeros(2,4))
	sde_sol = solve(sde_prob_ssa)
end

# ╔═╡ 46b15beb-6ee0-405e-96ca-6c0330a3c461
plot(sde_sol)

# ╔═╡ 1e7f2f9f-b8e6-44f8-8259-e2423368679d
md"_Gillespie Jump Process formulation for the same problem_"

# ╔═╡ f173fe40-ee7d-4513-8511-a18ec553ef5f
begin
	proda!(integrator) = integrator.u[1] += 1
	prodb!(integrator) = integrator.u[2] += 1
	dega!(integrator) = integrator.u[1] -= 1
	degb!(integrator) = integrator.u[2] -= 1
	jumps = Dict(proda => proda!, prodb => prodb!, dega => dega!, degb => degb!)
	constJumps = []
	for (rate, jump) in jumps
		push!(constJumps, ConstantRateJump((u, p, t) -> rate(u...), jump))
	end
	constJumps
end

# ╔═╡ 5c2c49e7-1948-41bb-82cc-f6a2460956fc
begin
	dprob_ssa = DiscreteProblem(u0, t)
	jprob_ssa = JumpProblem(dprob_ssa, Direct(), constJumps...)
end

# ╔═╡ 83b082e0-643e-4606-9db9-2957f3c6a470
ssa_sol = solve(jprob_ssa, SSAStepper())

# ╔═╡ 68c0fa0f-ff2d-4208-9a7d-415d0fa0abd1
stoch_tri = plot(ssa_sol, labels=[L"X" L"Y"], xlabel=L"t", legend=(0.9, 0.5))

# ╔═╡ 0715324c-8698-4a22-a146-9bc4c29c3a63
stoch_tri3 = plot(ssa_sol, vars=[(1, 2)], xlabel=L"X", ylabel=L"Y",legend=:none)

# ╔═╡ 73482204-dc5d-4fd5-aabb-72f284ca4416
pstoch = plot(stoch_tri, stoch_tri3, layout=(1, 2), size=(900, 500), title = ["A" "B"], titlelocation = :left)

# ╔═╡ 98664d78-abc5-45ca-a750-6184e0c285a8
# savefig(pstoch, "ch4_tri_stoch.pdf")

# ╔═╡ 2545abfb-df64-4a56-b217-d11fb7f2e95f
md"_Setup helper functions for running longer simulations_"

# ╔═╡ 0276796e-4f61-4d86-95b3-11df1b0ae242
function get_u0(i)
	offset = 1.0
	new_u0 = [0.0, 0.0]
	if i % 3 == 0
		new_u0 = [55, 55]
	elseif i % 3 == 1
		new_u0 .+= [5, 500]
	elseif i % 3 == 2
		new_u0 .+= [500, 5]
	end
	new_u0
end

# ╔═╡ a5d2a26e-5d33-480a-8712-8a38b43de1c9
md"_Storage for unnormalized probabilities that sum up the `dt`s in each bin_"

# ╔═╡ 31ed699a-77e2-4388-88a1-1dd740d56e40
md"_Rerun the following cell to add more values to the SSA simulation_"

# ╔═╡ 5250b601-7647-41cc-861a-393ad916cace
md"_Rerun the following cell to have the SDE simulation catch up to the SSA_"

# ╔═╡ d42b00d6-6f3f-44f9-8e5a-0ed8c8cc25ce
md"_the SDE probability distribution is complete non-sense_"

# ╔═╡ 076d3e3f-e3cd-4b60-8440-61c497b186a0
# savefig(ps, "ch4_ad_hoc.pdf")

# ╔═╡ a0bbd0f6-24d8-4705-b142-c1544346161c
md"## Appendix"

# ╔═╡ 125f4bfb-04b4-4fb3-8f44-413e8d8ee66c
md"_Helper function to generate -∇ ² to solve the Helmhotz equation_"

# ╔═╡ 13f7e53e-d2c9-49ec-b5c0-64e32299e758
begin
	eye(N, M) = sparse(I, N, M)
	eye(M) = eye(M, M)
	diff1(M) = [ [1.0 zeros(1, M - 1)]; diagm(1 => ones(M - 1)) - eye(M) ]
	sdiff1(M) = sparse(diff1(M))

	# make the discrete -Laplacian in 2d, with Dirichlet boundaries
	function laplacian2d(Nx, Ny, Lx, Ly)
		dx = Lx / (Nx + 1)
		dy = Ly / (Ny + 1)
		Dx = sdiff1(Nx) / dx
		Dy = sdiff1(Ny) / dy
		Ax = Dx' * Dx
		Ay = Dy' * Dy
		return Dx, Dy, kron(eye(Ny), Ax) + kron(Ay, eye(Nx))
	end
end

# ╔═╡ bcc69011-ae43-452e-b5bd-5c5eecaaae19
begin
	_, _, Δ = laplacian2d(nsample, nsample, maxval - minval, maxval - minval)
	sol = reshape(Δ \ reshape(rhs_new, length(rhs_new)),  size(rhs_new)...)
end

# ╔═╡ fdac6046-4801-434a-9d03-01e4c2e9403c
p_pot = surface(rhs_xs, rhs_ys, sol, xlabel=L"X", ylabel=L"Y", colorbar=:none)

# ╔═╡ 60be1d36-ebfa-4ab7-ba68-5c35b6d41aa0
md"_Helper function to compute 2D histograms for Probability distributions_"

# ╔═╡ f245c60d-84f8-4f24-b2a0-c55abd991641
begin
	nbins = 701
	bin_min = 0
	bin_max = 700
	xbins = range(bin_min, bin_max, length=nbins)
	ybins = range(bin_min, bin_max, length=nbins)

	function bucket_idx(val, vec)
		idx = 1
		while val > vec[idx]
			idx += 1
		end
		idx
	end

	function althist(x, y, dt, xedges, yedges)
	   counts = spzeros(length(yedges)-1, length(xedges)-1)
		for i=1:length(x) - 1
			r = bucket_idx(y[i], yedges)
			c = bucket_idx(x[i], xedges)
			counts[r,c] += dt[i]
		end
 		counts
	end
end

# ╔═╡ 18b7a92a-57f7-4b43-9d59-65b24889b676
begin
	unnorm_ssa = spzeros(length(xbins)-1, length(ybins)-1)
	unnorm_sde = spzeros(length(xbins)-1, length(ybins)-1)
	tspan = (0., 5e4)
end

# ╔═╡ 049020d7-1d6f-4603-9827-c1b1f22f7de6
prob_ssa = unnorm_ssa ./ sum(unnorm_ssa)

# ╔═╡ 6800f1bc-793e-49df-a60a-57305af8251a
prob_sde = unnorm_sde ./ sum(unnorm_sde)

# ╔═╡ f6a2e316-462e-4d3a-a6b2-37f02896599e
begin
	plotly()
	pl = surface(prob_sde, xlabel="X", ylabel="Y")
	gr()
	pl
end

# ╔═╡ 862c3e0a-261a-4e31-b986-f44d1056f66f
begin
	global n_points_sim
	for i in 1:3*10
		long_jprob_ssa = JumpProblem(
			DiscreteProblem(get_u0(i), tspan), Direct(), constJumps...)
		long_ssa_sol = solve(long_jprob_ssa, SSAStepper())
		unnorm_ssa .+= althist(long_ssa_sol[1,:], long_ssa_sol[2, :], diff(long_ssa_sol.t[:]), xbins, ybins)
	end
	unnorm_ssa
end

# ╔═╡ fe812819-697c-4840-a5df-f8ea1dd5e59b
begin
	i = 45
	num_sde_run_failures = 44
	while sum(unnorm_sde) < sum(unnorm_ssa)
		global i
		global num_sde_run_failures

		sde_run = SDEProblem(f, g, get_u0(i), tspan, noise_rate_prototype=zeros(2,4))
		sol = solve(sde_prob_ssa)

		if sol.t[end] < tspan[end]
			# the simulation did not complete succesfully.
			# This is a sign of divergence and a failing approximation
			num_sde_run_failures += 1
		end
		unnorm_sde .+= althist(sol[1,:], sol[2, :], diff(sol.t[:]), xbins, ybins)
		i += 1
	end
	unnorm_sde, num_sde_run_failures, i
end

# ╔═╡ d9665278-3e2d-47b9-b955-1ddee1fe61e9
begin
	xplot = range(bin_min, bin_max, length=nbins-1)
	yplot = range(bin_min, bin_max, length=nbins-1)
end

# ╔═╡ a26b07ac-5b4d-40c0-b84c-7b86c19a5362
 begin
	pp1 = surface(xplot, yplot,-log.(prob_ssa), c = :thermal, title="ad hoc Quasi-potential")
	pp2 =  surface(xplot, yplot, prob_ssa, c = :grayC, title="Prob_ssaability Distribution", zaxis=(formatter=y->string(round(Int, y / 10^-4))*"e-4" ) )
	ps = plot(pp2, pp1, layout =(1, 2), size=(800, 600), xlabel=L"X", ylabel=L"Y", colorbar=:none)
end

# ╔═╡ Cell order:
# ╟─13137d52-2ede-4a63-8eac-43b13def7466
# ╠═193ee450-7b32-4394-9401-7b91bce21f31
# ╠═385b2bda-24d8-464c-9257-c4393a0b1c29
# ╟─cf20c08b-5226-4c6c-8741-6857ede519a2
# ╠═1dc242f8-b0a1-11eb-2585-a121968b7921
# ╟─d7ef6de7-145c-4478-aa59-011562a5b8de
# ╠═812b4a91-8c80-4499-abd4-e65a464db31b
# ╠═4d905902-e730-4b79-a516-94ec9e236d59
# ╟─8416fdf7-1d4a-4af7-ba12-30f609c5dcdf
# ╠═1cd42539-b8db-4e2b-91a3-2038f299cf43
# ╠═bfb2c556-d3ec-449e-9f09-30d253bf3938
# ╠═bcc69011-ae43-452e-b5bd-5c5eecaaae19
# ╠═fdac6046-4801-434a-9d03-01e4c2e9403c
# ╠═3ffd11ce-5388-4dfa-81ec-535c633ff378
# ╟─c105a211-66fb-495d-83bc-8ecf5679c37a
# ╠═006d5d3c-1b9e-442e-8fca-d7aa3a3abd35
# ╠═da1ae333-aafb-4ca7-a332-5d17a287c3d9
# ╠═ee260481-d1ae-4a38-b922-6e10bbfe945e
# ╟─c077f2a7-6dd8-41a6-88ce-56ee38707081
# ╠═32a77522-b389-4435-b92f-d05449f8b8d1
# ╠═5989e134-f90b-43f0-bb2f-b2ffc5c74cb8
# ╠═46b15beb-6ee0-405e-96ca-6c0330a3c461
# ╟─1e7f2f9f-b8e6-44f8-8259-e2423368679d
# ╠═f173fe40-ee7d-4513-8511-a18ec553ef5f
# ╠═5c2c49e7-1948-41bb-82cc-f6a2460956fc
# ╠═83b082e0-643e-4606-9db9-2957f3c6a470
# ╠═68c0fa0f-ff2d-4208-9a7d-415d0fa0abd1
# ╠═0715324c-8698-4a22-a146-9bc4c29c3a63
# ╠═73482204-dc5d-4fd5-aabb-72f284ca4416
# ╠═98664d78-abc5-45ca-a750-6184e0c285a8
# ╟─2545abfb-df64-4a56-b217-d11fb7f2e95f
# ╠═0276796e-4f61-4d86-95b3-11df1b0ae242
# ╟─a5d2a26e-5d33-480a-8712-8a38b43de1c9
# ╠═18b7a92a-57f7-4b43-9d59-65b24889b676
# ╟─31ed699a-77e2-4388-88a1-1dd740d56e40
# ╠═862c3e0a-261a-4e31-b986-f44d1056f66f
# ╟─5250b601-7647-41cc-861a-393ad916cace
# ╠═fe812819-697c-4840-a5df-f8ea1dd5e59b
# ╠═049020d7-1d6f-4603-9827-c1b1f22f7de6
# ╠═6800f1bc-793e-49df-a60a-57305af8251a
# ╟─d42b00d6-6f3f-44f9-8e5a-0ed8c8cc25ce
# ╠═f6a2e316-462e-4d3a-a6b2-37f02896599e
# ╠═d9665278-3e2d-47b9-b955-1ddee1fe61e9
# ╠═a26b07ac-5b4d-40c0-b84c-7b86c19a5362
# ╠═076d3e3f-e3cd-4b60-8440-61c497b186a0
# ╟─a0bbd0f6-24d8-4705-b142-c1544346161c
# ╟─125f4bfb-04b4-4fb3-8f44-413e8d8ee66c
# ╠═13f7e53e-d2c9-49ec-b5c0-64e32299e758
# ╟─60be1d36-ebfa-4ab7-ba68-5c35b6d41aa0
# ╠═f245c60d-84f8-4f24-b2a0-c55abd991641


text:
module Tracer

using Cassette
import Cassette: overdub, prehook

Cassette.@context ArrayTraceCtx

_track(ctx::ArrayTraceCtx) = nothing
_track(ctx::ArrayTraceCtx, x, args...) = _track(ctx, args...)
function _track(ctx::ArrayTraceCtx, x::Array, args...)
    get!(ctx.metadata, objectid(x), length(x))
    _track(ctx, args...)
end

# Prehook everythig, grab all references to arrays
prehook(ctx::ArrayTraceCtx, f, args...) = _track(ctx, args...)

function track(f, args...)
    metadata = IdDict{UInt64, Int64}()
    ctx = ArrayTraceCtx(metadata = metadata)
    overdub(ctx, f, args...)
    return ctx.metadata
end

end


text:
using Missings

x1 = [1.0, 2.0, 3.0, 4.0, 5.0, missing]
ans1 = [-0.6885895075239652, -0.3442947537619826, 0.0, 0.3442947537619826, 0.6885895075239652, 0.0]

res, nNAentries, NAmed, zeroMAD = prepareColBicor(x1, 1.0, 1, false)
@test res == ans1
@test nNAentries == 1
@test NAmed == false
@test zeroMAD == false

# x2 = [0.0, 0.0, 0.0, 0.0, 0.0, missing]
#
# @test_throws  prepareColBicor(x2, 1.0, 1, false)

# x3 = [missing, missing, missing, missing, missing, missing]
# ans3 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
#
# @time res, nNAentries, NAmed, zeroMAD = prepareColBicor(x3, 1.0, 1, false)
# @test res == ans3
# @test nNAentries == 1
# @test NAmed == true
# @test zeroMAD == true


text:
"""
    ExactReach

ExactReach performs exact reachability analysis to compute the output reachable set for a network.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: HPolytope
3. Output: AbstractPolytope

# Return
`ReachabilityResult`

# Method
Exact reachability analysis.

# Property
Sound and complete.

# Reference
[W. Xiang, H.-D. Tran, and T. T. Johnson,
"Reachable Set Computation and Safety Verification for Neural Networks with ReLU Activations,"
*ArXiv Preprint ArXiv:1712.08163*, 2017.](https://arxiv.org/abs/1712.08163)
"""
struct ExactReach <: Solver end

function solve(solver::ExactReach, problem::Problem)
    reach = forward_network(solver, problem.network, problem.input)
    return check_inclusion(reach, problem.output)
end

# Forward_linear that coverts the input to HPolytope
forward_linear(solver::ExactReach, L::Layer, input) = forward_linear(solver, L, convert(HPolytope, input))

# Forward linear for a single polytope and a list of polytopes
forward_linear(solver::ExactReach, L::Layer, input::HPolytope) = affine_map(L, input)
forward_linear(solver::ExactReach, L::Layer, input::Vector{<:HPolytope}) = [affine_map(L, set) for set in input]

# Forward act for a single polytope and a list of polytopes
forward_act(solver::ExactReach, L::Layer, input::HPolytope) = forward_partition(L.activation, input)
function forward_act(solver::ExactReach, L::Layer, input::Vector{<:HPolytope})
    output = Vector{HPolytope}(undef, 0)
    for i in 1:length(input)
        append!(output, forward_partition(L.activation, input[i]))
    end
    return output
end


text:
module Pricing

using HTTP, JSON3, Dates

export getPrice, checkMarket, streamPrice, getLatestCandles

# ------------------------------------------------------------------------------------
# /accounts/{accountID}/pricing Endpoint
# ------------------------------------------------------------------------------------

"""
Ask / Bid pricing data

# Arguments
- 'price': Price of the ask / bid
- 'liquidity': Liquidity of the ask / bid
"""
mutable struct priceBucket
    price # Price of the ask / bid
    liquidity # liquidity of the ask / bid

    priceBucket() = new()
end

"""
Pricing data of an instrument

# Arguments
- 'type': Type
- 'instrument': Instrument of the price
- 'time': Time of the price update
- 'bids::Vector{priceBucket}': bid info
- 'asks::Vector{priceBucket}': ask info
- 'closeoutBid': Closeout bid price
- 'closeoutAsk': Closeout ask price
- 'tradeable': Can you trade this instrument
"""
mutable struct price
    type # Type
    instrument
    time # Time of the price update
    bids::Vector{priceBucket} # Bid information
    asks::Vector{priceBucket} # Ask information
    closeoutBid # Closeout bid price
    closeoutAsk # Closeout Ask price
    tradeable # Can you trade this instrument

    price() = new()
end

"Needed for JSON parsing"
mutable struct priceTopLayer
    prices::Vector{price} # Prices
    time

    priceTopLayer() = new()
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{priceBucket}) = JSON3.Mutable()
JSON3.StructType(::Type{price}) = JSON3.Mutable()
JSON3.excludes(::Type{price})=(:status,:quoteHomeConversionFactors,:unitsAvailable) #Ignore deprecated fields
JSON3.StructType(::Type{priceTopLayer}) = JSON3.Mutable()

"Coerce pricing data into proper types"
function coercePrice(price::price)
    # Coerce Asks
    for ask in price.asks
        ask.price = parse(Float32, ask.price)
    end
    # Coerce Bids
    for bid in price.bids
        bid.price = parse(Float32, bid.price)
    end
    price.time = DateTime(first(price.time, 23), Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sssssssssZ"))
    price.closeoutBid = parse(Float32, price.closeoutBid)
    price.closeoutAsk = parse(Float32, price.closeoutAsk)

    return price
end

"""
    function getPrice(config, instruments)

Get the most recent price update of an instrument

Returns an object of type 'priceData'

# Arguments
    - 'config::config': a valid struct with user configuracion data
    - 'instruments::Vector{String}': a vector of valid pairs (e.g. ["EUR_USD","EUR_JPY"])

#Example

    getPrice(userconfig, ["EUR_USD","EUR_JPY"])

"""
function getPrice(config, instruments::Vector{String})

    r = HTTP.get(string("https://", config.hostname, "/v3/accounts/",config.account,"/pricing?instruments=", join(instruments,",")),
                ["Authorization" => string("Bearer ", config.token)])

    data = JSON3.read(r.body, priceTopLayer)

    for priceData in data.prices
        coercePrice(priceData)
    end

    return data.prices #Does not return 'since' datetime=#
end

#TODO: move this to a utils file

"Exception thrown when the market is closed on the weekend"
struct ClosedMarketException <: Exception end

"Get the most recent price update of an instrument"
function checkMarket()
    dt = Dates.now()
    if Dates.dayofweek(dt) >= 5
        if Dates.dayofweek(dt) == 5 & Dates.hour(dt) < 4
        elseif Dates.dayofweek(dt) == 7 & Dates.hour(dt) >= 5
        else
            throw(ClosedMarketException())
        end
    end
end
# ------------------------------------------------------------------------------------
# /accounts/{accountID}/pricing/stream Endpoint
# ------------------------------------------------------------------------------------
"""
    function streamprice(f, config , instruments)

Returns a stream of price objects and apply a function to each one of them

# Arguments
    - 'f::Function': a function to apply to each price struct object. Streamprice accepts do block format
    - 'config::config': a valid struct with user configuracion data
    - 'instruments::Vector{String}': a vector of valid pairs (e.g. ["EUR_USD","EUR_JPY"])

#Example

    streamprice(userconfig, ["EUR_JPY"]) do price
        println(price)
    end
"""
function streamPrice(f::Function, config, instruments::Vector{String})

@async HTTP.open("GET",
        string("https://", config.streamingHostname, "/v3/accounts/",config.account,"/pricing/stream?instruments=", join(instruments,",")),
        ["Authorization" => string("Bearer ", config.token)]) do io
            for line in eachline(io)
                p = JSON3.read(line, price)
                p.type != "HEARTBEAT" && f(p) #Cleans HEARTBEAT ticks -> only prices are sent to f
            end
        end
end

# ------------------------------------------------------------------------------------
# /accounts/{accountID}/candles/latest Endpoint
# ------------------------------------------------------------------------------------

# All structs defined here except latestCandles are common with those in Instrument.jl
"""
Candlestick tick data

# Fields
- 'o': Opening price
- 'h': Highest price
- 'l': Lowest price
- 'c': Closing price
"""
mutable struct candlestickdata
    o   #open
    h   #high
    l   #low
    c   #close

    candlestickdata() = new()
end

"""
Candlesticks

# Fields
- 'time': The time the data corresponds to
- 'bid::candlestickdata': Candlestick tick data for bids
- 'ask::candlestickdata': Candlestick tick data for asks
- 'mid::candlestickdata': Candlestick tick data for the mid
- 'volume': Volume of trades
- 'complete': is this complete data?
"""
mutable struct candlestick
    complete::Bool
    volume
    time
    bid::candlestickdata
    ask::candlestickdata
    mid::candlestickdata

    candlestick() = new()
end

"""
Candles

# Fields
- 'instrument': Instrument the data relates to
- 'granularity': The granularity of the candle data
- 'candles::Vector{candlestick}': The candlestick data
"""
mutable struct candles
    instrument::String
    granularity::String
    candles::Vector{candlestick}

    candles() = new()
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{candlestickdata}) = JSON3.Mutable()
JSON3.StructType(::Type{candlestick}) = JSON3.Mutable()
JSON3.StructType(::Type{candles}) = JSON3.Mutable()

# Conversions to proper Julia types
function coerceCandleStick(config, candle::candlestick)

    RFC = Dates.DateFormat("yyyy-mm-ddTHH:MM:SS.sssssssssZ")

    candle.time = DateTime(candle.time, RFC)
    isdefined(candle, :bid) && (candle.bid = coerceCandleStickData(candle.bid))
    isdefined(candle, :ask) && (candle.ask = coerceCandleStickData(candle.ask))
    isdefined(candle, :mid) && (candle.mid = coerceCandleStickData(candle.mid))

    return candle
end

function coerceCandleStickData(candleData::candlestickdata)

    candleData.o = parse(Float32, candleData.o)
    candleData.h = parse(Float32, candleData.h)
    candleData.l = parse(Float32, candleData.l)
    candleData.c = parse(Float32, candleData.c)

    return candleData
end

struct latestCandles
    latestCandles::Vector{candles}
end

# Declaring JSON3 struct types
JSON3.StructType(::Type{latestCandles}) = JSON3.Struct()

"""

    function getLatestCandles(config, candleSpecs; kwargs...,)

Information includes: time, granularity, open, high, low, close, volume and a complete indicator

# Arguments
- 'config::config': a valid struct with user configuration data
- 'candleSpecs::Vector': A vector of tuples indicating the specifications for candle to retrieve.
    The tuple must have the instrument, granularity and price component in this order

Valid granularities: ["S5","S10","S15","S30","M1","M2","M4","M5","M10","M15","M30","H1","H2","H3","H4","H6","H8","H12","D","W","M"]
Valid prices: "A" for ask, "B" for bid, "M" for medium or a combination of them

# Keyword Arguments
- 'units::Number'
- 'smooth::Bool'
- 'dailyaligment::Int'
- 'alignmentTimezone::String'
- 'weeklyAlignment::String'

# Example

    getLatestCandles(userconfig, [("EUR_USD","M1","AB"),(("EUR_CHF","M5","M"))])
    getLatestCandles(sim02, [("EUR_USD","M1","AB"),(("EUR_CHF","M5","M"))])

"""

function getLatestCandles(config, candleSpecs::Vector, kwargs...,)

    cspec = join([join(i,":") for i in candleSpecs ],",")

    r = HTTP.get(string("https://",config.hostname,"/v3/accounts/",config.account,"/candles/latest"),
        ["Authorization" => string("Bearer ", config.token)];
        query = push!(Dict(),"candleSpecifications"=> cspec, kwargs...,),)

    temp = JSON3.read(r.body, latestCandles)

    #type coersions
    for latestc in temp.latestCandles
        for c in latestc.candles
        c = coerceCandleStick(config, c)
        end
    end

    return temp

end

# ------------------------------------------------------------------------------------
# /accounts/{accountID}/instrument/{instrument}/candle Endpoint
# ------------------------------------------------------------------------------------

"""
This endopoint is alsmost identical to /instruments/{instrument}/candles -> Not a priority

Only difference is the 'units' keyword argument

"""

end #module

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
@testset "Boykov Kolmogorov" begin
  @testset "Lattice graph" begin
    # image size
    sz = (9, 9)

    # number of pixels
    npix = prod(sz)

    # lattice graph
    G = DiGraph(npix+2)
    C = spzeros(npix+2, npix+2)

    # connect all pixels in the 9x9 image
    # with its immediate 4 neighbors
    for i in 1:sz[1]-1, j in 1:sz[2]
      u = LinearIndices(sz)[i,j]
      v = LinearIndices(sz)[i+1,j]
      add_edge!(G, u, v)
      add_edge!(G, v, u)
    end
    for i in 1:sz[1], j in 1:sz[2]-1
      u = LinearIndices(sz)[i,j]
      v = LinearIndices(sz)[i,j+1]
      add_edge!(G, u, v)
      add_edge!(G, v, u)
    end

    # create capacity for flow in the 4x4
    # subimage located at the top left
    for i in 1:3, j in 1:4
      u = LinearIndices(sz)[i,j]
      v = LinearIndices(sz)[i+1,j]
      C[u,v] = C[v,u] = 1
    end
    for i in 1:4, j in 1:3
      u = LinearIndices(sz)[i,j]
      v = LinearIndices(sz)[i,j+1]
      C[u,v] = C[v,u] = 1
    end

    # create capacity for flow in the 4x4
    # subimage located at the bottom right
    for i in 6:8, j in 6:9
      u = LinearIndices(sz)[i,j]
      v = LinearIndices(sz)[i+1,j]
      C[u,v] = C[v,u] = 1
    end
    for i in 6:9, j in 6:8
      u = LinearIndices(sz)[i,j]
      v = LinearIndices(sz)[i,j+1]
      C[u,v] = C[v,u] = 1
    end

    # create source node and connect it to the
    # leftmost column of the image
    # create target node and connect it to the
    # rightmost column of the image
    s = npix + 1
    t = npix + 2
    for i in 1:sz[1]
      u = LinearIndices(sz)[i,1]
      add_edge!(G, s, u)
      C[s,u] = C[u,s] = Inf
    end
    for i in 1:sz[1]
      u = LinearIndices(sz)[i,sz[2]]
      add_edge!(G, u, t)
      C[u,t] = C[t,u] = Inf
    end

    # now we are ready to start the flow
    flow, _, labels = maximum_flow(G, s, t, C, algorithm=BoykovKolmogorovAlgorithm())

    # because the two subimages are not connected
    # we must have zero flow from source to target
    @test flow == 0

    # the final cut represents the two disconnected
    # subimages filled with water of different color
    COLOR = reshape(labels[1:end-2], sz)
    @test COLOR == [
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  0  0  0  0  0  0  0  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
    ]

    # now let's create a bridge connecting the two
    # subimages to allow flow from source to target
    for (I1, I2) in [[(4,4), (5,4)], [(5,4), (5,5)],
                    [(5,5), (5,6)], [(5,6), (6,6)]]
      u = LinearIndices(sz)[I1...]
      v = LinearIndices(sz)[I2...]
      C[u,v] = C[v,u] = 1
    end

    flow, _, labels = maximum_flow(G, s, t, C, algorithm=BoykovKolmogorovAlgorithm())

    # because there is only one bridge,
    # the maximum flow allowed is one unit
    @test flow == 1

    # the final cut is unchanged compared to the previous one
    COLOR = reshape(labels[1:end-2], sz)
    @test COLOR == [
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  0  0  0  0  0  0  0  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
    ]

    # finally let's create a second bridge to increase
    # the maximum flow from one to two units
    for (I1, I2) in [[(4,4), (4,5)], [(4,5), (5,5)],
                    [(5,5), (6,5)], [(6,5), (6,6)]]
      u = LinearIndices(sz)[I1...]
      v = LinearIndices(sz)[I2...]
      C[u,v] = C[v,u] = 1
    end

    flow, _, labels = maximum_flow(G, s, t, C, algorithm=BoykovKolmogorovAlgorithm())

    # the maximum flow is now doubled
    @test flow == 2

    # the final cut is slightly different
    # near the corners of the two subimages
    COLOR = reshape(labels[1:end-2], sz)
    @test COLOR == [
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  1  1  1  0  0  0  0  2
      1  1  1  0  0  0  0  0  2
      1  0  0  0  0  0  0  0  2
      1  0  0  0  0  0  2  2  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
      1  0  0  0  0  2  2  2  2
    ]
  end

  @testset "Find path" begin
    # construct graph
    gg = lg.DiGraph(3)
    lg.add_edge!(gg, 1, 2)
    lg.add_edge!(gg, 2, 3)

    # source and sink terminals
    source, target = 1, 3

    for g in testdigraphs(gg)
      # default capacity
      capacity_matrix = LightGraphsFlows.DefaultCapacity(g)
      residual_graph = @inferred(LightGraphsFlows.residual(g))
      T = eltype(g)
      flow_matrix = zeros(T, 3, 3)
      TREE = zeros(T, 3)
      TREE[source] = T(1)
      TREE[target] = T(2)
      PARENT = zeros(T, 3)
      A = [T(source), T(target)]
      path = LightGraphsFlows.find_path!(
        residual_graph, source, target, flow_matrix,
        capacity_matrix, PARENT, TREE, A)

      @test path == [1, 2, 3]
    end
  end
end


text:
@testset "c`...`" begin
	@eval c``
	@test haskey(CBinding.CONTEXT_CACHE, @__MODULE__)
	ctx = CBinding.CONTEXT_CACHE[@__MODULE__]
	@test ctx isa CBinding.Context{:c}
	@test length(ctx.libs) == 1  # the julia process is the only lib

	@eval begin
		lib = startswith(Base.libm_name, "lib") ? Base.libm_name[4:end] : Base.libm_name
		libpath = normpath(joinpath(dirname(Base.julia_cmd().exec[1]), Base.LIBDIR, "julia"))
		c`-L$(libpath) -l$(lib)`
	end
	@test ctx !== CBinding.CONTEXT_CACHE[@__MODULE__]
	ctx = CBinding.CONTEXT_CACHE[@__MODULE__]
	@test ctx isa CBinding.Context{:c}
	@test length(ctx.libs) == 2
end


text:
abstract type RIPEMD_CTX end

mutable struct RIPEMD160_CTX <: RIPEMD_CTX
    state  :: Array{UInt32, 1} # length 5
    count  :: UInt64           # how many bytes we already ingested
    buffer :: Array{UInt8, 1}  # message is copied here, read as UInt32 in
                               # transform!
end

# length of the buffer
bytes_per_block(::Type{RIPEMD160_CTX}) = 64
words_per_block(::Type{RIPEMD160_CTX}) = 16

state_type(::Type{RIPEMD160_CTX}) = UInt32
digest_length(::Type{RIPEMD160_CTX}) = 20

function RIPEMD160_CTX()
    RIPEMD160_CTX(copy(INIT_STATE), 0 , zeros(UInt8, bytes_per_block(RIPEMD160_CTX)))
end