text:
export find_rotation_sets

struct rotation_element{F}
    axis::Vector{F}
    order::Int
end

function Base.:(==)(A::rotation_element, B::rotation_element)
    issame_axis(A.axis, B.axis) &&
    A.order == B.order
end

function issame_axis(a, b)
    """
    Checks if two axes are equivalent
    |a⋅b| = 1 assuming ||a|| = ||b|| = 1
    """
    A = normalize(a)
    B = normalize(b)
    d = abs(A ⋅ B)
    return isapprox(d, 1.0, atol=tol)
end

function rotation_set_intersection(rotation_set)
    """
    Finds the intersection of the elements in the sets
    within 'rotation_set'
    """
    out = deepcopy(rotation_set[1])
    len = size(rotation_set)[1]
    if len > 2
        for i = 1:size(rotation_set)[1]
            intersect!(out, rotation_set[i])
        end
    end
    return out
end

function find_rotation_sets(mol::Molecule, SEAs)
    """
    Loop through SEAs to label each set, and find principal axis and potential rotation orders
    as outlined in the flowchart in DOI: 10.1002/jcc.23493

    Returns a list of lists with each SEA's potential rotations
    """
    out_all_SEAs = []
    for sea in SEAs
        out_per_SEA = []
        len = size(sea.set)[1]
        if len < 2
            sea.label = "Single Atom"
        elseif len == 2
            sea.label = "Linear"
            sea.paxis = mol[sea.set[1]].xyz - mol[sea.set[2]].xyz
        else
            moit = eigenmoit(calcmoit([mol[i] for i in sea.set]))
            Ia, Ib, Ic = moit[1]
            Iav, Ibv, Icv = [moit[2][:,idx] for idx = 1:3]
            if isapprox(Ia, Ib, atol=tol) && isapprox(Ia, Ic, atol=tol)
                sea.label = "Spherical"
            elseif isapprox(Ia+Ib, Ic, atol=tol)
                paxis = Icv
                if isapprox(Ia, Ib, atol=tol)
                    sea.label = "Regular Polygon"
                    sea.paxis = paxis
                    for i = 2:len
                        if isfactor(len, i)
                            re = rotation_element(paxis, i)
                            push!(out_per_SEA, re)
                        end
                    end
                else
                    sea.label = "Irregular Polygon"
                    sea.paxis = paxis
                    for i = 2:len-1
                        if isfactor(len, i)
                            re = rotation_element(paxis, i)
                            push!(out_per_SEA, re)
                        end
                    end
                end
            else
                if !(isapprox(Ia, Ib, atol=tol) || isapprox(Ib, Ic, atol = tol))
                    sea.label = "Asymmetric Rotor"
                    for i in [Iav, Ibv, Icv]
                        re = rotation_element(i, 2)
                        push!(out_per_SEA, re)
                    end
                else
                    if Ia == Ib
                        paxis = Icv
                        sea.label = "Oblate Symmetric Top"
                        sea.paxis = paxis
                    else
                        paxis = Iav
                        sea.label = "Prolate Symmetric Top"
                        sea.paxis = paxis
                    end
                    k = div(len, 2)
                    for i = 2:k
                        if isfactor(k, i)
                            re = rotation_element(paxis, i)
                            push!(out_per_SEA, re)
                        end
                    end
                end
            end
            push!(out_all_SEAs, out_per_SEA)
        end
    end
    return out_all_SEAs
end

function isfactor(n, a)
    """
    Simple function, is a a factor of n? Return true/false
    """
    if n % a == 0
        return true
    else
        return false
    end
end

function find_rotations(mol::Molecule, rotation_set)
    """
    Find all actual rotation elements for the Molecule
    based off of potential rotation elements from 'find_rotation_sets'
    """
    molmoit = eigenmoit(calcmoit(mol))
    if molmoit[1][1] == 0.0 && molmoit[1][2] == molmoit[1][3]
        # This block executes if the Molecule is linear
        paxis = normalize(mol[1].xyz)
        re = rotation_element(paxis, 0)
        return [re]
    end
    rsi = rotation_set_intersection(rotation_set)
    out = []
    for i in rsi
        rmat = Molecules.Cn(i.axis, i.order)
        molB = Molecules.transform(mol, rmat)
        c = Molecules.isequivalent(molB, mol)
        if c
            push!(out, i)
        end
    end
    return out
end

function find_c2(mol::Molecule, SEAs)
    """
    Find C2 rotations as outlined in the paper DOI: 10.1002/jcc.23493
    Uses 'c2a', 'c2b', and 'c2c' to determine the presence of at least
    one C2 axis.
    """
    len = size(mol)[1]
    for sea in SEAs
        a = c2a(mol, sea)
        if a !== nothing
            return a
        else
            b = c2b(mol, sea)
            if b !== nothing
                return b
            else
                if sea.label == "Linear"
                    for sea2 in SEAs
                        if sea == sea2
                            continue
                        elseif sea2.label == "Linear"
                            c = c2c(mol, sea, sea2)
                            if c !== nothing
                                return c
                            end
                        end
                    end
                end
            end
        end
    end
    return nothing
end

function is_there_ortho_c2(mol::Molecule, cn_axis::Vector{T}, SEAs) where T
    """
    Finds C2 axes orthogonal to principal axis 'paxis' of Molecule
    When passing a principal axis to 'c2a', 'c2b', and 'c2c', returns
    true/false if there exists at least one C2 orthogonal to principal axis
    """
    len = size(mol)[1]
    for sea in SEAs
        if c2a(mol, cn_axis, sea)
            return true
        elseif c2b(mol, cn_axis, sea)
            return true
        elseif sea.label == "Linear"
            for sea2 in SEAs
                if sea == sea2
                    continue
                elseif sea2.label == "Linear"
                    if c2c(mol, cn_axis, sea, sea2)
                        return true
                    end
                end
            end
        end
    end
    return false
end

function num_C2(mol::Molecule, SEAs)
    """
    Counts number of unique C2 axes. Particularly useful for
    finding high symmetry groups. Uses C2 algorithms a and b
    but unlike previous functions, finds all unique C2 axes and
    returns the count.
    """
    axes = []
    for sea in SEAs
        a = all_c2a(mol, sea)
        if a !== nothing
            for i in a
                push!(axes, i)
            end
        end

        b = all_c2b(mol, sea)
        if b !== nothing
            for i in b
                push!(axes, i)
            end
        end
    end
    if size(axes)[1] < 1
        return nothing
    end
    unique_axes = [axes[1]]
    for i in axes
        check = true
        for j in unique_axes
            if issame_axis(i,j)
                check = false
                break
            end
        end
        if check
            push!(unique_axes, i)
        end
    end
    return size(unique_axes)[1]
end

function c2a(mol, sea)
    """
    Constructs a line going from the COM to the midpoint of every pair
    of atoms in a set of SEAs (assuming that midpoint is not the COM).
    That line is a potential C2 axis
    """
    len = size(sea.set)[1]
    for i = 1:len-1, j = i+1:len
        midpoint = mol[sea.set[i]].xyz + mol[sea.set[j]].xyz
        if midpoint == [0.0, 0.0, 0.0]
            continue
        else
            c2 = Molecules.Cn(midpoint, 2)
            molB = Molecules.transform(mol, c2)
            check = Molecules.isequivalent(mol, molB)
            if check
                return midpoint
            else
                continue
            end
        end
    end
    return nothing
end

function c2b(mol, sea)
    """
    Constructs a line from the COM to each atom in a set of SEAs.
    That line is a potential C2 axis.
    """
    len = size(sea.set)[1]
    for i = 1:len
        c2_axis = mol[sea.set[i]].xyz
        c2 = Molecules.Cn(c2_axis, 2)
        molB = Molecules.transform(mol, c2)
        check = Molecules.isequivalent(mol, molB)
        if check
            return c2_axis
        else
            continue
        end
    end
    return nothing
end

function c2c(mol, sea1, sea2)
    """
    Given two linear sets of SEAs (ij and kl), a potential C2 axis
    is the vector r_C2 = (i-j) x (k-l)
    """
    rij = mol[sea1.set[1]].xyz - mol[sea1.set[2]].xyz
    rkl = mol[sea2.set[1]].xyz - mol[sea2.set[2]].xyz
    c2_axis = cross(rij, rkl)
    c2 = Molecules.Cn(c2_axis, 2)
    molB = Molecules.transform(mol, c2)
    if Molecules.isequivalent(mol, molB)
        return c2_axis
    else
        return nothing
    end
end

function c2a(mol, cn_axis, sea)
    len = size(sea.set)[1]
    for i = 1:len-1, j = i+1:len
        midpoint = mol[sea.set[i]].xyz + mol[sea.set[j]].xyz
        if midpoint == [0.0, 0.0, 0.0]
            continue
        elseif issame_axis(midpoint, cn_axis)
            continue
        else
            c2 = Molecules.Cn(midpoint, 2)
            molB = Molecules.transform(mol, c2)
            check = Molecules.isequivalent(mol, molB)
            if check
                return true
            else
                continue
            end
        end
    end
    return false
end

function c2b(mol, cn_axis, sea)
    len = size(sea.set)[1]
    for i = 1:len
        c2_axis = mol[sea.set[i]].xyz
        if issame_axis(c2_axis, cn_axis)
            continue
        else
            c2 = Molecules.Cn(c2_axis, 2)
            molB = Molecules.transform(mol, c2)
            check = Molecules.isequivalent(mol, molB)
            if check
                return true
            else
                continue
            end
        end
    end
    return false
end

function c2c(mol, cn_axis, sea1, sea2)
    rij = mol[sea1.set[1]].xyz - mol[sea1.set[2]].xyz
    rkl = mol[sea2.set[1]].xyz - mol[sea2.set[2]].xyz
    c2_axis = cross(rij, rkl)
    if issame_axis(c2_axis, cn_axis)
        return false
    else
        c2 = Molecules.Cn(c2_axis, 2)
        molB = Molecules.transform(mol, c2)
        return Molecules.isequivalent(mol, molB)
    end
end

function all_c2a(mol, sea)
    out = []
    len = size(sea.set)[1]
    for i = 1:len-1, j = i+1:len
        midpoint = mol[sea.set[i]].xyz + mol[sea.set[j]].xyz
        if midpoint == [0.0, 0.0, 0.0]
            continue
        else
            c2 = Molecules.Cn(midpoint, 2)
            molB = Molecules.transform(mol, c2)
            check = Molecules.isequivalent(mol, molB)
            if check
                push!(out, midpoint)
            else
                continue
            end
        end
    end
    if size(out)[1] < 1
        return nothing
    end
    return out
end

function all_c2b(mol, sea)
    out = []
    len = size(sea.set)[1]
    for i = 1:len
        c2_axis = mol[sea.set[i]].xyz
        c2 = Molecules.Cn(c2_axis, 2)
        molB = Molecules.transform(mol, c2)
        check = Molecules.isequivalent(mol, molB)
        if check
            push!(out, c2_axis)
        else
            continue
        end
    end
    if size(out)[1] < 1
        return nothing
    end
    return out
end

function highest_ordered_axis(rotations::Vector{Any})
    # Returns the largest element of the set of rotation orders
    len = size(rotations)[1]
    ns = []
    for i = 1:len
        push!(ns, rotations[i].order)
    end
    return sort(ns)[len]
end

function is_there_sigmah(mol::Molecule, paxis)
    """
    Checks for reflection across plane with normal vector
    equal to principal rotation axis
    """
    σh = Molecules.reflection_matrix(paxis)
    molB = Molecules.transform(mol, σh)
    return Molecules.isequivalent(mol, molB)
end

function is_there_sigmav(mol, SEAs, paxis)
    """
    Checks for the presence of reflection planes that
    contain the principal c2_axis
    """
    axes = []
    for sea in SEAs
        len = size(sea.set,1)
        if len < 2
            continue
        end
        A = sea.set[1]
        for i = 2:len
            B = sea.set[i]
            #n = normalize(mol[A].xyz - mol[B].xyz)
            n = mol[A].xyz - mol[B].xyz
            σ = Molecules.reflection_matrix(n)
            molB = Molecules.transform(mol, σ)
            check = Molecules.isequivalent(mol, molB)
            if check
                push!(axes, n)
            else
                continue
            end
        end
    end
    if size(axes,1) < 1
        if mol_is_planar(mol)
            return true
        else
            return false
        end
    end
    unique_axes = [axes[1]]
    # This part is VERY inefficient!!!
    for i in axes
        check = true
        for j in unique_axes
            if issame_axis(i,j)
                check = false
                break
            end
            if check
                push!(unique_axes, i)
            end
        end
    end
    for i in unique_axes
        if issame_axis(i, paxis)
            continue
        else
            return true
        end
    end
    return false
end

function mol_is_planar(mol)
    """
    Checks if Molecule is planar
    """
    t = eltype(mol[1].xyz)
    len = size(mol,1)
    mat = zeros(Float64, (3,len))
    for i = 1:len
        mat[:,i] = mol[i].xyz
    end
    rank = LinearAlgebra.rank(mat)
    if rank < 3
        return true
    end
    return false
end


text:
import Base.convert
using MAT
using ForwardDiff
using Optim
using Polynomials
using JLD
using ROCAnalysis
using PyPlot
import Base.convert
using MAT
using ForwardDiff
using Optim
using Polynomials
using JLD

@everywhere path="/home/emma/github/fit_behaviour_julia/"

path_functions=path*"functions/"
path_save=path*"results/"

include(path_functions*"parameter_functions.jl")
include(path_functions*"functions_julia.jl")
include(path_functions*"simulations.jl")

Nstim=7000
Nframes=10
Tframe=.2
T=Tframe*Nframes
tau=1   #### this is always 1 in the theory.
dt_sim=tau/100.
NT=T/dt_sim
NTframes=NT/Nframes
sigma_a=.3

# model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:st_f=>st_logistic_time)
# param=Dict("c2"=>[0.7,0.05], "st"=>[2.0,0.5,0.5],"c4"=>[1.0],"sigma_a"=>sigma_a)

# model_num="2"
# model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:st_f=>st_logistic_time,:ro_f=>read_out_soft_sim)
# param=Dict("c2"=>[0.7,0.05], "st"=>[2.0,0.5,0.05],"c4"=>[1.0],"sigma_a"=>sigma_a,"ro"=>[3,0.3])

# model_num="3"
# model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_const,:st_f=>st_const,:ro_f=>read_out_TW_sim)
# param=Dict("c2"=>[-0.025,1],"c4"=>[-0.5],"c6"=>[0.4], "st"=>[1.],"sigma_a"=>sigma_a)

model_num="6"
model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect_sim)

param=Dict("c2"=>[1,0],"c4"=>[1],"x0"=>[0.0,0.0], "st"=>[2.,2,0.,0.],"bias_d"=>[0.],
"hbias_d"=>[0.0],"sigma_a"=>sigma_a)

Nsigmas=5
sigma_max=1

sigmas=exp10.(linspace(-1.17,0.,5))
stim=randn(Nstim,Nframes)
stim_test=randn(Nstim,Nframes)

internal_noise=sigma_a*randn(Nstim,Int(NT))
internal_noise_test=sigma_a*randn(Nstim,Int(NT))
stim_trans=zeros(Nstim,Nframes)

s=zeros(Nframes)

for istim in 1:Nstim
    stim[istim,:]=sigmas[(istim-1)%Nsigmas+1]*stim[istim,:]
    stim_test[istim,:]=sigmas[(istim-1)%Nsigmas+1]*stim_test[istim,:]
end

for istim in 1:Nstim
    # stim_trans[istim,:]=st_logistic_time(stim[istim,:],stim_trans[istim,:],param["st"])
    #st_logistic_time(stim[istim,:],s,param["st"])
    model_f[:st_f](stim[istim,:],s,param["st"])
    stim_trans[istim,:]=s
end

figure()
plot(stim[:],stim_trans[:],".")
#
#
d,xfinal=simulation_Nframe_general(param,stim,internal_noise,T,tau,model_f)
d_test,_=simulation_Nframe_general(param,stim_test,internal_noise_test,T,tau,model_f)
println("hola ",mean(0.5*(d+1)),mean(.5*(1+sign(xfinal))))
for isigma in 1:Nsigmas

    indices=isigma:Nsigmas:Nstim
end

#######Computing PR############
Ntrials=1000
Nstim_PR=100
d_pr=zeros(Ntrials,Nstim_PR)
d_pr_test=zeros(Ntrials,Nstim_PR)
x_pr=zeros(Ntrials)
for itrial in 1:Ntrials
    internal_noise=sigma_a*randn(Nstim_PR,Int(Nframes*NTframes))
    d_pr[itrial,:],x_pr=simulation_Nframe_general(param,stim[1:Nstim_PR,:],internal_noise,T,tau, model_f)
    internal_noise=sigma_a*randn(Nstim_PR,Int(Nframes*NTframes))
    d_pr_test[itrial,:],_=simulation_Nframe_general(param,stim_test[1:Nstim_PR,:],internal_noise,T,tau,model_f)
end

####compute trace ####
N_traces=20
istim_pdf=2
internal_noise=sigma_a*randn(N_traces,Int(Nframes*NTframes))
stim_pdf=transpose(repeat(stim[istim_pdf,:],1,N_traces))
d_traces,x=simulation_Nframe_general_trace(param,stim_pdf,internal_noise,T,tau, model_f)

PR_sim=zeros(Float64,Nstim_PR)
PR_sim_test=zeros(Float64,Nstim_PR)

for istim in 1:Nstim_PR
    PR_sim[istim]= mean( 0.5.*(d_pr[:,istim]+1))
    PR_sim_test[istim]= mean( 0.5.*(d_pr_test[:,istim]+1))
end

# path="/home/genis/model_fitting_julia/"
# filename="synthetic_data_DW_Nstim"*string(Nstim)*"_5sigmas_test_largeParam.jld"
filename="synthetic_data_DW_Nstim"*string(Nstim)*"_model"*model_num*".jld"

save_dict=Dict("d"=>d,"stim"=>stim,"param"=>param,"sigma_a"=>sigma_a,"Tframe"=>T,
"PR"=>PR_sim,"kernel"=>kernel_sigma, "PR_sim_test"=>PR_sim_test, "stim_test"=>stim_test,
"d_test"=>d_test,"Tframe"=>Tframe,"x_trace"=>x,"istim_pdf"=>istim_pdf)
JLD.save(path_save*filename,save_dict)

ext: jl

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
export hamming, profile, consensus, numseq

"""
	countgaps(s::BioSequence)
"""
function countgaps(s::BioSequence)
	return BioSequences.count(x->isgap(x)||isambiguous(x), s)
end
countgaps(s::AbstractStrain) = countgaps(s.seq)

function hamming_ignore_ambiguous(a::BioSequence, b::BioSequence)
	h = 0
	for (x,y) in zip(a,b)
		if x != y && !isambiguous(x) && !isambiguous(y)
			h += 1
		end
	end
	return h
end
function hamming_ambiguous_mismatch(a::BioSequence, b::BioSequence)
	h = 0
	for (x,y) in zip(a,b)
		if x != y
			h += 1
		end
	end
	return h
end
"""
	hamming(a::BioSequence, b::BioSequence, [normalization=false]; ambiguous=:ignore)
	hamming(a::AbstractStrain, b::AbstractStrain [, normalization=false]; ambiguous=:ignore)

`ambiguous` keyword:
- `:ignore` : ambiguous symbols are considered equal to anything .
- `:mismatch` : ambiguous symbols behave like regular symbols.
"""
function hamming(a::BioSequence, b::BioSequence,
	normalization=false;
	ambiguous=:ignore)
	if ambiguous == :ignore && normalization
		return hamming_ignore_ambiguous(a,b)/length(a)
	elseif ambiguous == :ignore && !normalization
		return Float64(hamming_ignore_ambiguous(a,b))
	elseif ambiguous == :mismatch && normalization
		return hamming_ambiguous_mismatch(a,b)/length(a)
	elseif ambiguous == :mismatch && !normalization
		return Float64(hamming_ambiguous_mismatch(a,b))
	else
		@error "Incorrect arguments."
	end
end
hamming(a::AbstractStrain, b::AbstractStrain, normalization=false; ambiguous=:ignore) = hamming(a.seq, b.seq, normalization, ambiguous=ambiguous)
"""
	hamming(P::Profile{A}, Q::Profile{A}, normalization=false;
				ambiguous=:ignore) where A <: BioSymbol
"""
function hamming(P::Profile{A}, Q::Profile{A}, normalization=false;
				ambiguous=:ignore) where A <: BioSymbol
	length(P) != length(Q) && @error "Profiles of different lengths"
	out = 0.
	for (i,(p,q)) in enumerate(zip(P,Q))
		val = 1.
		for a in BioSequences.alphabet(A)
			if ambiguous == :ignore && isambiguous(a)
				val -= get(p,a,0.) + get(q,a,0.) - get(p,a,0.)*get(q,a,0.)
			else
				val -= get(p,a,0.)*get(q,a,0.)
			end
		end
		out += val
	end
	if normalization
		out /= length(P)
	end
	return out
end
function hamming(P::Profile, s::BioSequence, normalization = false; ambiguous = :ignore)
	length(P) != length(s) && @error "Sequence and profile of different lengths"
	out = 0.
	for (i,(a,p)) in enumerate(zip(s,P))
		out += 1.
		for b in keys(p)
			if b == a || (ambiguous==:ignore && (isambiguous(a)||isambiguous(b)))
				out -= p[b]
			end
		end
	end
	if normalization
		out /= length(P)
	end
	return out
end
hamming(s::BioSequence, P::Profile) = hamming(P,s)
hamming(s::Strain, P::Profile) = hamming(s.seq, P)
hamming(P::Profile, s::Strain) = hamming(s.seq, P)

"""
	consensus(P::Profile{A}) where A<:BioSymbol
	consensus(X::Array{<:BioSequence,1})

Return the consensus *sequence*, of type `LongSequence{alphabet}` where `eltype(alphabet)==A`. Does **NOT** handle ties. In case of a tie, the chosen symbol is undetermined.

	consensus(X::Array{Strain} [, data = Dict(:strain=>"consensus")])

Return consensus of `X`, with field `data`. Does **NOT** handle ties. In case of a tie, the chosen symbol is undetermined.
"""
function consensus(P::Profile{A}) where A<:BioSymbol
	seq = LongSequence{symbol_to_alphabet(A)}(length(P))
	for (i,p) in enumerate(P)
		if !isempty(p)
			seq[i] = findmax(p)[2]
		else
			seq[i] = ambiguous(A)
		end
	end
	return seq
end
consensus(X::Array{<:BioSequence,1}) = consensus(Profile(X))
function consensus(X::Array{<:Strain,1}, data = Dict(:strain=>"consensus"))
	p = Profile(X)
	seq = consensus(p)
	return Strain(seq, data)
end

"""
	numseq(s::LongSequence{AminoAcidAlphabet})
	numseq(S::Array{LongSequence{AminoAcidAlphabet}})
"""
function numseq(s::LongSequence{AminoAcidAlphabet})
	σ = Array{Int64,1}(undef, length(s))
	for (i,a) in enumerate(s)
		σ[i] = haskey(aa2int, a) ? aa2int[a] : 1
	end
	return σ
end
function numseq(S::Array{LongSequence{AminoAcidAlphabet}})
	if isempty(S)
		out = Array{Int64,2}(undef, 0, 0)
	else
		out = Array{Int64,2}(undef, length(S), length(S[1]))
		for (i,s) in enumerate(S)
			out[i,:] .= numseq(s)
		end
	end
	return out
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
module Starlight

using Reexport
@reexport using Base: ReentrantLock, lock, unlock
@reexport using FileIO: load
@reexport using Pkg.Artifacts
@reexport using LazyArtifacts
@reexport using DataStructures: Queue, enqueue!, dequeue!
@reexport using DataFrames
@reexport using Colors
@reexport using SimpleDirectMediaLayer
@reexport using SimpleDirectMediaLayer.LibSDL2

export handleMessage, sendMessage, sendMessageTo, listenFor, listenForFrom, handleException, dispatchMessage
export System, App, awake!, shutdown!, system!, on, off, cat
export app, cfg, str_to_clrnt, get_env_int, get_env_str, get_env_clr, get_env_flt, get_env_bl

import DotEnv
cfg = DotEnv.config()

str_to_clrnt(s) = parse(Colorant, s)
get_env_int(n, d) = parse(Int, get(ENV, n, string(d)))
get_env_str(n, d) = get(ENV, n, string(d)) # no parse necessary since all env variables are strings by default
get_env_clr(n, d) = str_to_clrnt(get(ENV, n, string(d)))
get_env_flt(n, d) = parse(Float64, get(ENV, n, string(d)))
get_env_bl(n, d) = parse(Bool, get(ENV, n, string(d)))

const listeners = Dict{DataType, Set{Any}}()
const messages = Channel(Inf)
const from = Dict{Any, Dict{DataType, Set{Any}}}()

const listener_lock = ReentrantLock()

handleMessage(l, m) = nothing

function sendMessage(m)
  # drop if no one's listening
  if haskey(listeners, typeof(m))
    put!(messages, m)
  end
end

# allow direct invocation of message handlers
# for cases when we know all the recipients in
# advance and there would otherwise be too many
# listeners (i.e. collision notifications)
function sendMessageTo(m, targets...)
  for to in targets handleMessage(to, m) end
end

function add_listener!(e::Any, d::DataType)
  if !haskey(listeners, d) listeners[d] = Set{Any}() end
  push!(listeners[d], e)
end

function listenFor(e::Any, d::DataType)
  lock(listener_lock)
  add_listener!(e, d)
  unlock(listener_lock)
end

# used in cases where we care about the sender
# as well as the message
function listenForFrom(e::Any, d::DataType, e2::Any)
  lock(listener_lock)
  add_listener!(e, d)
  if !haskey(from, e) from[e] = Dict{DataType, Set{Any}}() end
  if !haskey(from[e], d) from[e][d] = Set{Any}() end
  push!(from[e][d], e2)
  unlock(listener_lock)
end

# TODO implement "unlisten" function and test
# NOTE this would only ever be used by the ECS
# when removing an entity, it's too much to ask
# users to write finalizers that do this, i.e.
# this is not terribly important right now

function handleException()
  for s in stacktrace(catch_backtrace())
    println(s)
  end
  rethrow()
end

# uses single argument to support
# being called as a job by a Clock,
# see Clock.jl for that interface
function dispatchMessage(arg)
  try
    m = take!(messages) # NOTE messages are fully processed in the order they are received
    d = typeof(m)
    @debug "dequeued message $(m) of type $(d)"
    if haskey(listeners, d)
      # Threads.@threads doesn't work on raw sets
      # because it needs to use indexing to split
      # up memory, i work around it this way
      Threads.@threads for l in Vector([listeners[d]...])
        # only call handler if the sender is
        # being listened for, i.e. matches one
        # the recipient cares about, or the
        # recipient doesn't care
        if haskey(from, l) && haskey(from[l], d)
          if !hasproperty(m, :from) || m.from ∉ from[l][d]
            continue
          end
        end
        handleMessage(l, m)
      end
    end
  catch
    handleException()
  end
end

abstract type System end

# these functions are supposed to return
# whether a value indicating whether the
# system is still running or not, except
# for App where it returns a vector of
# booleans indicating whether each system
# is still running
awake!(s::System) = true
shutdown!(s::System) = false

include("Clock.jl")
include("ECS.jl")
include("Entities.jl")
include("Scene.jl")
include("SDL.jl")
include("Audio.jl")

mutable struct App <: System
  systems::Dict{DataType, System}
  running::Vector{Bool}
  function App()
    # singleton pattern from Tom Kwong's book
    global app
    global app_lock
    lock(app_lock)
    try
      if !isassigned(app)

        a = new(Dict(), Vector{Bool}())

        system!(a, clk)
        system!(a, ecs)
        system!(a, sdl)
        system!(a, scn)

        a.running = [false for s in keys(a.systems)]

        app[] = finalizer(shutdown!, a)

      end

      # clock
      clk.freq = get_env_flt("CLOCK_FREQ", 0.01667) # 60hz

      fire_sec = get_env_bl("CLOCK_FIRE_SEC", false)
      fire_msec = get_env_bl("CLOCK_FIRE_MSEC", false)
      fire_usec = get_env_bl("CLOCK_FIRE_USEC", false)
      fire_nsec = get_env_bl("CLOCK_FIRE_NSEC", false)
      if fire_sec
        push!(Clock.periodic_jobs, (1e9, RT_SEC, "second"))
      end
      if fire_msec
        push!(Clock.periodic_jobs, (1e6, RT_MSEC, "millisecond"))
      end
      if fire_usec
        push!(Clock.periodic_jobs, (1e3, RT_USEC, "microsecond"))
      end
      if fire_nsec
        push!(Clock.periodic_jobs, (1.0, RT_NSEC, "nanosecond"))
      end

      # sdl
      sdl.bgrd = to_ARGB(get_env_clr("BACKGROUND_COLOR", "gray"))
      sdl.wdth = get_env_int("WINDOW_WIDTH", 800)
      sdl.hght = get_env_int("WINDOW_HEIGHT", 400)
      sdl.ttl = get_env_str("TITLE", "Starlight.jl")

    catch e
      rethrow()
    finally
      unlock(app_lock)
      return app[]
    end
  end
end

const app = Ref{App}()
const app_lock = ReentrantLock()

on(a::App) = all(a.running)
off(a::App) = all(!r for r in a.running)
# shrodinger's app is neither on nor off,
# i.e. some system is not synchronized with the others
cat(a::App) = !on(a) && !off(a)

system!(a::App, s::System) = a.systems[typeof(s)] = s
# note that if running from a script the app will
# still exit when julia exits, it will never block.
# figuring out whether/how to keep it alive is
# on the user.
function awake!(a::App)
  if !on(a)
    job!(a.systems[Clock], dispatchMessage) # this could be parallelized if not for mqueue_lock
    a.running = map(awake!, values(a.systems))
  end
  return a.running
end
function shutdown!(a::App)
  if !off(a)
    a.running = map(shutdown!, values(a.systems))
  end
  return a.running
end

awake!() = awake!(app[])
shutdown!() = shutdown!(app[])

end


text:
#!/usr/bin/env julia
using Gtk, Gtk.ShortNames

win = Window("Count Click")
v = Box(:v)
l = Label("You clicked 0 times.")
b = Button("Click!")
push!(win, v)
push!(v, l)
push!(v, b)
set_gtk_property!(v, :expand, l, true)

showall(win)

niter = 0
function click()
    global niter += 1
    set_gtk_property!(l, :label, "You clicked $niter times.")
    return nothing
end

signal_connect(x -> click(), b, "clicked")

if !isinteractive()
    c = Condition()
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
# ESCB-2, "efficient sampling for combinatorial bandits".
# Based on https://arxiv.org/abs/1502.03475.

abstract type ESCB2OptimisationAlgorithm end
function optimise_linear_sqrtlinear(instance::CombinatorialInstance{T}, ::ESCB2OptimisationAlgorithm,
                                    linear::Dict{T, Float64}, sqrtlinear::Dict{T, Float64}) where T end

mutable struct ESCB2 <: Policy
  algo::ESCB2OptimisationAlgorithm
end

mutable struct ESCB2Details <: PolicyDetails
  n_iterations::Int
  best_objective::Float64
  solver_time::Float64

  function ESCB2Details()
    return new(0, 0.0, 0.0)
  end
end

function _escb2_confidence_values(state::State{T}) where T
  t = state.arm_counts
  w = state.arm_average_reward
  n = state.round

  d = length(t) # n_arms
  fn = log(n) + 4 * d * log(log(n))
  fn = max(0, fn) # f is negative for some values of d and small values of n (sometimes as small as 2 or 3).

  return Dict(k => (fn / 2) / v for (k, v) in t)
end

# TODO: shouldn't it go to a more generic function? Many algorithm work with indices, so maybe there is something to do... but only for black-box algorithms?
function escb2_index(state::State{T}, solution::Vector{T}) where T
  w = state.arm_average_reward
  s2 = _escb2_confidence_values(state)
  return escb2_index(w, s2, solution)
end
function escb2_index(linear::Dict{T, Float64}, sqrtlinear::Dict{T, Float64}, solution::Vector{T}) where T
  if length(solution) == 0
    return 0.0
  end

  avg_reward = sum(linear[arm] for arm in keys(linear) if arm in solution)
  confidence = sum(sqrtlinear[arm] for arm in keys(sqrtlinear) if arm in solution)
  return avg_reward + sqrt(confidence)
end

function choose_action(instance::CombinatorialInstance{T}, policy::ESCB2, state::State{T}; with_trace::Bool=false) where T
  # For initialisation, if some arms have never been tried, force them to be tried.
  # Without this, the standard deviation term cannot exist.
  if any(v == 0 for v in values(state.arm_counts))
    weights = Dict(arm => iszero(state.arm_counts[arm]) ? 1.0 : 0.0 for arm in keys(state.arm_counts))

    t0 = time_ns()
    sol = solve_linear(instance, weights)
    t1 = time_ns()

    if with_trace
      run_details = ESCB2Details()
      run_details.n_iterations = 1
      run_details.solver_time = (t1 - t0) / 1_000_000_000
      return sol, run_details
    else
      return sol
    end
  end

  # If all arms have been tried, can compute the two parts of the objective function: the average reward and a standard deviation.
  w = state.arm_average_reward
  s2 = _escb2_confidence_values(state)

  # Solve the maximisation w^T x + \sqrt{s2^T x} with the required algorithm.
  return optimise_linear_sqrtlinear(instance, policy.algo, w, s2, with_trace=with_trace)
end


text:
using CSV
using DataFrames

function prepare_empty_col(sz, tp=Any)
    arr = Array{tp, 1}(undef, sz)
    for i in 1:sz
        arr[i] = ""
#         arr[i] = nothing
    end
    return arr
end

function add_empty_col!(dataframe::DataFrame, col_name)
#     col = prepare_empty_col(size(tbl, 1), String)
    col = prepare_empty_col(size(dataframe, 1))
    if typeof(col_name) == Symbol
        column_name = col_name
    else
        column_name = Symbol(col_name)
    end
#     tbl[column_name] = 1:size(tbl,1)
    dataframe[column_name] = col
#     display(tbl)
    return dataframe

end

function add_missing_cols!(dataframe::DataFrame, to_add::Dict)

    nms = [String(name) for name in names(dataframe)]
    # @info "names as string" nms #, nms
    ks = [String(name) for name in keys(to_add)]

    allkeys = collect(Set(vcat(nms, ks)))
    # @info "names as string" allkeys
    for key in allkeys
        if key in nms
            # println("key '$key' exists")
        else
            add_empty_col!(dataframe, key)
            # println("adding empty col for key '$key'")
        end

    end
    return dataframe
end

function prepare_empty_row(tbl)
    new_row = Dict()
    for nm in names(tbl)
        if eltype(tbl[nm]) <: Number
#             println("$nm is subtype Number")
            # new_row[nm] = 0
            new_row[nm] = missing
        elseif eltype(tbl[nm]) == String
#             println("$nm is String")
            new_row[nm] = ""
        else
            new_row[nm] = missing

        end
    end

    return new_row
end

function dict_string_keys_to_symbols(to_add)
    new_to_add = Dict()
    for (key, val) in to_add
    #     println(key)
        new_to_add[Symbol(key)] = val
    end
    return new_to_add
end

function extend_row_type!(tbl1::DataFrame, tbl2::Union{Dict, DataFrame})
    # Extend data type
    for name in names(tbl1)
    # name = :one
        tbtype = Array{Union{
            eltype(tbl1[name]),
            eltype(tbl2[name])
        },1}
        tbtype
        tbl1[name] = convert(tbtype, tbl1[name])

    end
    return tbl1
end

function add_row_to_csv(row::Tuple, filename)
    return add_row_to_csv(Dict(row), filename)
end

function add_row_to_csv(row::Dict, filename)
    if isfile(filename)
        tbl = CSV.read(filename)
        tbl = DataFrame(tbl)
        @debug "display tbl" tbl
        to_add = row
        add_missing_cols!(tbl, to_add)
        @debug "missing cols" tbl
        new_row = prepare_empty_row(tbl)
        new_to_add = dict_string_keys_to_symbols(to_add)
        # print("new row: $new_row\n")
        # print("new to add: $new_to_add\n")
         merge!(new_row, new_to_add)
        # print("tbl $tbl")
        # print("new_row: $new_row")
        extend_row_type!(tbl, new_row)
        push!(tbl, new_row)
    else
        tbl = DataFrame(row)
    end
    CSV.write(filename, tbl)
    return tbl
end


text:
const CompoundPeriodZero = Nanosecond(0)
const CompoundPeriodEmpty = CompoundPeriod(CompoundPeriodZero)

canonical(x::Nanosecond) = Dates.canonicalize(Microsecond(0)+x)
canonical(x::Microsecond) = Dates.canonicalize(x+Nanosecond(0))
canonical(x::Millisecond) = Dates.canonicalize(x+Microsecond(0))
canonical(x::Second) = Dates.canonicalize(x+Millisecond(0))
canonical(x::Minute) = Dates.canonicalize(x+Second(0))
canonical(x::Hour) = Dates.canonicalize(x+Minute(0))
canonical(x::Day) = Dates.canonicalize(x+Hour(0))
canonical(x::Week) = Dates.canonicalize(x+Day(0))
canonical(x::Month) = Dates.canonicalize(Year(0)+x)
canonical(x::Year) = Dates.canonicalize(x+Month(0))

function canonical(x::CompoundPeriod)
    c = Dates.canonicalize(x)
    length(c) == 0 && return CompoundPeriodZero
    res = c.periods[1]
    !signbit(res) && return c
    c = c - res
    for p in c.periods
        x, y = fldmod(p)
        res = res + x + y
    end
    return res
end

fldmod2(x::Int64, y::Month) = (Year(x), y)
fldmod2(x::Int64, y::Day) = (Week(x), y)
fldmod2(x::Int64, y::Hour) = (Day(x), y)
fldmod2(x::Int64, y::Minute) = (Hour(x), y)
fldmod2(x::Int64, y::Second) = (Minute(x), y)
fldmod2(x::Int64, y::Millisecond) = (Second(x), y)
fldmod2(x::Int64, y::Microsecond) = (Millisecond(x), y)
fldmod2(x::Int64, y::Nanosecond) = (Microsecond(x), y)

fldmod(x::Year)::Year = x
fldmod(x::Month)::Tuple{Year,Month} = fldmod2(fldmod(x, Month(MONTHS_PER_YEAR))...,)
fldmod(x::Week)::Week = x
fldmod(x::Day)::Tuple{Week,Day} = fldmod2(fldmod(x, Day(DAYS_PER_WEEK))...,)
fldmod(x::Hour)::Tuple{Day,Hour} = fldmod2(fldmod(x, Hour(HOURS_PER_DAY))...,)
fldmod(x::Minute)::Tuple{Hour,Minute} = fldmod2(fldmod(x, Minute(MINUTES_PER_HOUR))...,)
fldmod(x::Second)::Tuple{Minute,Second} = fldmod2(fldmod(x, Second(SECONDS_PER_MINUTE))...,)
fldmod(x::Millisecond)::Tuple{Second,Millisecond} = fldmod2(fldmod(x, Millisecond(MILLISECONDS_PER_SECOND))...,)
fldmod(x::Microsecond)::Tuple{Millisecond,Microsecond} = fldmod2(fldmod(x, Microsecond(MICROSECONDS_PER_MILLISECOND))...,)
fldmod(x::Nanosecond)::Tuple{Microsecond,Nanosecond} = fldmod2(fldmod(x, Nanosecond(NANOSECONDS_PER_MICROSECOND))...,)


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

 ## Basic functions ##

isreal(x::AbstractArray) = all(isreal,x)
iszero(x::AbstractArray) = all(iszero,x)
isreal(x::AbstractArray{<:Real}) = true
all(::typeof(isinteger), ::AbstractArray{<:Integer}) = true

## Constructors ##

"""
    vec(a::AbstractArray) -> Vector

Reshape array `a` as a one-dimensional column vector.

```jldoctest
julia> a = [1 2 3; 4 5 6]
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

julia> vec(a)
6-element Array{Int64,1}:
 1
 4
 2
 5
 3
 6
```
"""
vec(a::AbstractArray) = reshape(a,_length(a))
vec(a::AbstractVector) = a

_sub(::Tuple{}, ::Tuple{}) = ()
_sub(t::Tuple, ::Tuple{}) = t
_sub(t::Tuple, s::Tuple) = _sub(tail(t), tail(s))

"""
    squeeze(A, dims)

Remove the dimensions specified by `dims` from array `A`.
Elements of `dims` must be unique and within the range `1:ndims(A)`.
`size(A,i)` must equal 1 for all `i` in `dims`.

```jldoctest
julia> a = reshape(collect(1:4),(2,2,1,1))
2×2×1×1 Array{Int64,4}:
[:, :, 1, 1] =
 1  3
 2  4

julia> squeeze(a,3)
2×2×1 Array{Int64,3}:
[:, :, 1] =
 1  3
 2  4
```
"""
function squeeze(A::AbstractArray, dims::Dims)
    for i in 1:length(dims)
        1 <= dims[i] <= ndims(A) || throw(ArgumentError("squeezed dims must be in range 1:ndims(A)"))
        size(A, dims[i]) == 1 || throw(ArgumentError("squeezed dims must all be size 1"))
        for j = 1:i-1
            dims[j] == dims[i] && throw(ArgumentError("squeezed dims must be unique"))
        end
    end
    d = ()
    for i = 1:ndims(A)
        if !in(i, dims)
            d = tuple(d..., size(A, i))
        end
    end
    reshape(A, d::typeof(_sub(size(A), dims)))
end

squeeze(A::AbstractArray, dim::Integer) = squeeze(A, (Int(dim),))

## Unary operators ##

conj(x::AbstractArray{<:Real}) = x
conj!(x::AbstractArray{<:Real}) = x

real(x::AbstractArray{<:Real}) = x
imag(x::AbstractArray{<:Real}) = zero(x)

+(x::AbstractArray{<:Number}) = x
*(x::AbstractArray{<:Number,2}) = x

# index A[:,:,...,i,:,:,...] where "i" is in dimension "d"

"""
    slicedim(A, d::Integer, i)

Return all the data of `A` where the index for dimension `d` equals `i`. Equivalent to
`A[:,:,...,i,:,:,...]` where `i` is in position `d`.

```jldoctest
julia> A = [1 2 3 4; 5 6 7 8]
2×4 Array{Int64,2}:
 1  2  3  4
 5  6  7  8

julia> slicedim(A,2,3)
2-element Array{Int64,1}:
 3
 7
```
"""
function slicedim(A::AbstractArray, d::Integer, i)
    d >= 1 || throw(ArgumentError("dimension must be ≥ 1"))
    nd = ndims(A)
    d > nd && (i == 1 || throw_boundserror(A, (ntuple(k->Colon(),nd)..., ntuple(k->1,d-1-nd)..., i)))
    A[setindex(indices(A), i, d)...]
end

function flipdim(A::AbstractVector, d::Integer)
    d == 1 || throw(ArgumentError("dimension to flip must be 1"))
    reverse(A)
end

"""
    flipdim(A, d::Integer)

Reverse `A` in dimension `d`.

```jldoctest
julia> b = [1 2; 3 4]
2×2 Array{Int64,2}:
 1  2
 3  4

julia> flipdim(b,2)
2×2 Array{Int64,2}:
 2  1
 4  3
```
"""
function flipdim(A::AbstractArray, d::Integer)
    nd = ndims(A)
    1 ≤ d ≤ nd || throw(ArgumentError("dimension $d is not 1 ≤ $d ≤ $nd"))
    if isempty(A)
        return copy(A)
    end
    inds = indices(A)
    B = similar(A)
    nnd = 0
    for i = 1:nd
        nnd += Int(length(inds[i])==1 || i==d)
    end
    indsd = inds[d]
    sd = first(indsd)+last(indsd)
    if nnd==nd
        # flip along the only non-singleton dimension
        for i in indsd
            B[i] = A[sd-i]
        end
        return B
    end
    alli = [ indices(B,n) for n in 1:nd ]
    for i in indsd
        B[[ n==d ? sd-i : alli[n] for n in 1:nd ]...] = slicedim(A, d, i)
    end
    return B
end

function circshift(a::AbstractArray, shiftamt::Real)
    circshift!(similar(a), a, (Integer(shiftamt),))
end
circshift(a::AbstractArray, shiftamt::DimsInteger) = circshift!(similar(a), a, shiftamt)
"""
    circshift(A, shifts)

Circularly shift the data in an array. The second argument is a vector giving the amount to
shift in each dimension.

```jldoctest
julia> b = reshape(collect(1:16), (4,4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> circshift(b, (0,2))
4×4 Array{Int64,2}:
  9  13  1  5
 10  14  2  6
 11  15  3  7
 12  16  4  8

julia> circshift(b, (-1,0))
4×4 Array{Int64,2}:
 2  6  10  14
 3  7  11  15
 4  8  12  16
 1  5   9  13
```

See also [`circshift!`](@ref).
"""
function circshift(a::AbstractArray, shiftamt)
    circshift!(similar(a), a, map(Integer, (shiftamt...,)))
end

# Uses K-B-N summation
function cumsum_kbn{T<:AbstractFloat}(v::AbstractVector{T})
    r = similar(v)
    if isempty(v); return r; end

    inds = indices(v, 1)
    i1 = first(inds)
    s = r[i1] = v[i1]
    c = zero(T)
    for i=i1+1:last(inds)
        vi = v[i]
        t = s + vi
        if abs(s) >= abs(vi)
            c += ((s-t) + vi)
        else
            c += ((vi-t) + s)
        end
        s = t
        r[i] = s+c
    end
    return r
end

# Uses K-B-N summation
# TODO: Needs a separate IndexCartesian method, this is only fast for IndexLinear

"""
    cumsum_kbn(A, [dim::Integer=1])

Cumulative sum along a dimension, using the Kahan-Babuska-Neumaier compensated summation
algorithm for additional accuracy. The dimension defaults to 1.
"""
function cumsum_kbn{T<:AbstractFloat}(A::AbstractArray{T}, axis::Integer=1)
    dimsA = size(A)
    ndimsA = ndims(A)
    axis_size = dimsA[axis]
    axis_stride = 1
    for i = 1:(axis-1)
        axis_stride *= size(A,i)
    end

    if axis_size <= 1
        return A
    end

    B = similar(A)
    C = similar(A)

    for i = 1:length(A)
        if div(i-1, axis_stride) % axis_size == 0
            B[i] = A[i]
            C[i] = zero(T)
        else
            s = B[i-axis_stride]
            Ai = A[i]
            B[i] = t = s + Ai
            if abs(s) >= abs(Ai)
                C[i] = C[i-axis_stride] + ((s-t) + Ai)
            else
                C[i] = C[i-axis_stride] + ((Ai-t) + s)
            end
        end
    end

    return B + C
end

## Other array functions ##

"""
    repmat(A, m::Integer, n::Integer=1)

Construct a matrix by repeating the given matrix (or vector) `m` times in dimension 1 and `n` times in
dimension 2.

```jldoctest
julia> repmat([1, 2, 3], 2)
6-element Array{Int64,1}:
 1
 2
 3
 1
 2
 3

julia> repmat([1, 2, 3], 2, 3)
6×3 Array{Int64,2}:
 1  1  1
 2  2  2
 3  3  3
 1  1  1
 2  2  2
 3  3  3
```
"""
function repmat(a::AbstractVecOrMat, m::Int, n::Int=1)
    o, p = size(a,1), size(a,2)
    b = similar(a, o*m, p*n)
    for j=1:n
        d = (j-1)*p+1
        R = d:d+p-1
        for i=1:m
            c = (i-1)*o+1
            b[c:c+o-1, R] = a
        end
    end
    return b
end

function repmat(a::AbstractVector, m::Int)
    o = length(a)
    b = similar(a, o*m)
    for i=1:m
        c = (i-1)*o+1
        b[c:c+o-1] = a
    end
    return b
end

@inline repmat(a::AbstractVecOrMat, m::Integer, n::Integer=1) = repmat(a, Int(m), Int(n))
@inline repmat(a::AbstractVector, m::Integer) = repmat(a, Int(m))

"""
    repeat(A::AbstractArray; inner=ntuple(x->1, ndims(A)), outer=ntuple(x->1, ndims(A)))

Construct an array by repeating the entries of `A`. The i-th element of `inner` specifies
the number of times that the individual entries of the i-th dimension of `A` should be
repeated. The i-th element of `outer` specifies the number of times that a slice along the
i-th dimension of `A` should be repeated. If `inner` or `outer` are omitted, no repetition
is performed.

```jldoctest
julia> repeat(1:2, inner=2)
4-element Array{Int64,1}:
 1
 1
 2
 2

julia> repeat(1:2, outer=2)
4-element Array{Int64,1}:
 1
 2
 1
 2

julia> repeat([1 2; 3 4], inner=(2, 1), outer=(1, 3))
4×6 Array{Int64,2}:
 1  2  1  2  1  2
 1  2  1  2  1  2
 3  4  3  4  3  4
 3  4  3  4  3  4
```
"""
function repeat(A::AbstractArray;
                inner=ntuple(n->1, Val{ndims(A)}),
                outer=ntuple(n->1, Val{ndims(A)}))
    return _repeat(A, rep_kw2tup(inner), rep_kw2tup(outer))
end

rep_kw2tup(n::Integer) = (n,)
rep_kw2tup(v::AbstractArray{<:Integer}) = (v...)
rep_kw2tup(t::Tuple) = t

rep_shapes(A, i, o) = _rshps((), (), size(A), i, o)

_rshps(shp, shp_i, ::Tuple{}, ::Tuple{}, ::Tuple{}) = (shp, shp_i)
@inline _rshps(shp, shp_i, ::Tuple{}, ::Tuple{}, o) =
    _rshps((shp..., o[1]), (shp_i..., 1), (), (), tail(o))
@inline _rshps(shp, shp_i, ::Tuple{}, i, ::Tuple{}) = (n = i[1];
    _rshps((shp..., n), (shp_i..., n), (), tail(i), ()))
@inline _rshps(shp, shp_i, ::Tuple{}, i, o) = (n = i[1];
    _rshps((shp..., n * o[1]), (shp_i..., n), (), tail(i), tail(o)))
@inline _rshps(shp, shp_i, sz, i, o) = (n = sz[1] * i[1];
    _rshps((shp..., n * o[1]), (shp_i..., n), tail(sz), tail(i), tail(o)))
_rshps(shp, shp_i, sz, ::Tuple{}, ::Tuple{}) =
    (n = length(shp); N = n + length(sz); _reperr("inner", n, N))
_rshps(shp, shp_i, sz, ::Tuple{}, o) =
    (n = length(shp); N = n + length(sz); _reperr("inner", n, N))
_rshps(shp, shp_i, sz, i, ::Tuple{}) =
    (n = length(shp); N = n + length(sz); _reperr("outer", n, N))
_reperr(s, n, N) = throw(ArgumentError("number of " * s * " repetitions " *
    "($n) cannot be less than number of dimensions of input ($N)"))

@propagate_inbounds function _repeat(A::AbstractArray, inner, outer)
    shape, inner_shape = rep_shapes(A, inner, outer)

    R = similar(A, shape)
    if any(iszero, shape)
        return R
    end

    # fill the first inner block
    if all(x -> x == 1, inner)
        R[indices(A)...] = A
    else
        inner_indices = [1:n for n in inner]
        for c in CartesianRange(indices(A))
            for i in 1:ndims(A)
                n = inner[i]
                inner_indices[i] = (1:n) + ((c[i] - 1) * n)
            end
            R[inner_indices...] = A[c]
        end
    end

    # fill the outer blocks along each dimension
    if all(x -> x == 1, outer)
        return R
    end
    src_indices  = [1:n for n in inner_shape]
    dest_indices = copy(src_indices)
    for i in 1:length(outer)
        B = view(R, src_indices...)
        for j in 2:outer[i]
            dest_indices[i] += inner_shape[i]
            R[dest_indices...] = B
        end
        src_indices[i] = dest_indices[i] = 1:shape[i]
    end

    return R
end


text:
const CACHE_FLAG = Ref{Bool}(true)

"""
    enable_cache()

Enable command compile cache. See also [`disable_cache`](@ref).
"""
function enable_cache()
    CACHE_FLAG[] = true
    return
end

"""
    disable_cache()

Disable command compile cache. See also [`enable_cache`](@ref).
"""
function disable_cache()
    CACHE_FLAG[] = false
    return
end

"""
    @cast <function expr>
    @cast <module expr>
    @cast <function name>
    @cast <module name>

Cast a Julia object to a command object. `@cast` will
always execute the given expression, and only create and
register an command object via [`command`](@ref) after
analysing the given expression.

# Example

The most basic way is to use `@cast` on a function as your command
entry.

```julia
@cast function your_command(arg1, arg2::Int)
    # your processing code
end
```

This will create a [`LeafCommand`](@ref) object and register it
to the current module's `CASTED_COMMANDS` constant. You can access
this object via `MODULE_NAME.CASTED_COMMANDS["your_command"]`.

Note this will not create a functional CLI, to create a function
CLI you need to create an entry point, which can be declared by
[`@main`](@ref).
"""
macro cast(ex)
    if CACHE_FLAG[] && iscached()
        return esc(ex)
    else
        return esc(cast_m(__module__, QuoteNode(__source__), ex))
    end
end

function cast_m(m::Module, line::QuoteNode, ex)
    ret = Expr(:block)
    pushmaybe!(ret, create_casted_commands(m))

    if ex isa Symbol
        push!(ret.args, parse_module(m, line, ex))
        return ret
    end

    def = splitdef(ex; throw = false)
    if def === nothing
        push!(ret.args, parse_module(m, line, ex))
        return ret
    end

    push!(ret.args, parse_function(m, line, ex, def))
    return ret
end

function parse_function(m::Module, line::QuoteNode, ex, def)
    haskey(def, :name) || error("command entry cannot be annoymous")
    def[:name] isa Symbol || error("command name should be a Symbol")

    return quote
        $ex
        Core.@__doc__ $(def[:name])
        $(xcall(
            set_cmd!,
            casted_commands(m),
            xcall(command, def[:name], parse_args(def), parse_kwargs(def); line = line),
        ))
    end
end

function parse_module(m::Module, line::QuoteNode, ex::Expr)
    ex.head === :module ||
        throw(Meta.ParseError("invalid expression, can only cast functions or modules"))

    return quote
        $ex
        $(xcall(set_cmd!, casted_commands(m), xcall(command, ex.args[2]; line = line)))
    end
end

function parse_module(m::Module, line::QuoteNode, ex::Symbol)
    return xcall(set_cmd!, casted_commands(m), xcall(command, ex; line = line))
end

function parse_args(def)
    types = Expr(:vect)
    haskey(def, :args) || return types

    # (name, type, require)
    for each in def[:args]
        push!(types.args, Expr(:tuple, to_argument(def, each)...))
    end
    return types
end

function parse_kwargs(def)
    types = Expr(:vect)
    haskey(def, :kwargs) || return types

    # (name, type, isflag)
    for each in def[:kwargs]
        push!(types.args, Expr(:tuple, to_option_or_flag(each)...))
    end
    return types
end

function to_argument(def, ex)
    # (name, type, require, vararg, default_value)
    @smatch ex begin
        ::Symbol => (string(ex), Any, true, false, nothing)
        :($name::$type) => (string(name), wrap_type(def, type), true, false, nothing)
        :($name::$type...) => (string(name), wrap_type(def, type), true, true, nothing)
        Expr(:kw, :($name::$type), value) =>
            (string(name), wrap_type(def, type), false, false, string(value))
        :($name...) => (string(name), Any, true, true, nothing)
        Expr(:kw, name::Symbol, value) => (string(name), Any, false, false, string(value))
        _ => throw(Meta.ParseError("invalid syntax for command line entry: $ex"))
    end
end

function to_option_or_flag(ex)
    @smatch ex begin
        Expr(:kw, name::Symbol, value) => (string(name), Any, false, string(value))
        Expr(:kw, :($name::Bool), false) => (string(name), Bool, true, false)
        Expr(:kw, :($name::$type), value) => (string(name), type, false, string(value))
        Expr(:kw, :($name::Bool), true) =>
            throw(Meta.ParseError("Boolean options must use false as default value, and will be parsed as flags. got $name"))
        ::Symbol || :($name::$type) =>
            throw(Meta.ParseError("options should have default values or make it a positional argument"))
        _ => throw(Meta.ParseError("invalid syntax: $ex"))
    end
end

function wrap_type(def, type)
    if haskey(def, :whereparams)
        return Expr(:where, type, def[:whereparams]...)
    else
        return type
    end
end

function create_casted_commands(m)
    if isdefined(m, :CASTED_COMMANDS)
        return
    else
        return :(const CASTED_COMMANDS = Dict{String,Any}())
    end
end

casted_commands(m) = GlobalRef(m, :CASTED_COMMANDS)

# Entry
"""
    @main <function expr>
    @main [options...]

Create an `EntryCommand` and use it as the entry of the entire CLI.
If you only have one function to cast to a CLI, you can use `@main`
instead of `@cast` so it will create both the command object and
the entry.

If you have declared commands via `@cast`, you can create the entry
via `@main [options...]`, available options are:

- `name`: default is the current module name in lowercase.
- `version`: default is the current project version or `v"0.0.0"`. If
    it's `v"0.0.0"`, the version will not be printed.
- `doc`: a description of the entry command.
"""
macro main(xs...)
    return esc(main_m(__module__, QuoteNode(__source__), xs...))
end

function main_m(m::Module, line::QuoteNode, ex::Expr)
    if CACHE_FLAG[] && iscached()
        return quote
            Core.@__doc__ $ex
            include($(cachefile()[1]))
        end
    end

    ret = Expr(:block)
    def = splitdef(ex; throw = false)
    var_cmd, var_entry = gensym(:cmd), gensym(:entry)
    push!(ret.args, ex)

    if def === nothing
        ex.head === :(=) && return create_entry(m, line, ex)
        ex.head === :module ||
            throw(Meta.ParseError("invalid expression, can only cast functions or modules"))
        cmd = xcall(command, ex.args[2]; line = line)
        push!(ret.args, :($var_cmd = $cmd))
    else
        push!(ret.args, :(Core.@__doc__ $(def[:name])))
        cmd = xcall(command, def[:name], parse_args(def), parse_kwargs(def); line = line)
        push!(ret.args, :($var_cmd = $cmd))
    end

    push!(ret.args, :($var_entry = $(xcall(Types, :EntryCommand, var_cmd; version = get_version(m)))))
    push!(ret.args, precompile_or_exec(m, var_entry))
    return ret
end

function main_m(m::Module, line::QuoteNode, ex::Symbol)
    CACHE_FLAG[] && iscached() && return :(include($(cachefile()[1])))
    var_cmd, var_entry = gensym(:cmd), gensym(:entry)
    quote
        $var_cmd = $(xcall(command, ex; line = line))
        $var_entry = $(xcall(Types, :EntryCommand, var_cmd; line = line))
        $(precompile_or_exec(m, var_entry))
    end
end

function main_m(m::Module, line::QuoteNode, kwargs...)
    CACHE_FLAG[] && iscached() && return :(include($(cachefile()[1])))
    return create_entry(m, line, kwargs...)
end

function create_entry(m::Module, line::QuoteNode, kwargs...)
    configs = Dict{Symbol,Any}(:name => default_name(m), :version => get_version(m), :doc => "")
    for kw in kwargs
        for key in [:name, :version, :doc]
            if kw.args[1] === key
                configs[key] = kw.args[2]
            end
        end
    end

    ret = Expr(:block)

    var_cmd, var_entry = gensym(:cmd), gensym(:entry)
    cmd = xcall(
        Types,
        :NodeCommand,
        configs[:name],
        line,
        configs[:doc],
        :(collect(values($m.CASTED_COMMANDS))),
    )
    entry = xcall(Types, :EntryCommand, var_cmd, configs[:version], line)

    push!(ret.args, :($var_cmd = $cmd))
    push!(ret.args, :($var_entry = $entry))
    push!(ret.args, precompile_or_exec(m, var_entry))
    return ret
end

function precompile_or_exec(m::Module, entry)
    if m == Main && CACHE_FLAG[]
        return quote
            $create_cache($entry)
            include($(cachefile()[1]))
        end
    elseif m == Main
        return quote
            $(xcall(m, :eval, xcall(CodeGen, :codegen, entry)))
            command_main()
        end
    else
        quote
            $(create_casted_commands(m))
            $(xcall(set_cmd!, casted_commands(m), entry, "main"))
            $(xcall(m, :eval, xcall(CodeGen, :codegen, entry)))

            """
                comonicon_build([sysimg=true]; kwargs...)

            Build the CLI manually. This will set system image build to be
            `incremental=true` and `filter_stdlibs=false` to get better compile
            speed locally. For more detailed reference, please
            refer to [Comonicon documentation](https://rogerluo.me/Comonicon.jl/).
            """
            comonicon_build(sysimg = true; kwargs...) =
                $(GlobalRef(Comonicon, :build))($m, sysimg; kwargs...)

            """
                comonicon_install(;kwargs...)

            Install the CLI manually. This will use the default configuration in `Comonicon.toml`,
            if it exists. See also [`comonicon_build`](@ref). For more detailed reference, please
            refer to [Comonicon documentation](https://rogerluo.me/Comonicon.jl/).
            """
            comonicon_install(; kwargs...) = $(GlobalRef(Comonicon, :install))($m; kwargs...)

            """
                comonicon_install_path([quiet=false])

            Install the `PATH` and `FPATH` to your shell configuration file. You can use `comonicon_install_path(true)`,
            to skip interactive prompt.
            For more detailed reference, please refer to [Comonicon documentation](https://rogerluo.me/Comonicon.jl/).
            """
            comonicon_install_path() = $(GlobalRef(Comonicon.BuildTools, :install_env_path))()

            precompile(Tuple{typeof($m.command_main),Array{String,1}})
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE

using FEMBase, MortarContact2D, Test
using MortarContact2D: get_slave_dofs, get_master_dofs

# Original problem:
#
# using JuliaFEM
# using JuliaFEM.Preprocess
# using JuliaFEM.Testing
#
# mesh = Mesh()
# add_node!(mesh, 1, [0.0, 0.0])
# add_node!(mesh, 2, [1.0, 0.0])
# add_node!(mesh, 3, [1.0, 0.5])
# add_node!(mesh, 4, [0.0, 0.5])
# gap = 0.1
# add_node!(mesh, 5, [0.0, 0.5+gap])
# add_node!(mesh, 6, [1.0, 0.5+gap])
# add_node!(mesh, 7, [1.0, 1.0+gap])
# add_node!(mesh, 8, [0.0, 1.0+gap])
# add_element!(mesh, 1, :Quad4, [1, 2, 3, 4])
# add_element!(mesh, 2, :Quad4, [5, 6, 7, 8])
# add_element!(mesh, 3, :Seg2, [1, 2])
# add_element!(mesh, 4, :Seg2, [7, 8])
# add_element!(mesh, 5, :Seg2, [4, 3])
# add_element!(mesh, 6, :Seg2, [6, 5])
# add_element_to_element_set!(mesh, :LOWER, 1)
# add_element_to_element_set!(mesh, :UPPER, 2)
# add_element_to_element_set!(mesh, :LOWER_BOTTOM, 3)
# add_element_to_element_set!(mesh, :UPPER_TOP, 4)
# add_element_to_element_set!(mesh, :LOWER_TOP, 5)
# add_element_to_element_set!(mesh, :UPPER_BOTTOM, 6)
# upper = Problem(Elasticity, "UPPER", 2)
# upper.properties.formulation = :plane_stress
# upper.elements = create_elements(mesh, "UPPER")
# update!(upper.elements, "youngs modulus", 288.0)
# update!(upper.elements, "poissons ratio", 1/3)
# lower = Problem(Elasticity, "LOWER", 2)
# lower.properties.formulation = :plane_stress
# lower.elements = create_elements(mesh, "LOWER")
# update!(lower.elements, "youngs modulus", 288.0)
# update!(lower.elements, "poissons ratio", 1/3)
# bc_upper = Problem(Dirichlet, "UPPER_TOP", 2, "displacement")
# bc_upper.elements = create_elements(mesh, "UPPER_TOP")
# #update!(bc_upper.elements, "displacement 1", -17/90)
# #update!(bc_upper.elements, "displacement 1", -17/90)
# update!(bc_upper.elements, "displacement 1", -0.2)
# update!(bc_upper.elements, "displacement 2", -0.2)
# bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 2, "displacement")
# bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
# update!(bc_lower.elements, "displacement 1", 0.0)
# update!(bc_lower.elements, "displacement 2", 0.0)
# contact = Problem(Contact2D, "LOWER_TO_UPPER", 2, "displacement")
# contact_slave_elements = create_elements(mesh, "LOWER_TOP")
# contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
# add_slave_elements!(contact, contact_slave_elements)
# add_master_elements!(contact, contact_master_elements)
# solver = Solver(Nonlinear)
# push!(solver, upper, lower, bc_upper, bc_lower, contact)
# solver()
# master = first(contact_master_elements)
# slave = first(contact_slave_elements)
# um = master("displacement", (0.0,), 0.0)
# us = slave("displacement", (0.0,), 0.0)
# la = slave("lambda", (0.0,), 0.0)
# info("um = $um, us = $us, la = $la")
# @test isapprox(um, [-0.20, -0.15])
# @test isapprox(us, [0.0, -0.05])
# @test isapprox(la, [0.0, 30.375])

# Equivalent test problem
slave = Element(Seg2, [1, 2])
update!(slave, "geometry", ([0.0, 0.5], [1.0, 0.5]))
update!(slave, "displacement", 0.0 => (zeros(2), zeros(2)))
update!(slave, "lambda", 0.0 => (zeros(2), zeros(2)))
update!(slave, "displacement", 1.0 => ([-0.01875, -0.05], [0.01875, -0.05]))
update!(slave, "lambda", 1.0 => ([0.0, 30.375], [0.0, 30.375]))
master = Element(Seg2, [3, 4])
update!(master, "geometry", ([1.0, 0.6], [0.0, 0.6]))
update!(master, "displacement", 0.0 => (zeros(2), zeros(2)))
update!(master, "displacement", 1.0 => ([-0.18125, -0.15], [-0.21875, -0.15]))
problem = Problem(Contact2D, "LOWER_TO_UPPER", 2, "displacement")
push!(problem.properties.store_fields, "complementarity condition")
push!(problem.properties.store_fields, "active nodes")
push!(problem.properties.store_fields, "inactive nodes")
push!(problem.properties.store_fields, "stick nodes")
push!(problem.properties.store_fields, "slip nodes")
add_slave_elements!(problem, [slave])
add_master_elements!(problem, [master])
um = master("displacement", (0.0,), 1.0)
us = slave("displacement", (0.0,), 1.0)
la = slave("lambda", (0.0,), 1.0)
@test isapprox(um, [-0.20, -0.15])
@test isapprox(us, [0.0, -0.05])
@test isapprox(la, [0.0, 30.375])

assemble!(problem, 0.0)
@test iszero(sparse(problem.assembly.K))
@test iszero(sparse(problem.assembly.C1))
@test iszero(sparse(problem.assembly.C2))
@test iszero(sparse(problem.assembly.D))
@test iszero(sparse(problem.assembly.f))
@test iszero(sparse(problem.assembly.g))

empty!(problem.assembly)
assemble!(problem, 1.0)
@test iszero(sparse(problem.assembly.K))
@test iszero(sparse(problem.assembly.f))
C1 = Matrix(sparse(problem.assembly.C1, 4, 8))
C2 = Matrix(sparse(problem.assembly.C2, 4, 8))
D = Matrix(sparse(problem.assembly.D, 4, 8))
g = Vector(sparse(problem.assembly.g, 4, 1)[:])
C1_expected = [
    0.5  0.0  0.0  0.0   0.0   0.0  -0.5   0.0
    0.0  0.5  0.0  0.0   0.0   0.0   0.0  -0.5
    0.0  0.0  0.5  0.0  -0.5   0.0   0.0   0.0
    0.0  0.0  0.0  0.5   0.0  -0.5   0.0   0.0]
C2_expected = [
    0.0  0.5  0.0  0.0  0.0   0.0  0.0  -0.5
    0.0  0.0  0.0  0.0  0.0   0.0  0.0   0.0
    0.0  0.0  0.0  0.5  0.0  -0.5  0.0   0.0
    0.0  0.0  0.0  0.0  0.0   0.0  0.0   0.0]
D_expected = [
    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
    1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
    0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
    0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0]
g_expected = zeros(4)
@test isapprox(C1, C1_expected)
@test isapprox(C2, C2_expected)
@test isapprox(D, D_expected)
@test isapprox(g, g_expected; atol=1.0e-12)

@test_throws Exception add_elements!(problem, [slave])
@test get_slave_dofs(problem) == [1, 2, 3, 4]
@test get_master_dofs(problem) == [5, 6, 7, 8]


text:
struct TestConfig
    atol::Float64 # absolute tolerance for ...
    rtol::Float64 # relative tolerance for ...
    solve::Bool # optimize and test result
    query::Bool # can get objective function, and constraint functions, and constraint sets
    modify_lhs::Bool # can modify function of a constraint
    duals::Bool # test dual solutions
    infeas_certificates::Bool # check for infeasibility certificates when appropriate
    # The expected "optimal" status returned by the solver. Either
    # MOI.OPTIMAL or MOI.LOCALLY_SOLVED.
    optimal_status::MOI.TerminationStatusCode
    basis::Bool # can get variable and constraint basis status
    function TestConfig(;atol::Float64 = 1e-8, rtol::Float64 = 1e-8, solve::Bool = true, query::Bool = true,
                        modify_lhs::Bool = true, duals::Bool = true, infeas_certificates::Bool = true,
                        optimal_status = MOI.OPTIMAL, basis::Bool = false)
        new(
            atol,
            rtol,
            solve,
            query,
            modify_lhs,
            duals,
            infeas_certificates,
            optimal_status,
            basis
            )
    end
end

"""
    @moitestset setname subsets

Defines a function `setnametest(model, config, exclude)` that runs the tests defined in the dictionary `setnametests`
with the model `model` and config `config` except the tests whose dictionary key is in `exclude`.
If `subsets` is `true` then each test runs in fact multiple tests hence the `exclude` argument is passed
as it can also contains test to be excluded from these subsets of tests.
"""
macro moitestset(setname, subsets=false)
    testname = Symbol(string(setname) * "test")
    testdict = Symbol(string(testname) * "s")
    if subsets
        runtest = :( f(model, config, exclude) )
    else
        runtest = :( f(model, config) )
    end
    esc(:(
        function $testname(model::$MOI.ModelLike, config::$MOI.Test.TestConfig, exclude::Vector{String} = String[])
            for (name,f) in $testdict
                if name in exclude
                    continue
                end
                @testset "$name" begin
                    $runtest
                end
            end
        end
    ))
end


text:
struct Disk
    color
    radius::Vector{Float64}
    parameters::Vector{Complex}
end

"""
    draw_disk(to_draw, origin, dash)
Basics for drawing a disk. The origin is what the disk will be relative to, and dash will be used for circles with outlines
"""
function draw_disk(to_draw::Disk, position::Tuple{Float64,Float64}, dash::Int, relative_size::Float64 = 0.5)
    return compose(context(), circle(position[1], position[2], to_draw.radius[1]*relative_size), fill(to_draw.color))
end

"""
    disk_compose_single(root, origin, leaf_only::Int)
"""
function disk_compose_single(command::Expr, disks::Dict{String, Disk}, center_point, leaf_only::Int)
    tree = []

    center_disk = disks[string(command.args[1])]
    relative_radius = center_disk.radius[1] * 0.5

    if size(command.args)[1] == 1
        return [draw_disk(center_disk, center_point, leaf_only)]
    else

        for i = 2:size(command.args)[1]
            new_center = (0.5 + (0.5*real(center_disk.parameters[i-1])), 0.5 + (0.5*imag(center_disk.parameters[i-1])))
            tree = [tree; compose(context(center_point[1]-relative_radius, center_point[2]-relative_radius, 2*relative_radius, 2*relative_radius),
                disk_compose_single(command.args[i], disks, new_center, leaf_only)...)]
        end

        if leaf_only != 2
            tree = [tree; [draw_disk(center_disk, center_point, leaf_only)]]
        end
    end
    return tree
end

"""
    disk_compose_single_base(root, leaf_only::Int)
"""
function disk_compose_single_base(command::Expr, disks::Dict{String,Disk}, leaf_only::Int, )
    tree = []

    center_point = (0.5, 0.5)
    center_disk = disks[string(command.args[1])]
    relative_radius = center_disk.radius[1] * 0.5

    for i = 2:size(command.args)[1]
        new_center = (0.5 + (0.5*real(center_disk.parameters[i-1])), 0.5 + (0.5*imag(center_disk.parameters[i-1])))
        tree = [tree; compose(context(center_point[1]-relative_radius, center_point[2]-relative_radius, 2*relative_radius, 2*relative_radius),
            disk_compose_single(command.args[i], disks, new_center, leaf_only)...)]
    end

    return (context(),
    tree...,
    compose(context(), circle(center_point[1], center_point[2], center_disk.radius[1]*0.5), fill(center_disk.color)),
    compose(context(), rectangle()), fill(colorant"rgba(0,0,0,0)"))
end

max_stars_repo_path: src/Assisting Files/Disk.jl

max_issues_repo_path: src/Assisting Files/Disk.jl

max_forks_repo_path: src/Assisting Files/Disk.jl


text:
using DiffEqFlux, CUDA, OrdinaryDiffEq, StochasticDiffEq, Test

CUDA.allowscalar(false)

mp = Chain(Dense(2,2)) |> gpu
x = Float32[2.; 0.] |> gpu
xs = Float32.(hcat([0.; 0.], [1.; 0.], [2.; 0.])) |> gpu
tspan = (0.0f0,25.0f0)
dudt = Chain(Dense(2,50,tanh),Dense(50,2)) |> gpu
dudt_tracker = Chain(Dense(2,50,CUDA.tanh),Dense(50,2)) |> gpu

NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)(x)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)(x)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())(x)

NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)(xs)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1)(xs)
NeuralODE(dudt,tspan,Tsit5(),saveat=0.1,sensealg=TrackerAdjoint())(xs)

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt_tracker,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=TrackerAdjoint())
grads = Zygote.gradient(()->sum(Array(node(x))),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false,sensealg=BacksolveAdjoint())
grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
@test ! iszero(grads[x])
@test ! iszero(grads[node.p])

grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
@test ! iszero(grads[xs])
@test ! iszero(grads[node.p])

# Adjoint
@testset "adjoint mode" begin
    node = NeuralODE(dudt,tspan,Tsit5(),save_everystep=false,save_start=false)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=0.0:0.1:10.0)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])

    node = NeuralODE(dudt,tspan,Tsit5(),saveat=1f-1)
    grads = Zygote.gradient(()->sum(node(x)),Flux.params(x,node))
    @test ! iszero(grads[x])
    @test ! iszero(grads[node.p])

    grads = Zygote.gradient(()->sum(node(xs)),Flux.params(xs,node))
    @test ! iszero(grads[xs])
    @test ! iszero(grads[node.p])
end

NeuralDSDE(dudt,mp,(0.0f0,2.0f0),SOSRI(),saveat=0.0:0.1:2.0)(x)
sode = NeuralDSDE(dudt,mp,(0.0f0,2.0f0),SOSRI(),saveat=Float32.(0.0:0.1:2.0),dt=1f-1)
@test_broken grads = Zygote.gradient(()->sum(sode(x)),Flux.params(x,sode)) isa Tuple

@test ! iszero(grads[x])
@test ! iszero(grads[sode.p])

grads = Zygote.gradient(()->sum(sode(xs)),Flux.params(xs,sode))
@test ! iszero(grads[xs])
@test ! iszero(grads[sode.p])


text:
const internalFrac = 0.5
const internalV = 623.0
const fgMgConv = 135.2
const Nspecies = 55

# Mark surface species
# Variables with C on end are only for one compartment
surface = vcat(ones(27), zeros(28))
pYc = vcat(zeros(4), ones(2), zeros(4), ones(2), zeros(4), ones(2), ones(9))
pY = vcat(pYc, pYc * internalFrac, 0)
ligPiece = [0, 1, 1, 2, 1, 2]
hetPiece = [1, 1, 2]
boundLigC = vcat(ligPiece, ligPiece, ligPiece, hetPiece, hetPiece, hetPiece)
boundLig = vcat(boundLigC, boundLigC * internalFrac, 0)
totalPiece = [1, 1, 1, 1, 2, 2]
totalC = vcat(totalPiece, totalPiece, totalPiece, 2 * ones(9))
total = vcat(totalC, totalC * internalFrac, 0)
recpSpecificC = [
    vcat(ones(6), zeros(12), ones(3), zeros(3), ones(3)),  # AXL
    vcat(zeros(6), ones(6), zeros(6), ones(6), zeros(3)),  # MerTK
    vcat(zeros(6), zeros(6), ones(6), zeros(3), ones(6)),
] # Tyro3
recpSpecific = [
    vcat(recpSpecificC[1], recpSpecificC[1] * internalFrac, 0),  # AXL
    vcat(recpSpecificC[2], recpSpecificC[2] * internalFrac, 0),  # MerTK
    vcat(recpSpecificC[3], recpSpecificC[3] * internalFrac, 0),
] # Tyro3

" This makes a rate for a receptor. "
function paramTAMrate(params::AbstractVector{T}) where {T}
    fBnd = 0.6
    AXL = TAMrates{T}([1.2, 0.042, fBnd, fBnd * params[4]], zeros(6), params[1], 1.50) # From Kariolis et al
    MerTK = TAMrates{T}([fBnd, fBnd * params[5], fBnd, fBnd * params[6]], zeros(6), params[2], 0.0)
    Tyro3 = TAMrates{T}([fBnd, fBnd * params[7], fBnd, fBnd * params[8]], zeros(6), params[3], 0.0)

    TAM = TAMsType{T}(AXL, MerTK, Tyro3)

    hetR = hetRates{T}(zeros(10), 0.0, 0.0)
    hetRs = hetRType{T}(hetR, deepcopy(hetR), deepcopy(hetR))

    TAM.Axl.xRev[5] = 0.0144 # From Kariolis et al

    return TAM, hetRs
end

" Setup the parameters for the full TAM receptor model. "
function param(params::Vector{T})::Rates{T} where {T}
    @assert all(params .>= 0.0)
    @assert params[3] <= 1.0

    TAM, hetRs = paramTAMrate(view(params, 8:15))

    out = Rates{T}(TAM, params[1], params[2], params[3], params[4], params[5], params[6], params[7], hetRs)

    return detailedBalance(out)
end

" Tracks the hetero-receptor interactions between two receptors. "
function het_module(Rone, Rtwo, dRone, dRtwo, hetR, hetDim, dhetDim, tr, Gas, dLi)
    dRr1 = tr.xFwd * Rone[3] * Rtwo[3] - hetR.xRev[1] * hetDim[3]
    dRr2 = tr.xFwd * Rone[1] * Rtwo[4] - hetR.xRev[2] * hetDim[3]
    dRr3 = tr.xFwd * Rone[2] * Rtwo[2] - hetR.xRev[3] * hetDim[3]
    dRr4 = tr.xFwd * Rone[4] * Rtwo[1] - hetR.xRev[4] * hetDim[3]
    dRr5 = tr.xFwd * Rone[3] * Rtwo[1] - hetR.xRev[5] * hetDim[2]
    dRr6 = tr.xFwd * Rone[1] * Rtwo[2] - hetR.xRev[6] * hetDim[2]
    dRr7 = tr.xFwd * Rone[2] * Rtwo[1] - hetR.xRev[7] * hetDim[1]
    dRr8 = tr.xFwd * Rone[1] * Rtwo[3] - hetR.xRev[8] * hetDim[1]
    dRr9 = hetR.xFwd15 * Gas * hetDim[1] - hetR.xRev[9] * hetDim[3]
    dRr10 = hetR.xFwd16 * Gas * hetDim[2] - hetR.xRev[10] * hetDim[3]

    dRone[1] += -dRr2 - dRr6 - dRr8
    dRone[2] += -dRr3 - dRr7
    dRone[3] += -dRr1 - dRr5
    dRone[4] += -dRr4

    dRtwo[1] += -dRr4 - dRr5 - dRr7
    dRtwo[2] += -dRr3 - dRr6
    dRtwo[3] += -dRr1 - dRr8
    dRtwo[4] += -dRr2

    dhetDim[1] += dRr7 + dRr8 - dRr9 # AMD1
    dhetDim[2] += dRr5 + dRr6 - dRr10 # MAD1
    dhetDim[3] += dRr1 + dRr2 + dRr3 + dRr4 + dRr9 + dRr10 # AMD2

    if !(dLi isa Nothing)
        dLi[1] += -dRr9 - dRr10
    end

    nothing
end

" Handles trafficking of receptor and ligand. "
function trafFunc(du, u, r, pYcIn)
    for ii = 1:length(pYcIn)
        if pYcIn[ii] == 1
            fElse = 1.0
            intRate = r.pYinternalize
        else
            fElse = r.fElse
            intRate = r.internalize
        end

        du[ii] += -u[ii] * intRate + r.kRec * (1.0 - fElse) * u[ii + length(pYcIn)] * internalFrac # Endocytosis, recycling
        du[ii + length(pYcIn)] +=
            u[ii] * intRate / internalFrac - r.kRec * (1.0 - fElse) * u[ii + length(pYcIn)] - r.kDeg * fElse * u[ii + length(pYcIn)] # ", degradation
    end
end

function compartmentReact(u, du, Gas, dLi, r)
    for (ii, aa) in enumerate((1, 7, 13))
        dR = view(du, aa:(aa + 5))
        R = view(u, aa:(aa + 5))

        dRr1 = r.TAMs[ii].binding[1] * R[1] * Gas - r.TAMs[ii].binding[2] * R[2]
        dRr2 = r.TAMs[ii].binding[3] * R[1] * Gas - r.TAMs[ii].binding[4] * R[3]
        dRr3 = r.TAMs[ii].binding[3] * R[2] * Gas - r.TAMs[ii].binding[4] * R[4]
        dRr4 = r.TAMs[ii].binding[1] * R[3] * Gas - r.TAMs[ii].binding[2] * R[4]
        dRr5 = r.xFwd * R[1] * R[2] - r.TAMs[ii].xRev[1] * R[5]
        dRr6 = r.xFwd * R[1] * R[3] - r.TAMs[ii].xRev[2] * R[5]
        dRr7 = r.xFwd * R[1] * R[4] - r.TAMs[ii].xRev[3] * R[6]
        dRr8 = r.xFwd * R[2] * R[2] - r.TAMs[ii].xRev[4] * R[6]
        dRr9 = r.xFwd * R[3] * R[3] - r.TAMs[ii].xRev[5] * R[6]
        dRr10 = r.TAMs[ii].xFwd6 * Gas * R[5] - r.TAMs[ii].xRev[6] * R[6]

        dR[1] += -dRr7 - dRr6 - dRr5 - dRr1 - dRr2 # AXL
        dR[2] += -2 * (dRr8) - dRr5 + dRr1 - dRr3 # AXLgas1
        dR[3] += -2 * (dRr9) - dRr6 + dRr2 - dRr4 # AXLgas2
        dR[4] += -dRr7 + dRr3 + dRr4 # AXLgas12
        dR[5] += -dRr10 + dRr6 + dRr5 # AXLdimer1
        dR[6] += dRr10 + dRr9 + dRr8 + dRr7 # AXLdimer2

        if !(dLi isa Nothing)
            dLi[1] += -dRr10 - dRr1 - dRr2 - dRr3 - dRr4
        end

        dR[1] += r.TAMs[ii].expression
    end

    het_module(view(u, 1:6), view(u, 7:12), view(du, 1:6), view(du, 7:12), r.hetR.AM, view(u, 19:21), view(du, 19:21), r, Gas, dLi)
    het_module(view(u, 7:12), view(u, 13:18), view(du, 7:12), view(du, 13:18), r.hetR.MT, view(u, 22:24), view(du, 22:24), r, Gas, dLi)
    het_module(view(u, 1:6), view(u, 13:18), view(du, 1:6), view(du, 13:18), r.hetR.AT, view(u, 25:27), view(du, 25:27), r, Gas, dLi)
end

function TAMreact(du, u, r::Rates, t)
    fill!(du, 0.0)

    compartmentReact(view(u, 1:27), view(du, 1:27), r.gasCur, nothing, r)
    compartmentReact(view(u, 28:54), view(du, 28:54), u[end] / internalV, view(du, Nspecies), r)

    trafFunc(du, u, r, pYc)

    du[end] = -r.kDeg * u[end] # Gas6 degradation

    nothing
end

function detailedBalance(out)
    for T in out.TAMs
        KD1 = T.binding[2] / T.binding[1]
        KD2 = T.binding[4] / T.binding[3]

        if T.binding[4] <= T.binding[2]
            T.xRev[1] = T.binding[4]
            T.xRev[2] = T.xRev[1] * KD1 / KD2
        else
            T.xRev[2] = T.binding[2]
            T.xRev[1] = T.xRev[2] * KD2 / KD1
        end
    end

    for ii = 2:3
        out.TAMs[ii].xFwd6 = max(out.TAMs[ii].binding[1], out.TAMs[ii].binding[3])
        out.TAMs[ii].xRev[6] =
            out.TAMs[1].xRev[6] * out.TAMs[ii].binding[2] * out.TAMs[ii].binding[4] / out.TAMs[1].binding[4] / out.TAMs[1].binding[2]
    end

    for T in out.TAMs
        KD1 = T.binding[2] / T.binding[1]
        KD2 = T.binding[4] / T.binding[3]

        T.xRev[5] = T.xRev[6] * KD1 * T.xRev[1] / KD2 / KD2 / T.xFwd6
        T.xRev[4] = T.xRev[5] * KD2 * KD2 / KD1 / KD1
        T.xRev[3] = T.xRev[4] * KD1 / KD2
    end

    # ligand binding to the one ligand dimer
    for ii = 1:3
        x = [:Axl, :MerTK, :Axl][ii]
        y = [:MerTK, :Tyro3, :Tyro3][ii]

        KD11 = out.TAMs[x].binding[2] / out.TAMs[x].binding[1]
        KD12 = out.TAMs[y].binding[2] / out.TAMs[y].binding[1]
        KD21 = out.TAMs[x].binding[4] / out.TAMs[x].binding[3]
        KD22 = out.TAMs[y].binding[4] / out.TAMs[y].binding[3]

        out.hetR[ii].xFwd15 = max(out.TAMs[y].binding[1], out.TAMs[x].binding[3])
        out.hetR[ii].xFwd16 = max(out.TAMs[x].binding[1], out.TAMs[y].binding[3])

        out.hetR[ii].xRev[1] = out.TAMs[x].xRev[5] * KD12 / KD11
        out.hetR[ii].xRev[2] = out.hetR[ii].xRev[1] * KD21 / KD12
        out.hetR[ii].xRev[3] = KD22 * KD21 / KD11 / KD12 * out.hetR[ii].xRev[1]
        out.hetR[ii].xRev[4] = out.hetR[ii].xRev[1] * KD22 / KD11

        if out.TAMs[x].binding[2] <= out.TAMs[x].binding[4]
            out.hetR[ii].xRev[8] = out.TAMs[x].binding[2]
            out.hetR[ii].xRev[6] = out.hetR[ii].xRev[8] * KD11 / KD21
        else
            out.hetR[ii].xRev[6] = out.TAMs[x].binding[4]
            out.hetR[ii].xRev[8] = out.hetR[ii].xRev[6] * KD21 / KD11
        end

        if out.TAMs[y].binding[2] <= out.TAMs[y].binding[4]
            out.hetR[ii].xRev[5] = out.TAMs[y].binding[2]
            out.hetR[ii].xRev[7] = out.hetR[ii].xRev[5] * KD22 / KD12
        else
            out.hetR[ii].xRev[7] = out.TAMs[y].binding[4]
            out.hetR[ii].xRev[5] = out.hetR[ii].xRev[7] * KD12 / KD22
        end

        out.hetR[ii].xRev[9] = out.hetR[ii].xRev[1] * KD21 / KD12
        out.hetR[ii].xRev[10] = out.hetR[ii].xRev[3] * KD11 / KD22
    end

    return out
end

function swapIgs(out::Rates{TT})::Rates{TT} where {TT <: Real}
    out = detailedBalance(deepcopy(out))

    for T in out.TAMs
        T.binding[1], T.binding[3] = T.binding[3], T.binding[1]
        T.binding[2], T.binding[4] = T.binding[4], T.binding[2]
        T.xRev[4], T.xRev[5] = T.xRev[5], T.xRev[4]
    end

    return detailedBalance(out)
end


text:
f = x -> x^3
f_dash = x -> 3 * x^2
partition = collect(-1.0:1.0:1.0)

@testset "sampling tests with x^3 for error tolerance" begin
    PR.silence()
    p = deepcopy(partition)
    for i = 1:10
        lb = -1.0
        ub = 1.0
        λ = rand()
        x_val = λ * lb + (1 - λ) * ub
        m = Model(cbc_optimizer)
        @variable(m, -1.0 <= x <= 1.0)
        @variable(m, y)
        construct_univariate_relaxation!(m, f, x, y, p, true, error_tolerance = 1e-2)
        @constraint(m, x == x_val)

        # solve for min y
        @objective(m, Min, y)
        optimize!(m)
        y_min = objective_value(m)

        # solve for max y
        @objective(m, Max, y)
        optimize!(m)
        y_max = objective_value(m)
        @test abs(y_max - y_min) <= 1e-2
    end
end

@testset "sampling tests with x^3 for binary variable budget" begin
    num_vars = 10
    tol = 1e-5
    for l in [0.1, 0.25, 0.5, 1.0]
        base_partition = collect(-1.0:l:1.0)
        num_base = length(base_partition) - 1
        m = Model(cbc_optimizer)
        @variable(m, -1.0 <= x <= 1.0)
        @variable(m, y)
        formulation_info = construct_univariate_relaxation!(
            m,
            f,
            x,
            y,
            base_partition,
            true,
            error_tolerance = tol,
            num_additional_partitions = num_vars,
        )
        @test length(formulation_info.variables[:z]) == num_base + num_vars
    end
end

@testset "sampling tests with x^3 for length tolerance" begin
    err_tol = 1e-7
    len_tol = 0.1
    p = deepcopy(partition)
    m = Model(cbc_optimizer)
    @variable(m, -1.0 <= x <= 1.0)
    @variable(m, y)
    construct_univariate_relaxation!(
        m,
        f,
        x,
        y,
        p,
        true,
        error_tolerance = err_tol,
        length_tolerance = len_tol,
    )
    for i = 1:length(p)-1
        @test p[i+1] - p[i] >= len_tol
    end
end

@testset "sampling tests with x^3 for derivative tolerance" begin
    err_tol = 1e-7
    der_tol = 0.1
    p = deepcopy(partition)
    m = Model(cbc_optimizer)
    @variable(m, -1.0 <= x <= 1.0)
    @variable(m, y)
    construct_univariate_relaxation!(
        m,
        f,
        x,
        y,
        p,
        true,
        error_tolerance = err_tol,
        derivative_tolerance = der_tol,
    )

    for i = 1:length(p)-1
        d = f_dash(p[i])
        d_next = f_dash(p[i+1])
        @test abs(d_next - d) >= der_tol
    end
end


text:
struct ROCCompilerParams <: AbstractCompilerParams
    device::RuntimeDevice{HSAAgent}
    global_hooks::NamedTuple
end

ROCCompilerJob = CompilerJob{GCNCompilerTarget,ROCCompilerParams}

GPUCompiler.runtime_module(::ROCCompilerJob) = AMDGPU

# filter out functions from device libs
GPUCompiler.isintrinsic(job::ROCCompilerJob, fn::String) =
    invoke(GPUCompiler.isintrinsic,
           Tuple{CompilerJob{GCNCompilerTarget}, typeof(fn)},
           job, fn) ||
    startswith(fn, "rocm")

function GPUCompiler.process_module!(job::ROCCompilerJob, mod::LLVM.Module)
    invoke(GPUCompiler.process_module!,
           Tuple{CompilerJob{GCNCompilerTarget}, typeof(mod)},
           job, mod)
    # Run this early (before optimization) to ensure we link OCKL
    emit_exception_user!(mod)
end
function GPUCompiler.finish_module!(job::ROCCompilerJob, mod::LLVM.Module)
    invoke(GPUCompiler.finish_module!,
           Tuple{CompilerJob{GCNCompilerTarget}, typeof(mod)},
           job, mod)
    delete_exception_user!(mod)
end

function GPUCompiler.link_libraries!(job::ROCCompilerJob, mod::LLVM.Module,
                                     undefined_fns::Vector{String})
    invoke(GPUCompiler.link_libraries!,
           Tuple{CompilerJob{GCNCompilerTarget}, typeof(mod), typeof(undefined_fns)},
           job, mod, undefined_fns)

    link_device_libs!(job.target, mod)
end

GPUCompiler.ci_cache(::ROCCompilerJob) = ci_cache

GPUCompiler.method_table(::ROCCompilerJob) = method_table


text:
using DataFrames:rename!

"""
    fastby(fn, b[, v])

Group by `b` then apply `fn` to `v` once grouped. If `v` is not provided then
`fn` is applied `b`
"""
fastby!(fn, b::AbstractVector) = begin
    sort!(b, alg = RadixSort)
    _contiguousby_rle(fn, b, b)
end

fastby(fn, b) = fastby!(fn, copy(b))

"""
Fast Group By algorithm
"""
fastby(fn::Function, byvec, valvec) =  length(byvec) == length(valvec) == 0 ? throw(error("length of byvec and valvec can not be 0")) : fastby!(fn, copy(byvec), copy(valvec))

# fastby(fn::Function, byvec, valvec) = length(byvec) == length(valvec) == 0 ? throw(error("length of byvec and valvec can not be 0")) : fastby!(fn, copy(byvec), copy(valvec))

"""
    fast(fn, df, bycol; name = [])

"""
fastby(fn::Function, df::AbstractDataFrame, bycol::Symbol; name = Symbol(string(fn)*string(bycol))) = begin
    res = fastby(fn, df[!, bycol])

    resdf = DataFrame(res)
    rename!(resdf, :x1=>bycol, :x2=>name)
end

function fastby(fn::Function, df::DF, bycol::Symbol, valcol::Symbol) where DF <: AbstractDataFrame
    res_vec = DataFrame(fastby!(fn, copy(df[!, bycol]), copy(df[!, valcol])))
    rename!(res_vec, :x1 => bycol, :x2 => :V1)
    res_vec
end

# fastby(fn::NTuple{N, Function}, df::AbstractDataFrame, bycol::Symbol, valcol::NTuple{N,Symbol}) where N =
#     DataFrame(fastby(fn, df[bycol], ((df[vc] for vc in valcol)...)) |> collect, vcat(bycol, valcol...))

function fastby(fn::Function, x::Vector{Bool}, y)
    # TODO: fast path for sum and mean
    ([true, false], [fn(@view(y[x])), fn(@view(y[.!x]))])
end

function fastby!(fn::Function, byvec::AbstractVector{T}, valvec::AbstractVector{S}) where {T, S}
    length(byvec) == length(valvec) || throw(DimensionMismatch())
    outType = typeof(fn(valvec[1:1]))
    if issorted(byvec)
        h = _contiguousby_vec(fn, byvec, valvec)
    else
        h = _fastby!(fn, byvec, valvec)
    end
    return h
end

"""
Internal: single-function fastby, one by, one val
"""
function _fastby!(fn::Function, byvec::AbstractVector{T}, valvec::AbstractVector{S}) where {T <: Union{BaseRadixSortSafeTypes, String}, S}
    # l = length(byvec)
    #grouptwo!(byvec, valvec)
    SortingLab.sorttwo!(byvec, valvec)
    #return _contiguousby(fn, byvec, valvec)
    return _contiguousby_vec(fn, byvec, valvec)
end

function _fastby!(fn::Function, byvec::AbstractVector{Bool}, valvec::AbstractVector{S}) where S
    grouptwo!(byvec, valvec)
    return _contiguousby_vec(fn, byvec, valvec)
end

function fastby(fn::Function, df::DataFrame, byvec::Union{AbstractVector{Symbol}, NTuple{N, Symbol}}, valsymbol::Symbol) where N
    indexes = fcollect(size(df,1))
    for bv in reverse(byvec)
        cdfbv = df[bv][indexes]
        grouptwo!(cdfbv, indexes)
    end

    # by now all the groups are done
    dfiv = df[indexes, byvec]
    lo = 1
    lastrow = dfiv[lo,:]
    valvec = df[valsymbol]

    # df1 = DataFrame()
    for i in 2:size(dfiv,1)
        if dfiv[i,:] != lastrow
            lastrow[valsymbol] = fn(valvec[lo:i-1])

            # df1 = lastrow
            lo = i
            lastrow = dfiv[lo,byvec]
            break
        end
    end

    for i in lo:size(dfiv,1)
        if dfiv[i,:] != lastrow
            lastrow[valsymbol] = fn(valvec[lo:i-1])
            # df1 = vcat(df1, lastrow)
            lo = i
            lastrow = dfiv[lo,byvec]
        end
    end

    lastrow[valsymbol] = fn(valvec[lo:size(dfiv,1)])
    lastrow
    # df1 = vcat(df1, lastrow)
    # df1
end

"""
Internal multi-function fastby, one by, one val
"""
# function _fastby!(fn::Union{Vector{Function},NTuple{N, Function}}, byvec::AbstractVector{T}, valvec::AbstractVector{S}) where {N, T <: BaseRadixSortSafeTypes, S}
#     l = length(byvec)
#     grouptwo!(byvec, valvec)
#     lastby = byvec[1]
#
#     res = Dict{T}()
#
#     j = 1
#
#     for i = 2:l
#         @inbounds byval = byvec[i]
#         if byval != lastby
#             viewvalvec = @view valvec[j:i-1]
#             @inbounds res[lastby] = ((fn1(viewvalvec) for fn1 in fn)...)
#             j = i
#             @inbounds lastby = byvec[i]
#         end
#     end
#
#     viewvalvec = @view valvec[j:l]
#     @inbounds res[byvec[l]] = ((fn1(viewvalvec) for fn1 in fn)...)
#     return res
# end

"""
Internal multi-function fastby, 1 categorical by, one val
"""
# function _fastby!(fn::Vector{Function}, byvec::CategoricalVector, valvec::AbstractVector{S}) where {S}
#     l = length(byvec)
#     grouptwo!(byvec, valvec)
#     lastby = byvec[1]
#
#     res = Dict{T}()
#
#     j = 1
#
#     for i = 2:l
#         @inbounds byval = byvec[i]
#         if byval != lastby
#             viewvalvec = @view valvec[j:i-1]
#             @inbounds res[lastby] = ((fn1(viewvalvec) for fn1 in fn)...)
#             j = i
#             @inbounds lastby = byvec[i]
#         end
#     end
#
#     viewvalvec = @view valvec[j:l]
#     @inbounds res[byvec[l]] = ((fn1(viewvalvec) for fn1 in fn)...)
#     return res
# end


text:
using thesis, PGFPlotsX
setup_pgfplotsx()

ts = LinRange(thesis.example_trange..., 21)
tts = LinRange(thesis.example_trange..., 201)

## Digital neuron
num_bits = 4
thresholds_dig = LinRange(thesis.example_yrange..., 2^num_bits+1)
levels_dig = 0.5*(thresholds_dig[1:end-1]+thresholds_dig[2:end])
ys_dig = thesis.example_input.(ts)
ys_bits_dig = thesis.discretize.(ys_dig, Ref(thresholds_dig), Ref(digits.(eachindex(thresholds_dig).-1, base=2, pad=num_bits)))
ys_num_dig = thesis.discretize.(ys_dig, Ref(thresholds_dig), Ref(levels_dig))
ys_trans_dig = thesis.discretize.(ys_dig, Ref(levels_dig), Ref(string.(0:2^num_bits,base=2^num_bits)))
y_rec_dig(t) = ys_num_dig[searchsortedlast(ts, t)]
τ = minimum(diff(ts))

β = 2.0
δ = 0.3
## LIF neuron
sol_lif,spikes_lif=lif_pair_get_solution(thesis.example_trange, thesis.example_input;α=β,β=β,δ=δ)
y_lif=thesis.example_input_integral.(spikes_lif)
y_rec_lif = t->sol_lif(t,idxs=2)
y_filtered = t->sol_lif(t,idxs=3)

δ = 0.15
## LNP neuron
sol_lnp,spikes_lnp=lnp_pair_get_solution(thesis.example_trange, thesis.example_input ;β=β,δ=δ)
y_lnp=thesis.example_input_integral.(spikes_lnp)
y_rec_lnp = t->sol_lnp(t,idxs=1)
y_filtered = t->sol_lnp(t,idxs=2)

## Build the monstrosity

@pgf gp = TikzPicture(
    GroupPlot(
    {
        group_style={group_size="3 by 4"},
        major_tick_style={draw="none"},
    },
    # input for digital
    {title="digital encoding", height="4cm",width="6cm", grid = "minor", ylabel=raw"input $s(t)$", xticklabel="\\empty", minor_xtick=ts, minor_ytick=levels_dig, ymin=thesis.example_yrange[1], ymax=thesis.example_yrange[2], xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick},
        Table(tts, thesis.example_input.(tts))
    ),
    (Plot(
        {no_markers, color="Set2-B", very_thick},
        Coordinates([(t,0),(t,thesis.example_input(t))])
        ) for (i,t) ∈ enumerate(ts))...,
    PlotInc(
        {only_marks, color="Set2-B"},
        Table(ts, ys_dig)
    ),
    PlotInc(
        {only_marks, mark="+", color="black"},
        Table(ts, ys_num_dig)
    ),
    [raw"\coordinate (bot11) at (axis cs:5,0);"],
    # input for lif
    {title="lIF encoding", height="4cm",width="6cm", grid = "minor", yminorgrids="false",xticklabel="\\empty", minor_xtick=spikes_lif, ymin=thesis.example_yrange[1], ymax=thesis.example_yrange[2], xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick, fill="Set2-B", fill_opacity=0.5},
        Table([thesis.example_trange[1];tts;thesis.example_trange[2]], [thesis.example_yrange[1];thesis.example_input.(tts);thesis.example_yrange[1]])
    ),
    # (Plot(
    #     {no_markers, color="black", very_thick},
    #     Coordinates([(t,0),(t,thesis.example_input(t))])
    #     ) for (i,t) ∈ enumerate(ts))...,
    [raw"\coordinate (bot12) at (axis cs:5,0);"],
    # input for lnp
    {title="LNP encoding", height="4cm",width="6cm", grid = "minor", yminorgrids="false",xticklabel="\\empty", minor_xtick=spikes_lnp, ymin=thesis.example_yrange[1], ymax=thesis.example_yrange[2], xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick, fill="Set2-B", fill_opacity=0.5},
        Table([thesis.example_trange[1];tts;thesis.example_trange[2]], [thesis.example_yrange[1];thesis.example_input.(tts);thesis.example_yrange[1]])
    ),
    # (Plot(
    #     {no_markers, color="black", very_thick},
    #     Coordinates([(t,0),(t,thesis.example_input(t))])
    #     ) for (i,t) ∈ enumerate(ts))...,
    [raw"\coordinate (bot13) at (axis cs:5,0);"],
    # mechanism for digital
    {height="4cm",width="6cm", xticklabel="\\empty", "group/empty plot"},
    # mechanism for lif
    {xtick_pos="bottom", yticklabel_pos="left", ytick_pos="right", tick_align="outside", color="black", minor_tick_length="2mm", height="4cm",width="6cm", grid = "minor", ylabel=raw"{$\int_{0}^t s(t)dt$}", xticklabel="\\empty", minor_xtick=spikes_lif, minor_ytick=y_lif, ymin=0, xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick, color="Set2-B"},
        Table(tts, thesis.example_input_integral.(tts))
    ),
    Plot(
        {only_marks, mark="+", color="black"},
        Table(spikes_lif, y_lif)
    ),
    [raw"\coordinate (top22) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot22) at (axis cs:5,0);"],
    # mechanism for lnp
    {xtick_pos="bottom", yticklabel_pos="left", ytick_pos="right", tick_align="outside", minor_tick_length="2mm", height="4cm",width="6cm", grid = "minor", xticklabel="\\empty", minor_xtick=spikes_lnp, minor_ytick=y_lnp, ymin=0, xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick, color="Set2-B"},
        Table(tts, thesis.example_input_integral.(tts))
    ),
    Plot(
        {only_marks, mark="+", color="black"},
        Table(spikes_lnp, y_lnp)
    ),
    [raw"\coordinate (top23) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot23) at (axis cs:5,0);"],
    # message for digital
    {height="2cm",width="6cm", grid = "minor", ylabel="message", xticklabel="\\empty", minor_xtick=ts, ytick="\\empty", ymin=0, ymax=1, xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        Table([],[])
    ),
    ["\\node[] at (axis cs:$t,0.5) {\\strut $y};" for (t,y) ∈ zip(0.5*(ts[1:end-1].+ts[2:end]), ys_trans_dig)],
    ["\\coordinate (label31) at (\$(axis cs:10,0)+(0,-4pt)\$);", raw"\coordinate (top31) at (axis cs:5,1);\coordinate (bot31) at (axis cs:5,0);"],
    # message for lif
    {height="2cm",width="6cm", grid = "minor", ylabel="spikes", xticklabel="\\empty", xtick="\\empty", ytick="\\empty", ymin=0, ymax=1, xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    (
        Plot(
            {no_markers, color="black", very_thick},
            Coordinates([(t,0),(t,1)])
        ) for (i,t) ∈ enumerate(spikes_lif)
    )...,
    ["\\coordinate (label32) at (\$(axis cs:10,0)+(0,-4pt)\$);", raw"\coordinate (top32) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot32) at (axis cs:5,0);"],
    # message for lnp
    {height="2cm",width="6cm", grid = "minor", xticklabel="\\empty", xtick="\\empty", ytick="\\empty", ymin=0, ymax=1, xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    (
        Plot(
            {no_markers, color="black", very_thick},
            Coordinates([(t,0),(t,1)])
        ) for (i,t) ∈ enumerate(spikes_lnp)
    )...,
    ["\\coordinate (label33) at (\$(axis cs:10,0)+(0,-4pt)\$);", raw"\coordinate (top33) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot33) at (axis cs:5,0);"],
    # reconstruction for digital
    {height="4cm",width="6cm", grid = "minor", xlabel=raw"time $t$", ylabel="reconstruction", minor_xtick=ts, minor_ytick=levels_dig, ymin=thesis.example_yrange[1], ymax=thesis.example_yrange[2], xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick},
        Table(tts, thesis.example_input.(tts))
    ),
    Plot(
        {no_marks, very_thick, color="black"},
        Table(tts, y_rec_dig.(tts))
    ),
    PlotInc(
        {no_marks, color="Set2-B", very_thick},
        Table(tts, (thesis.example_input_integral.(tts).-thesis.example_input_integral.(tts.-τ))/τ)
    ),
    [raw"\coordinate (top41) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot42) at (axis cs:5,0);"],
    # reconstruction for lif
    {height="4cm",width="6cm", grid = "minor", yminorgrids="false",xlabel=raw"time $t$", minor_xtick=spikes_lif, ymin=thesis.example_yrange[1], ymax=thesis.example_yrange[2], xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick},
        Table(tts, thesis.example_input.(tts))
    ),
    Plot(
        {no_marks, very_thick, color="black"},
        Table(tts, y_rec_lif.(tts))
    ),
    PlotInc(
        {no_marks, color="Set2-B", very_thick},
        Table(tts, y_filtered.(tts))
    ),
    [raw"\coordinate (top42) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot42) at (axis cs:5,0);"],
    # reconstruction for lnp
    {height="4cm",width="6cm", grid = "minor", yminorgrids="false",xlabel=raw"time $t$", minor_xtick=spikes_lnp, ymin=thesis.example_yrange[1], ymax=thesis.example_yrange[2], xmin=thesis.example_trange[1], xmax=thesis.example_trange[2]},
    PlotInc(
        {no_marks, very_thick},
        Table(tts, thesis.example_input.(tts))
    ),
    Plot(
        {no_marks, very_thick, color="black"},
        Table(tts, y_rec_lnp.(tts))
    ),
    PlotInc(
        {no_marks, color="Set2-B", very_thick},
        Table(tts, y_filtered.(tts))
    ),
    [raw"\coordinate (top43) at (axis cs:5,\pgfkeysvalueof{/pgfplots/ymax});", raw"\coordinate (bot43) at (axis cs:5,0);"],
    ),
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot11) -- node[draw, fill=white, rounded corners, align=center] {analog-to-digital\\conversion} (top31);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot31) -- (top41);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot12) -- (top22);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot22) -- (top32);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot32) -- (top42);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot13) -- (top23);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot23) -- (top33);",
    raw"\draw[->, shorten >= 4pt, shorten <= 4pt] (bot33) -- (top43);",
    "\\node[anchor=north east] at (label31) {($(sum(sum.(ys_bits_dig))) bits)};",
    "\\node[anchor=north east] at (label32) {($(length(spikes_lif)) spikes)};",
    "\\node[anchor=north east] at (label33) {($(length(spikes_lnp)) spikes)};",
)
#pgfsave("fig/encoding_schemes.pdf",gp)

max_stars_repo_name: jleugeri/phd

max_issues_repo_name: jleugeri/phd

max_forks_repo_name: jleugeri/phd