text:
# from https://github.com/chengchingwen/Transformers.jl/tree/master/src/fix

for typ ∈ atomictypes
    lt = llvmtypes[typ]
    ilt = llvmtypes[inttype(typ)]
    # Note: atomic_cas! succeeded (i.e. it stored "new") if and only if the result is "cmp"
    if typ <: Integer
        @eval atomic_cas!(x::Ptr{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                     %rs = cmpxchg $lt* %ptr, $lt %1, $lt %2 acq_rel acquire
                     %rv = extractvalue { $lt, i1 } %rs, 0
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     x, cmp, new)
    else
        @eval atomic_cas!(x::Ptr{$typ}, cmp::$typ, new::$typ) =
            llvmcall($"""
                     %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                     %icmp = bitcast $lt %1 to $ilt
                     %inew = bitcast $lt %2 to $ilt
                     %irs = cmpxchg $ilt* %iptr, $ilt %icmp, $ilt %inew acq_rel acquire
                     %irv = extractvalue { $ilt, i1 } %irs, 0
                     %rv = bitcast $ilt %irv to $lt
                     ret $lt %rv
                     """, $typ, Tuple{Ptr{$typ},$typ,$typ},
                     x, cmp, new)
    end

    arithmetic_ops = [:add, :sub]
    for rmwop in [arithmetic_ops..., :xchg, :and, :nand, :or, :xor, :max, :min]
        rmw = string(rmwop)
        fn = Symbol("atomic_", rmw, "!")
        if (rmw == "max" || rmw == "min") && typ <: Unsigned
            # LLVM distinguishes signedness in the operation, not the integer type.
            rmw = "u" * rmw
        end
        if rmwop in arithmetic_ops && !(typ <: ArithmeticTypes) continue end
        if typ <: Integer
            @eval $fn(x::Ptr{$typ}, v::$typ) =
                llvmcall($"""
                         %ptr = inttoptr i$WORD_SIZE %0 to $lt*
                         %rv = atomicrmw $rmw $lt* %ptr, $lt %1 acq_rel
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, x, v)
        else
            rmwop == :xchg || continue
            @eval $fn(x::Ptr{$typ}, v::$typ) =
                llvmcall($"""
                         %iptr = inttoptr i$WORD_SIZE %0 to $ilt*
                         %ival = bitcast $lt %1 to $ilt
                         %irv = atomicrmw $rmw $ilt* %iptr, $ilt %ival acq_rel
                         %rv = bitcast $ilt %irv to $lt
                         ret $lt %rv
                         """, $typ, Tuple{Ptr{$typ}, $typ}, x, v)
        end
    end
end

const opnames = Dict{Symbol, Symbol}(:+ => :add, :- => :sub, :* => :mul, :/ => :div)
for op in [:+, :-, :max, :min, :*, :/]
    opname = get(opnames, op, op)
    for typ ∈ atomictypes
        typ == Bool && continue
        !(op == :/ || op == :*) && typ <: Integer && continue
        @eval function $(Symbol("atomic_", opname, "!"))(var::Ptr{T}, val::T) where T<: $typ
            IT = inttype(T)
            old = unsafe_load(var)
            while true
                new = $op(old, val)
                cmp = old
                old = atomic_cas!(var, cmp, new)
                reinterpret(IT, old) == reinterpret(IT, cmp) && return new
                # Temporary solution before we have gc transition support in codegen.
                ccall(:jl_gc_safepoint, Cvoid, ())
            end
        end
    end
end

for op = [:add, :sub, :max, :min, :mul, :div]
    @eval function $(Symbol("scatter_", op, "!"))(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple},
                                                  s::Int=size(ys,1)) where T
        Threads.@threads for i = 1:s
            @inbounds for ind = CartesianIndices(xs)
                $(Symbol("atomic_", op, "!"))(
                    pointer(ys, Base._to_linear_index(ys, i, xs[ind]...)),
                    us[i, ind]
                )
            end
        end
        ys
    end
end

function scatter_mean!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple},
                       s::Int=size(ys,1)) where T
    Ns = zero(ys)
    ys_ = zero(ys)
    scatter_add!(Ns, one.(us), xs, s)
    scatter_add!(ys_, us, xs, s)
    ys .+= map((x,y) -> ifelse(iszero(y), x, x/y), ys_, Ns)
    return ys
end

@adjoint function scatter_add!(ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    ys_ = copy(ys)
    scatter_add!(ys_, us, xs)
    ys_, Δ -> (Δ, gather(zero(Δ)+Δ, xs), nothing)
end

@adjoint function scatter_sub!(ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    ys_ = copy(ys)
    scatter_sub!(ys_, us, xs)
    ys_, Δ -> (Δ, -gather(zero(Δ)+Δ, xs), nothing)
end

@adjoint function scatter_mul!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where T
    ys_ = copy(ys)
    scatter_mul!(ys_, us, xs)
    ys_, function (Δ)
        Δy = Δ .+ zero(ys)
        scatter_mul!(Δy, us, xs)
        rev_xs = gather_indices(xs)
        Δu = gather(ys, xs) .* gather(zero(Δ)+Δ, xs)
        @inbounds for ind = CartesianIndices(xs)
            inds = filter(x -> x != ind, rev_xs[xs[ind]])
            for i = 1:size(us, 1)
                Δu[i, ind] *= prod(j -> us[i, j], inds)
            end
        end
        (Δy, Δu, nothing)
    end
end

@adjoint function scatter_div!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where T
    ys_ = copy(ys)
    scatter_div!(ys_, us, xs)
    ys_, function (Δ)
        Δy = Δ .+ zero(ys)
        scatter_div!(Δy, us, xs)
        rev_xs = gather_indices(xs)
        Δu = - gather(ys, xs) .* gather(zero(Δ)+Δ, xs) ./ us.^2
        @inbounds for ind = CartesianIndices(xs)
            inds = filter(x -> x != ind, rev_xs[xs[ind]])
            for i = 1:size(us, 1)
                Δu[i, ind] /= prod(j -> us[i, j], inds)
            end
        end
        (Δy, Δu, nothing)
    end
end

function gather_indices(X::Array{T}) where T
    Y = DefaultDict{T,Vector{CartesianIndex}}(CartesianIndex[])
    @inbounds for (ind, val) = pairs(X)
        push!(Y[val], ind)
    end
    Y
end

@adjoint function scatter_max!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where T
    max = copy(ys)
    scatter_max!(max, us, xs)
    max, function (Δ)
       Δy = (ys .== max) .* Δ
       Δu = (us .== gather(max, xs)) .* gather(zero(Δ)+Δ, xs)
       (Δy, Δu, nothing)
    end
end

@adjoint function scatter_min!(ys::Array{T}, us::Array{T}, xs::Array{<:IntOrTuple}) where T
    min = copy(ys)
    scatter_min!(min, us, xs)
    min, function (Δ)
       Δy = (ys .== min) .* Δ
       Δu = (us .== gather(min, xs)) .* gather(zero(Δ)+Δ, xs)
       (Δy, Δu, nothing)
    end
end

@adjoint function scatter_mean!(ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    ys_ = copy(ys)
    scatter_mean!(ys_, us, xs)
    ys_, function (Δ)
        Δu = gather(zero(Δ)+Δ, xs)
        counts = zero.(xs)
        @inbounds for i = 1:size(ys, 2)
            counts += sum(xs.==i) * (xs.==i)
        end
        @inbounds for ind = CartesianIndices(counts)
            Δu[:, ind] ./= counts[ind]
        end
        (Δ, Δu, nothing)
    end
end

function scatter!(op::Symbol, ys::AbstractArray, us::AbstractArray, xs::AbstractArray)
    if op == :add
        return scatter_add!(ys, us, xs)
    elseif op == :sub
        return scatter_sub!(ys, us, xs)
    elseif op == :mul
        return scatter_mul!(ys, us, xs)
    elseif op == :div
        return scatter_div!(ys, us, xs)
    elseif op == :max
        return scatter_max!(ys, us, xs)
    elseif op == :min
        return scatter_min!(ys, us, xs)
    elseif op == :mean
        return scatter_mean!(ys, us, xs)
    end
end


text:
# using AuditoryFilters     # seems to break dependencies, find alternative

# include("/home/dahlbom/.julia/dev/DelNetExperiment/src/ExternalPeripheralModels.jl")
# include("/home/dahlbom/research/carfac/JCARFAC/src/JCARFAC.jl")
# EPM = ExternalPeripheralModels

################################################################################
# Interface to Zilany and Meddis Models
################################################################################

# function zilanymodel(sig::Array{Float64,1}, fs, db, f_lo, f_hi, numchannels;
#                      anf_num = (0, 10, 0), species="human", seed=1)
#     params = EPM.ZilanyParams(fs, db, anf_num, (f_lo, f_hi, numchannels),
#                               species, seed)
#     cfs_an, spikes = EPM.sigtospikedata(sig, params) |> EPM.spikedatatochannels
# end

################################################################################
# Tools for simple sub-cortical model
################################################################################

function inhomogenouspoisson(λs, fs; rt = 0.0, densityfactor = 5000.0)
	probs = (λs ./ fs) * densityfactor
	ts = collect(0.0:1/fs:(length(probs)-1)/fs)
	(maximum(probs) > 1.0) && println("WARNING: Maximum ($(maximum(probs))) greater than 1.0")
	(minimum(probs) < 0.0) && println("WARNING: Minimum prob. less than 0.0")
	spiketimes = Array{Float64,1}(undef,0)
	lastfired = -1.0
	for (k,t) ∈ enumerate(ts)
		if rand() < probs[k]
			if t > lastfired + rt
				push!(spiketimes,t)
				lastfired = t
			end
		end
	end
	spiketimes
end

# function auditoryspikes(signal::Vector{Float64}, numchannels, nervesperchannel, fs;
# 						rt=0.0, densityfactor=1000.0, lowf=100.0)
# 	periout = simpleperiphery(signal, numchannels, fs; lowf=lowf)
# 	numsamples, numchannels = size(periout)
# 	n = 1
# 	out = Array{Spike,1}(undef,0)
# 	for chan ∈ numchannels:-1:1
# 		#k = numchannels-chan+1
# 		k = chan
# 		for nerve ∈ 1:nervesperchannel
# 			spiketimes = inhomogenouspoisson(periout[:,k], fs;
# 										 rt=rt, densityfactor=densityfactor)
# 			for time ∈ spiketimes
# 				push!(out, Spike(n, time))
# 			end
# 			n += 1
# 		end
# 	end
# 	out
# end
#
# function auditoryspikes(signal::Vector{Float64}, numchannels, nervesperchannel,
# 						outputchannels, fs;
# 						rt=0.0, densityfactor=100.0, lowf=100.0)
# 	spikes = auditoryspikes(signal, numchannels, nervesperchannel, fs;
# 							rt=rt, densityfactor=densityfactor, lowf=lowf)
# 	@assert length(outputchannels) >= numchannels * nervesperchannel
# 	[Spike(outputchannels[s.n], s.t) for s ∈ spikes]
# end
#
# function simpleperiphery(signal::Vector{Float64}, numchannels, fs; lowf=100.0)
# 	fb = make_erb_filterbank(fs, numchannels, lowf)
# 	fbout = filt(fb, signal)
# 	map(x -> x < 0.0 ? 0.0 : x, fbout)
# end

function plotperiphery(periout::Array{T,2}, fs) where T <: Number
	p = plot(; legend=:none)
	numsamples, numchannels = size(periout)
	ts = 0:1/fs:(numsamples-1)/fs
	for chan ∈ numchannels:-1:1
		plot!(p, ts, periout[:,chan] .+ 0.25*(numchannels - chan + 1))
	end
	p
end

function hc(f0, nh, dur, fs;
			norm=true,
			weights=[1.0 for k ∈ 1:nh],
			ϕs=[0.0 for k ∈ 1:nh])
	ts = 0.0:1/fs:dur-1.0/fs
	sig = zeros(length(ts))
	for (k,w) ∈ enumerate(weights)
		sig .+= w .* sin.(2π*f0*k .* ts .+ ϕs[k])
	end
	if norm
		sig ./= maximum(sig)
	end
	sig
end

ext: jl

max_stars_repo_name: analogouscircuit/DelayNetTools

max_issues_repo_name: analogouscircuit/DelayNetTools

max_forks_repo_name: analogouscircuit/DelayNetTools


text:
using LinearAlgebra, QuadGK, StatsBase, Distributions, Statistics
import Base: rand, +, /, -
import Distributions: cdf, ccdf

mutable struct MAPHDist
    α::Adjoint{Float64, Vector{Float64}}
    T::Matrix{Float64}
    T0::Matrix{Float64}
end

model_size(maph::MAPHDist) = (p = size(maph.T,1), q = size(maph.T0,2)) #transient is p and abosrbing is q

"""

Returns parameters of a two phase hyper-exponential fitting a mean and an SCV.

"""

function hyper_exp_fit(mean::Float64, scv::Float64)
    scv < 1.0 && error("SCV must be greater than 1")
    μ1 = 1/(scv+1)
    p = (scv-1)/(scv+1+2/μ1^2-4/μ1)
    μ2 = 2*(1-μ1)/(2-μ1*(scv+1))
    α = zeros(2)'
    α[1] = 1-p
    α[2] = p

    T = zeros(2,2)
    T[1,1] = -1/μ1
    T[2,2] = -(1-p)/(1-2*p)

    return (α, (1/mean)*T)
end

"""

Returns parameters of a hypo-exponential (generalized erlang) dist which is a sum of n exponentials with the last one different

"""

function hypo_exp_fit(mean::Float64,scv::Float64)

    scv ≥ 1.0 && error("SCV must be less than 1")

    n = Int(ceil(1/scv))

    ν1 = n/(1+sqrt((n-1)*(n*scv-1)))
    ν2 = ν1*(n-1)/(ν1-1)

    α = zeros(n)'
    α[1] = 1
    T = zeros(n,n)
    T[1,1] = -ν1
    T[1,2] = ν1

    for i = 2:(n-1)
        T[i,i] = -ν2
        T[i,i+1] = ν2
    end

    T[n,n] = -ν2

    return (α, (1/mean)*T)

end

function scv(data::Vector{Float64})
    return (std(data)/mean(data))^2
end

"""
Merge PH distributions into a MAPH distribution
"""
function Merge_Dist(probs::Vector{Float64},dist::Vector{Any})
    merged_α = (probs./sum(probs)).*[reshape(dist[i][1],length(dist[i][1])) for i = 1:length(dist)]
    α_temp = reduce(vcat,merged_α)'
    merged_T = (probs./sum(probs)).*[dist[i][2] for i = 1:length(dist)]

    T_temp = cat(merged_T...,dims = (1,2))

    merged_T0 = -(probs./sum(probs)).*[sum(dist[i][2],dims=2) for i = 1:length(dist)]

    T0_temp = cat(merged_T0...,dims=(1,2))

    return α_temp, T_temp, T0_temp
end

"""

Create an MAPHDist of dimensions pxq where q is the length of `probs`, `means`, and `scvs` and p is specified.

This tries to do a best "moment fit" for the probability of absorbitions, means, and scvs

"""
function MAPHDist(p::Int, probs::Vector{Float64}, means::Vector{Float64}, scvs::Vector{Float64})
    q = length(probs)
    length(means) != q && error("Dimension mismatch")
    length(scvs) != q && error("Dimension mismatch")

    πhat_order = sortperm(probs)
    sorted_πhat = probs[πhat_order]
    sorted_scvs = scvs[πhat_order]
    sorted_means = means[πhat_order]

    num_phases = [sorted_scvs[i] ≥ 1 ? 2 : ceil(1/sorted_scvs[i]) for i in 1:q]

    required_phases = sum(num_phases)
    @show required_phases
    # selected_cases= [required_phases - sum(num_phases[1:k]) - p for k=1:q]

    K=0

    K = findfirst((x)->x ≤ p, [sum(num_phases[k:end]) for k=1:q])
    @show K
    @show num_phases, πhat_order,[sum(num_phases[k:q]) for k=1:q]

    K ∉ (1:q) && error("Not enough phases to start")

    # α = zeros(p)'
    # T = ones(p,p).*eps()
    # T0 = ones(p,q).*eps()

    dist = []
    for k = 1:q
        if sorted_scvs[k]≥1
            push!(dist,hyper_exp_fit(sorted_means[k],sorted_scvs[k]))
        end
        if sorted_scvs[k]<1
            push!(dist,hypo_exp_fit(sorted_means[k],sorted_scvs[k]))
        end
    end

    # if K>0
    #     for i = 1:K
    #         dist[i] = (ones(1,1),ones(1,1)).*eps()
    #     end
    # end
    α = zeros(p)'
    T = ones(p,p).*eps()
    T0 = ones(p,q).*eps()
    @show typeof(dist)

    reverse_order = sortperm(πhat_order)

    reversed_dist = dist[reverse_order]

    # reversed_α = (probs./sum(probs)).*[reshape(reversed_dist[i][1],length(reversed_dist[i][1])) for i = 1:q]
    # α_temp = reduce(vcat,reversed_α)'

    # reversed_T = (probs./sum(probs)).*[reversed_dist[i][2] for i = 1:q]

    # T_temp = cat(reversed_T...,dims = (1,2))

    # reversed_T0 = -(probs./sum(probs)).*[sum(reversed_dist[i][2],dims=2) for i = 1:q]

    # T0_temp = cat(reversed_T0...,dims=(1,2))

    if p ≥ required_phases
        α_temp, T_temp, T0_temp = Merge_Dist(probs,reversed_dist)
        m,n = model_size(MAPHDist(α_temp, T_temp, T0_temp))
        α[1:m] = α_temp
        T[1:m,1:m] = T_temp
        T0[1:m,1:n] = T0_temp
    else
        idx_drop = πhat_order[1:K-1]
        @show idx_drop
        probs = probs[eachindex(probs) .∉ Ref(idx_drop)]
        reversed_dist = reversed_dist[eachindex(reversed_dist) .∉ Ref(idx_drop)]
        α_temp, T_temp, T0_temp = Merge_Dist(probs,reversed_dist)
        m,n = model_size(MAPHDist(α_temp, T_temp, T0_temp))
        α[1:m] =α_temp
        T[1:m,1:m] = T_temp
        T0[1:m,1:n]=T0_temp

    end

    maph = MAPHDist(α,T,T0)

    return maph
end

# function MAPHInit(p::Int=1,q::Int=1)
#     α = rand(p)
#     α = (α/sum(α))'
#     T = rand(1:0.01:20,p,p)
#     T0 = rand(1:0.01:20,p,q)
#     for i = 1:p
#         temp = T[i,i]
#         T[i,i]= -sum(T[i,:])-sum(T0[i,:])+temp
#     end
#     return MAPHDist(α,T,T0)
# end

# MAPHDist(;p::Int=1,q::Int=1) = MAPHInit(p::Int=1,q::Int=1)
# MAPHDist(;p::Int=1,q::Int=1) = MAPHDist(
#                         Vector{Float64}(undef,p)',
#                         Matrix{Float64}(undef,p,p),
#                         Matrix{Float64}(undef,p,q))

# ccdf(maph::MAPHDist, y::Float64, a::Int) = exp(maph.α*maph.T*y)
# cdf = 1-ccdf
# V_j///

function q_matrix(d::MAPHDist)::Matrix{Float64}
    p, q = model_size(d)
    return [zeros(q,p) zeros(q,q) ; d.T0 d.T]
end

function p_matrix(d::MAPHDist)::Matrix{Float64}
    p, q = model_size(d)
    PT = I-inv(Diagonal(d.T))*d.T
    PT0 = -inv(Diagonal(d.T))*d.T0
    return [I zeros(q,p); PT0 PT]
end

function rand(d::MAPHDist; full_trace = false)
    p, q = model_size(d)
    transient_states = (q+1):(q+p)
    all_states = 1:(q+p)
    Pjump = p_matrix(d)
    Λ = vcat(zeros(q), -diag(d.T))

    if full_trace
        states = Int[]
        sojourn_times = Float64[]
    end

    state = sample(transient_states,weights(d.α)) #initial state
    t = 0.0
    while state ∈ transient_states
        sojourn_time = rand(Exponential(1/Λ[state]))
        if full_trace
            push!(states,state)
            push!(sojourn_times,sojourn_time)
        end
        t += sojourn_time
        state = sample(all_states,weights(Pjump[state,:]))
    end

    if full_trace
        push!(states,state)
        return (sojourn_times, states)
    else
        return (y = t, a = state)
    end
end

"""
Create a named tuple observation.
"""
observation_from_full_traj(times::Vector{Float64}, states::Vector{Int64}) = (y = sum(times),a = last(states))

mutable struct MAPHSufficientStats
    B::Vector{Float64} #initial starts
    Z::Vector{Float64} #time spent
    N::Matrix{Float64} #transitions
    MAPHSufficientStats(B::Vector{Float64}, Z::Vector{Float64}, N::Matrix{Float64}) = new(B,Z,N)
    function MAPHSufficientStats(maph::MAPHDist)
        p, q = model_size(maph)
        new(zeros(p), zeros(p), zeros(p,p+q))
    end
end

+(ss1::MAPHSufficientStats, ss2::MAPHSufficientStats) = MAPHSufficientStats(ss1.B+ss2.B, ss1.Z+ss2.Z, ss1.N+ ss2.N)
/(ss::MAPHSufficientStats,n::Real) = MAPHSufficientStats(ss.B/n,ss.Z/n,ss.N/n)
/(ss1::MAPHSufficientStats,ss2::MAPHSufficientStats) = MAPHSufficientStats(ss1.B ./ ss2.B, ss1.Z ./  ss2.Z, ss1.N ./  ss2.N)
-(ss1::MAPHSufficientStats, ss2::MAPHSufficientStats) = MAPHSufficientStats(ss1.B-ss2.B, ss1.Z-ss2.Z, ss1.N - ss2.N)

function sufficient_stat_from_trajectory(d::MAPHDist, sojourn_times::Array{Float64}, states::Array{Int})::MAPHSufficientStats
    p,q = model_size(d)
    transient_states = (q+1):(q+p)

    ss =  MAPHSufficientStats(d)

    for s = 1:p
        if states[1] == transient_states[s]
            ss.B[s] +=1
        end
    end

    for i = 1:(length(states)-1)
        ss.Z[states[i]-q] += sojourn_times[i]
        ss.N[states[i]-q, states[i+1]] += 1
    end

    return ss
end

# model_size(maph::MAPHDist) = (p = size(maph.T,1), q = size(maph.T0,2)) #transient is p and abosrbing is q

SingleObs = NamedTuple{(:y, :a), Tuple{Float64, Int64}}
MAPHObsData = Vector{SingleObs}

#temp
function very_crude_c_solver(y::Float64,i::Int,j::Int,maph::MAPHDist)
    quadgk(u -> (maph.α*exp(maph.T*u))[i]*exp(maph.T*(y-u))*maph.T0[:,j] , 0, y, rtol=1e-8) |> first
end

function sufficient_stats(observation::SingleObs, maph::MAPHDist; c_solver = very_crude_c_solver)::MAPHSufficientStats
    stats = MAPHSufficientStats(maph)

    p, q = model_size(maph)

    a(y::Float64) = maph.α*exp(maph.T*y)
    b(y::Float64,j::Int) = exp(maph.T*y)*maph.T0[:,j]
    c(y::Float64,i::Int,j::Int) = very_crude_c_solver(y,i,j,maph)

    D = Diagonal(maph.T)
    PT = I-inv(Diagonal(maph.T))*maph.T
    PT0 = -inv(Diagonal(maph.T))*maph.T0
    A = inv(I-PT)*PT0
    PA = maph.α*A

    # display(PT)
    # display(PT0)
    # display(PA)
    # display(sum(PA))

    EB(y::Float64, i::Int, j::Int) = maph.α[i] * b(y,j)[i]*PA[j]/ (maph.α*b(y,j))
    EZ(y::Float64, i::Int, j::Int) = c(y,i,j)[i]*PA[j]/(maph.α*b(y,j))
    # ENT(y::Float64,i::Int,j::Int) = i+q != j ?  maph.T[i,:].*c(y,i,j)*PA[j]/(maph.α*b(y,j)) : zeros(p)
    # ENT2(y::Float64,i::Int,j::Int) =  maph.T[i,:].*c(y,i,j)*PA[j]/(maph.α*b(y,j))

    ENT(y::Float64,i::Int,k::Int,j::Int) = i !=k ? maph.T[i,:].*c(y,i,j)*PA[j]/(maph.α*b(y,j)) : zeros(p)

    # ENA(y::Float64,i::Int,j::Int) = a(y)[i]*maph.T0[:,j][i]/(maph.α*b(y,j))
    ENA(y::Float64,i::Int,j::Int) = PA[j]*a(y)[i]*maph.T0[i,j]/(maph.α*b(y,j))

    stats.B = [sum([EB(observation.y, i, j) for j = 1:q]) for i =1:p]
    stats.Z = [sum([EZ(observation.y,i,j) for j =1:q]) for i = 1:p]

    for i = 1:p
        for k = (q+1):(q+p)
            V = sum([ENT(observation.y,i,k-q,j) for j in 1:q])
            stats.N[i,k] = V[k-q]
        end

        for j = 1:q
            stats.N[i,j] = ENA(observation.y,i,j)
        end

    end

    # for i= 1:p
    #     V = sum([ENT(observation.y,i,j) for j in 1:q])
    #     for k = q+1:q+p
    #         stats.N[i,k] = V[k-q]
    #     end

    #     for j = 1:q
    #         stats.N[i,j] = ENA(observation.y,i,j)
    #     end
    # end

    return stats
end

sufficient_stats(data::MAPHObsData, maph::MAPHDist; c_solver = very_crude_c_solver)::MAPHSufficientStats = mean([sufficient_stats(d, maph) for d in data])

# mutable struct ProbabilityMatrix
#     P::Matrix{Float64} # transient to transient
#     P0::Matrix{Float64}#transient to absorbing
#     function ProbabilityMatrix(maph::MAPHDist)
#         p, q = model_size(maph)
#         new(zeros(p,p), zeros(p,q))
#     end
# end

# function update_probability_matrix(maph::MAPHDist, probmatrix::ProbabilityMatrix)
#     probmatrix.P = I-inv(Diagonal(maph.T))*maph.T
#     probmatrix.P0 = -inv(Diagonal(maph.T))*maph.T0
# end

function stats_to_dist(maph::MAPHDist,ss::MAPHSufficientStats)::MAPHDist
    p,q = model_size(maph)
    α = ss.B'
    T = ss.N[:,(q+1):(q+p)]./ss.Z
    T0 = ss.N[:,1:q]./ss.Z
    for i = 1:p
        temp = T[i,i]
        T[i,i]= -sum(T[i,:])-sum(T0[i,:])+temp
    end
    return MAPHDist(α,T,T0)
end
"""
Fits ... QQQQ
"""
function fit!(maph::MAPHDist,data::MAPHObsData)::MAPHDist
    #EM Loop

    p,q = model_size(maph)

    for k in 1:100
        ss = sufficient_stats(data[1],maph)
        for i in 2:10^2
            ss = ss+sufficient_stats(data[i],maph)
        end
        ss = ss/10^2
                maph = stats_to_dist(maph,ss)

    end

    return maph
end

"""
Returns an array of data filtered according to absorbing state being the index of the array.
"""
function absorb_filter_data(data, maph::MAPHDist)
    p, q = model_size(maph)
    filter_data = []

    for i = 1:q
        temp_data = filter(data) do obs first(obs).a == i end
        push!(filter_data,temp_data)
    end

    return filter_data
end

function time_filter_data(data, n::Int64)

    max_time = maximum(((d)->first(d).y).(data) )#[data[i].y for i = 1:length(data)])
    time_vec = Array(LinRange(0, max_time, n)) #QQQQ - cleanup later

    filter_data = []

    for i = 1:(n-1)
        temp_data = filter(data) do obs first(obs).y ≥ time_vec[i] && first(obs).y < time_vec[i+1] end
        push!(filter_data,temp_data)
    end

    return(filter_data)
end

function test_example()
    Λ₄, λ45, λ54, Λ₅ = 5, 2, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]

    maph = MAPHDist([0.5,0.5]',T_example, T0_example)
    stats = MAPHSufficientStats(maph)
    # @show maph
    # @show model_size(maph)

    data = [(y=2.3,a=2),(y=5.32,a=1),(y=15.32,a=2)]
    # update_sufficient_stats(maph, data,stats)

    # @show stats

    # QQQQ update_sufficient_stats(maph,data,stats)

    test_stats = MAPHSufficientStats(maph)

    for _ in 1:10^5
        times, states = rand(maph, full_trace = true)
        ss = sufficient_stat_from_trajectory(maph, times, states)
        test_stats.N += ss.N
        test_stats.Z += ss.Z
        test_stats.B += ss.B
    end

    test_stats.B = test_stats.B/sum(test_stats.B)
    test_stats.Z = test_stats.Z/sum(test_stats.Z)
    test_stats.N = test_stats.N/sum(test_stats.N)

    # @show(stats.N,test_stats.N)
    @show test_stats
end

# test_example()

function test_example2()
    Λ₄, λ45, λ54, Λ₅ = 5, 2, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]

    maph = MAPHDist([0.5,0.5]', T_example, T0_example)

    data = []
    full_trace =[]

    println("starting generating data")
    for i in 1:2*10^4
        times, states = rand(maph, full_trace = true)
        push!(full_trace, (times,states))
        push!(data, (observation_from_full_traj(times,states),i))
        i % 10^5 == 0 && print(".")
    end
    println("\nfinished generating data")

    absorb = absorb_filter_data(data, maph)
    time_bin = time_filter_data(absorb[1], 1000)

    #loop over all bins
    println("\n start initialization...")
    for i in 1:length(time_bin)
        ss_i = MAPHSufficientStats[]
        for trace in full_trace[last.(time_bin[i])]
            ss = sufficient_stat_from_trajectory(maph, trace[1], trace[2])
            push!(ss_i, ss)
        end

        if !isempty(ss_i)
            mean_observed_ss = mean(ss_i)
            obs = first(data[last(last.(time_bin[i]))])
            computed_ss = sufficient_stats(obs, maph)
            errs_N = (mean_observed_ss.N - computed_ss.N) ./ computed_ss.N #./ computed_ss
            # @show obs
            # sufficient_stats()
        end
    end

    ab1 = filter(x->x.a==1,first.(data))
    ab2 = filter(x->x.a==2,first.(data))
    ab3 = filter(x->x.a==3,first.(data))

    probs = [length(ab1),length(ab2),length(ab3)]./length(data)
    means = [mean(first.(ab1)),mean(first.(ab2)),mean(first.(ab3))]
    scvs = [scv(first.(ab1)),scv(first.(ab2)),scv(first.(ab3))]

    dist = MAPHDist(10,probs,means,scvs)
    println("finish initialization")

    println("starting simulations")
    fit!(dist,first.(data))
    println("\nfinished simulations")

    return dist

    # @show first(data)
    # # @show m[1]
    # @show maximum(data.y)
    # @show maximum[data[i].y for i = 1:length(data)]
end

# test_example2();

# filter(ret) do obs obs.y ≥ 0.0 && obs.y < 0.005 && obs.a == 1 end

# a = (filter(ret) do obs obs.a ==1 end);

#QQQQ - An example with a deterministic path....
function test_example3()

    ϵ = 0.00000001

    T_example = [-(1.0+4ϵ) 1.0-4ϵ ϵ
                ϵ -(1.0 + 4ϵ) 1.0-4ϵ
                ϵ ϵ -(1 + 4ϵ) ]

    T0_example = [ϵ ϵ ϵ;
                  ϵ ϵ ϵ;
                  1-4ϵ ϵ ϵ]

    T_example2 = [-(1.0+3ϵ) 1.0 ϵ
                 ϵ -(1.0+3ϵ) 1.0
                 ϵ ϵ -(1+3ϵ)]
    T0_example2 = [ϵ ϵ;
                   ϵ ϵ;
                   1-3ϵ ϵ]

    display(sum(T0_example2[:,1]))
    maph = MAPHDist([0.5,0.5, 0.0]',T_example, T0_example)

    obs = (y=100.0, a=2)
    sufficient_stats(obs, maph)
end

# ss = test_example2()

function test_example4()
    Λ₄, λ45, λ54, Λ₅ = 5, 2, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]

    maph = MAPHDist([0.5,0.5]', T_example, T0_example)

    println("starting simulations")
    data = [rand(maph) for _ in 1:10^2]

    p,q = model_size(maph)
    est_maph = MAPHInit(p,q)
    # est_maph = MAPHDist(;model_size(maph)...)
    @assert model_size(est_maph) == model_size(maph)

    #fit! gets an MAPHDist object for two possible reasons:
        #Reason #1 (always) - to know p and q.
        #Reason #2 (sometimes) - to have a starting guess. QQQQ - later give it a flag to say
    fit!(est_maph,data)

    println("\nfinished simulations")
end

#ss=test_example4()

# a, T = hypo_exp_fit(0.23)

# @show T

function test_init()
    out = MAPHDist(4,[0.3,0.5,0.1,0.08,0.02],[2.0,3.1,2.3,4.5,0.2],[1.1,2.3,2.0,0.33,1.5])
    @show model_size(out)
    data = [rand(out) for _ in 1:10^5]
    pi_est = [count((x)->x.a == i, data)/length(data) for i in 1:5]
    μ_est = [mean(first.(filter((x)->x.a == i, data))) for i in 1:5]
    @show pi_est
    @show μ_est
    return out
end

maph = test_init()

@show model_size(maph)


text:
using Documenter, Soss

DocMeta.setdocmeta!(Soss, :DocTestSetup,
    quote
        using Soss
        import Random
        Random.seed!(3)
    end; recursive=true)

makedocs(;
    modules=[Soss],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/cscherrer/Soss.jl/blob/{commit}{path}#L{line}",
    sitename="Soss.jl",
    authors="Chad Scherrer",
    assets=String[],
)

deploydocs(;
    repo="github.com/cscherrer/Soss.jl",
    push_preview=true
)


text:
include("../src/QPnorm.jl")
include("../examples/subproblems.jl")
using Main.QPnorm
using Random, Test

function optimality_metrics(P, q, A, b, r_min, r_max, x, λ)
    m, n = size(A)

    f = dot(x, P*x)/2 + dot(x, q)

    grad_residual = norm(P*x + q + A'*λ[1:end-1] + λ[end]*x, Inf)
    infeasibility = max(maximum(A*x - b), norm(x)^2 - r_max^2, r_min^2 - norm(x)^2, 0)

    # Find active r, if any
    if abs(norm(x) - r_min) < abs(norm(x) - r_max)
        r = r_min
    else
        r = r_max
    end

    complementarity = maximum(minimum([abs.(λ) abs.([A*x - b; norm(x)^2 - r^2])], dims=2))
    dual_infeasibility = max(-minimum(λ), 0.0)

    active_set = λ .>= 1e-8
    A_active = [A; x'][active_set, :]
    if length(A_active) > 0
        V = nullspace([A; x'][active_set, :])
    else
        V = diagm(0 => ones(n))
    end
    if length(V) > 0
        min_eig = minimum(eigvals(Symmetric(V'*(P + I*λ[end])*V)))
    else
        min_eig = 0.0
    end

    return f, grad_residual, infeasibility, dual_infeasibility, complementarity, min_eig
end

rng = MersenneTwister(123)
tol = 1e-7
@testset "Optimality Conditions of a random problem" begin
    n = 100
    m = 200
    P = randn(rng, n, n); P = (P + P')/2;
    q = randn(rng, n)
    A = randn(rng, m, n)
    b = randn(rng, m)
    r_min = 10.0; r_max = 20.0;
    x_init = find_feasible_point(A, b, r_min, r_max)
    x, λ = Main.QPnorm.solve(P, q, A, b, x_init, r_min=r_min, r_max=r_max, printing_interval=100)
    f, grad, inf, dinf, compl, min_eig = optimality_metrics(P, q, A, b, r_min, r_max, x, λ)
    @test grad < tol
    @test inf < tol
    @test dinf < tol
    @test compl < tol
    @test -min_eig < tol
end


text:
using Documenter
using ReinforcementLearning

"filter concrete subtypes"
is_concrete_type_of(t) = x -> begin
    println(x)
    x isa Type && x <: t && !isabstracttype(x)
end

makedocs(
    modules = [
        ReinforcementLearning,
        ReinforcementLearningBase,
        ReinforcementLearningCore,
        ReinforcementLearningEnvironments,
        ReinforcementLearningZoo,
    ],
    format = Documenter.HTML(
        prettyurls = true,
        analytics = "UA-149861753-1",
        canonical = "https://juliareinforcementlearning.github.io/ReinforcementLearning.jl/latest/",
        assets = ["assets/favicon.ico", "assets/custom.css"],
    ),
    sitename = "ReinforcementLearning.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Overview" => "overview.md",
        "A Quick Example" => "a_quick_example.md",
        "Manual" => [
            "RLBase" => "rl_base.md",
            "RLCore" => "rl_core.md",
            "RLEnvs" => "rl_envs.md",
            "RLZoo" => "rl_zoo.md",
        ],
        "Tips for Developers" => "tips_for_developers.md",
        "Experiments" => ["Play Atari Games with DQN" => "experiments/atari_dqn.md"],
    ],
)

deploydocs(
    repo = "github.com/JuliaReinforcementLearning/ReinforcementLearning.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
using BulkSMS: BulkSMSClient, send
using BulkSMS: multiple, shorten

include("config.jl")

client = BulkSMSClient(CONFIG["MSISDN"], CONFIG["USERNAME"], CONFIG["PASSWORD"])
message = "My first message"
# response = send(client, message, msisdn="+336........")

@info "Send '$message'"
response = send(client, message)  # use client MSISDN by default

# Action when long message can be set using:
# response = send(client, message, action_when_long=shorten)  # shorten long message
# response = send(client, message, action_when_long=multiple)  # multiple messages

@info "response: $response"


text:
function translation(x::T,y::T,z::T)::Array{T,2} where {T<:Real}
    return [1 0 0 x; 0 1 0 y; 0 0 1 z; 0 0 0 1]
end

function translation(x::Real,y::Real,z::Real)
    return translation(promote(x,y,z)...)
end

function scaling(x::T,y::T,z::T)::Array{T,2} where {T<:Real}
    return [x 0 0 0; 0 y 0 0; 0 0 z 0; 0 0 0 1]
end

function scaling(x::Real,y::Real,z::Real)
    return scaling(promote(x,y,z)...)
end

function rotation_x(angle::T)::Array{T,2} where {T<:Real}
    return [1 0 0 0; 0 cos(angle) -sin(angle) 0; 0 sin(angle) cos(angle) 0; 0 0 0 1]
end

function rotation_y(angle::T)::Array{T,2} where {T<:Real}
    return [cos(angle) 0 sin(angle) 0; 0 1 0 0; -sin(angle) 0 cos(angle) 0; 0 0 0 1]
end

function rotation_z(angle::T)::Array{T,2} where {T<:Real}
    return [cos(angle) -sin(angle) 0 0; sin(angle) cos(angle) 0 0; 0 0 1 0; 0 0 0 1]
end

function shearing(x_y::T,x_z::T,y_x::T,y_z::T,z_x::T,z_y::T)::Array{T,2} where {T<:Real}
    return [1 x_y x_z 0; y_x 1 y_z 0; z_x z_y 1 0; 0 0 0 1]
end


text:
using DataStructures: Queue, enqueue!, dequeue!, isempty

include("../LoadGraph/LoadGraph.jl")
using .LoadGraph: load_adjacency_list, AdjacencyList

function main(dir_path="../../cleaned_data", dtype::Type=UInt32)
    for file in filter(x->endswith(x, ".txt"), readdir(dir_path))
        if file != "friendster.txt"
        println("=========== ", file, " ===========")
        graph = load_adjacency_list("$dir_path/$file", dtype)

        @time all_cc = all_connected_components(graph)
        println("Number of connected components: ", length(all_cc))
        fraction = maximum(length(cc) for cc in all_cc) / length(graph)
        println("Fraction of nodes in the largest CC: ",
                round(100*fraction, digits=2), "%")

        @time diam_lb = diameter_lb(graph)
        println("Diameter lower bound : ", diam_lb)
        println()
        end
    end
end

function BFS(graph::AdjacencyList{T}, node::T) where T
    fifo = Queue{T}()
    enqueue!(fifo, node)
    marks = Dict{T, Bool}([(node, true)])
    distances = Dict{T, Int}([(node, 0)])

    connected_component = Vector{T}()
    while ! isempty(fifo)
        u = dequeue!(fifo)
        push!(connected_component, u)
        for v in graph[u]
            if ! haskey(marks, v)
                enqueue!(fifo, v)
                marks[v] = true
                distances[v] = distances[u] + 1
            end
        end
    end
    return [connected_component, distances]
end

function all_connected_components(graph::AdjacencyList{T}) where T
    connected_components = Vector{T}[]
    marks = Dict{T, Bool}()
    for key in keys(graph)
        if length(marks) == length(graph)
            break
        elseif ! haskey(marks, key) # We check that the node is not marked yet
            # We mark this new node
            marks[key] = true
            cc, _ = BFS(graph, key)

            # We mark the nodes that are in the connected component returned by BFS
            for node in cc
                marks[node] = true
            end
            push!(connected_components, cc)
        end
    end

    return connected_components
end

"""Computes a lower bound of the diameter of a graph.

The distance is the number of edges between two nodes.
"""
function diameter_lb(graph::AdjacencyList)
    node0 = highest_degree_node(graph)
    node1, _ = furthest_node(graph, node0)
    _, dist = furthest_node(graph, node1)
    return maximum(values(dist))
end

function highest_degree_node(graph::AdjacencyList{T}) where T
    max_deg = 0
    arg_max_degree = collect(Iterators.take(keys(graph), 1))[1]
    for (node, neighbour) in graph
        if length(neighbour) > max_deg
            arg_max_degree = node
        end
    end
    return arg_max_degree
end

function furthest_node(graph::AdjacencyList{T}, node::T) where T
    #get the connected components of node and distances from node
    cc_node, dist = BFS(graph, node)
    furthest = cc_node[end]#get the node that is the furthest from node
    return furthest, dist
end

main()

max_stars_repo_name: quentin-burthier/GrandsGraphes

max_issues_repo_name: quentin-burthier/GrandsGraphes

max_forks_repo_name: quentin-burthier/GrandsGraphes


text:
struct UnitCell
    unit::Molecule
    nextunit::Molecule
    bonds::Set{Bond}
end

function hamiltonian(T::Type, f, uc::UnitCell)
    mol = uc.unit
    mol′ = uc.nextunit
    N = countorbitals(mol)
    H0 = fill(zero(T), N, N)
    t1 = fill(zero(T), N, N)

    i = 1
    for (atom1, atom1′) in zip(mol.atoms, mol′.atoms)
        Ni = countorbitals(atom1)

        j = 1
        for (atom2, atom2′) in zip(mol.atoms, mol′.atoms)
            Nj = countorbitals(atom2)

            if atom1 === atom2
                H0[i:i+Ni-1, j:j+Nj-1] += hamiltonian(T, f, atom1)
            elseif Bond(atom1, atom2) in mol.bonds
                H0[i:i+Ni-1, j:j+Nj-1] += hamiltonian(T, f, atom1, atom2)
            end

            if Bond(atom1, atom2′) in uc.bonds
                t1[i:i+Ni-1, j:j+Nj-1] += hamiltonian(T, f, atom1, atom2′)
            end

            j += Nj
        end

        i += Ni
    end
    return H0, t1
end

function overlap(T::Type, f, uc::UnitCell)
    mol = uc.unit
    mol′ = uc.nextunit
    N = countorbitals(mol)
    S0 = fill(zero(T), N, N)
    S1 = fill(zero(T), N, N)

    i = 1
    for (atom1, atom1′) in zip(mol.atoms, mol′.atoms)
        Ni = countorbitals(atom1)

        j = 1
        for (atom2, atom2′) in zip(mol.atoms, mol′.atoms)
            Nj = countorbitals(atom2)

            if atom1 === atom2
                S0[i:i+Ni-1, j:j+Nj-1] += overlap(T, f, atom1)
            elseif Bond(atom1, atom2) in mol.bonds
                S0[i:i+Ni-1, j:j+Nj-1] += overlap(T, f, atom1, atom2)
            end

            if Bond(atom1, atom2′) in uc.bonds
                S1[i:i+Ni-1, j:j+Nj-1] += overlap(T, f, atom1, atom2′)
            end

            j += Nj
        end

        i += Ni
    end
    return S0, S1
end


text:
# Put in a separate page so it can be used by SciMLDocs.jl

pages=[
    "Home" => "index.md",
    "ode_def.md",
]

max_stars_repo_name: JuliaDiffEq/ParameterizedFunctions

max_issues_repo_name: JuliaDiffEq/ParameterizedFunctions

max_forks_repo_name: JuliaDiffEq/ParameterizedFunctions


text:
module BitemporalPostgresDDL
include("../../src/DDL.jl")
using .DDL
export up, down
end


text:
const rcParams =
    @compat Dict{Symbol,Any}(:plotlib => "PyPlot",
                            :storeresults => false,
                            :chunk_defaults =>
                                Dict{Symbol,Any}(
                                :echo=> true,
                                :results=> "markup",
                                :fig=> true,
                                :include=> true,
                                :eval => true,
                                :fig_cap=> nothing,
                                #Size in inches
                                :fig_width => 6,
                                :fig_height => 4,
                                :fig_path=> "figures",
                                :dpi => 96,
                                :term=> false,
                                :name=> nothing,
                                :wrap=> true,
                                :engine=> "julia",
                                :option_string=> "",
                                #Defined in formats
                                :fig_ext => nothing,
                                :fig_pos=> nothing,
                                :fig_env=> nothing,
                                :out_width=> nothing,
                                :out_height=> nothing,
                                )
                            )

# Working towards Knitr compatible options, implemented options are
# added to defaultoptions dictionary above and work in progress stays here,
# options from https://github.com/yihui/knitr/blob/master/R/defaults.R
# If you need a particular options, consider implementing it and making a
# pull request.

#tidy = FALSE,
#tidy.opts = NULL,
#collapse = FALSE
#prompt = FALSE
#highlight = TRUE
#strip.white = TRUE
#size = 'normalsize'
#background = '#F7F7F7',
#cache = FALSE
#cache.path = 'cache/'
#cache.vars = NULL
#cache.lazy = TRUE,
#dependson = NULL
#autodep = FALSE,
#fig.keep = 'high'
#fig.show = 'asis'
#fig.align = 'default'
#dev = NULL
#dev.args = NULL
#fig.ext = NULL
#fig.scap = NULL
#fig.lp = 'fig:'
#fig.subcap = NULL,
#out.extra = NULL
#fig.retina = 1,
#external = TRUE
#sanitize = FALSE
#interval = 1
#aniopts = 'controls,loop',
#warning = TRUE
#error = TRUE
#message = TRUE,
#render = NULL,
#ref.label = NULL
#child = NULL
#split = FALSE
#purl = TRUE


text:
using IterativeSolvers

export bicgstabl, bicgstabl!

bicgstabl(A, b, Pl=1; kwargs...) =  bicgstabl!(IterativeSolvers.zerox(A,b), A, b, Pl; kwargs...)

function bicgstabl!(x, A, b, Pl=1; tol::Real=size(A,2)*eps(), maxiter::Int=size(A,2))
    K = KrylovSubspace(A, length(b), 1, Vector{IterativeSolvers.Adivtype(A,b)}[])
    IterativeSolvers.init!(K, x)
    bicgstabell(K; tol=tol, l=1)
end

function bicgstabell(K::KrylovSubspace; tol::Real=eps(), l::Int=1)
    k=-l
    r = zeros(3)
    IterativeSolvers.initrand!(K)
    #XXX Choose r̃[0]
    r[0]=b-IterativeSolvers.nextvec(K)
    u[-1]=0
    x[0]=K.v0
    ρ₀=1
    α=0
    ω=1

    while norm(r[k+l]) >= tol
        k+=l
        û[0]=u[k-1]; r̂[0]=r[k]; x̂[0]=x[k]
        ρ₀*=-ω
        for j=0:l-1 #BiCG part
            ρ₁=r̂[j]⋅r̃[0]; β=β[k+j]=α*ρ1/ρ₀; ρ₀=ρ₁
            for i=0:j
                û[i]=r̂[i]-β*û[i]
            end
            û[j+1]=K.A*û[j]
            γ=û[j+1]⋅r̃[0]; α=α[k+j]=ρ0/γ
            for i=0:j
                r̂[i]-=α*û[i+1]
            end
            r̂[j+1]=K.A*r̂[j]; x̂0+=α*û[0]
        end
        for j=1:l #Mod G-S; MR part
            for i=1:j-1
                τ[i,j]=r̂[j]⋅r̂[i]/σ[i]
                r̂[j]-=τ[i,j]*r̂[i]
            end
            σ[j]=r̂[j]⋅r̂[j];γ′[j]=r̂[0]⋅r̂[j]/σ[j]
        end
        ω=γ[l]=γ′[l]
        for j=l-1:-1:1 #γ=T\γ′
            γ[j]=γ′[j]-sum([τ[j,i]*γ[i] for i=j+1:l])
        end
        for j=1:l-1 #γ″=TSγ
            γ″[j]=γ[j+1]-sum([τ[j,i]*γ[i+1] for i=j+1:l-1])
        end
        x̂[0]+=γ[1]*r̂[0]; r̂[0]-=γ′[l]*r̂[l]; û[0]-=γ[l]*û[l] #update
        for j=1:l-1
            û[0]-=γ[j]*û[j]
            x̂[0]+=γ″[j]*r̂[j]
            r̂[0]-=γ′[j]*r̂[j]
        end
        u[k+l-1]=û[0]; r[k+l]=r̂[0]; x[k+l]=x̂[0]
    end
    #TODO What to return?
end


text:
module Trixi2Txt

using EllipsisNotation
using Glob: glob
using Printf: @printf
using HDF5: h5open, attrs, exists
using Tullio: @tullio
using LoopVectorization

include("../src/solvers/dg/interpolation.jl")

function trixi2txt(filename::AbstractString...;
                   variables=[], output_directory=".", nvisnodes=nothing, max_supported_level=11)
  # Convert filenames to a single list of strings
  if isempty(filename)
    error("no input file was provided")
  end
  filenames = String[]
  for pattern in filename
    append!(filenames, glob(pattern))
  end

  # Iterate over input files
  for (index, filename) in enumerate(filenames)
    # Check if data file exists
    if !isfile(filename)
      error("file '$filename' does not exist")
    end

    # Make sure it is a data file
    if !is_solution_restart_file(filename)
      error("file '$filename' is not a data file")
    end

    # Get mesh file name
    meshfile = extract_mesh_filename(filename)

    # Check if mesh file exists
    if !isfile(meshfile)
      error("mesh file '$meshfile' does not exist")
    end

    # Read mesh
    center_level_0, length_level_0, leaf_cells, coordinates, levels = read_meshfile(meshfile)

    # Read data
    labels, data, n_elements, n_nodes, element_variables, time = read_datafile(filename)

    # Check if dimensions match
    if length(leaf_cells) != n_elements
      error("number of elements in '$(filename)' do not match number of leaf cells in " *
            "'$(meshfile)' " *
            "(did you forget to clean your 'out/' directory between different runs?)")
    end

    # Determine resolution for data interpolation
    max_level = maximum(levels)
    if max_level > max_supported_level
      error("Maximum refinement level in data file $max_level is higher than " *
            "maximum supported level $max_supported_level")
    end
    max_available_nodes_per_finest_element = 2^(max_supported_level - max_level)
    if nvisnodes == nothing
      max_nvisnodes = 2 * n_nodes
    elseif nvisnodes == 0
      max_nvisnodes = n_nodes
    else
      max_nvisnodes = nvisnodes
    end
    nvisnodes_at_max_level = min(max_available_nodes_per_finest_element, max_nvisnodes)
    resolution = nvisnodes_at_max_level * 2^max_level
    nvisnodes_per_level = [2^(max_level - level)*nvisnodes_at_max_level for level in 0:max_level]

    # Interpolate data
    structured_data = unstructured2structured(data, levels, resolution, nvisnodes_per_level)

    # Interpolate cell-centered values to node-centered values
    node_centered_data = cell2node(structured_data)

    # Determine x coordinates
    xs = collect(range(-1, 1, length=resolution+1)) .* length_level_0/2 .+ center_level_0[1]

    # Check that all variables exist in data file
    if isempty(variables)
      append!(variables, labels)
    else
      for var in variables
        if !(var in labels)
          error("variable '$var' does not exist in the data file $filename")
        end
      end
    end

    # Create output directory if it does not exist
    mkpath(output_directory)

    # Determine output file name
    base, _ = splitext(splitdir(filename)[2])
    output_filename = joinpath(output_directory, "$(base).txt")

    # Write to file
    open(output_filename, "w") do io
      # Header
      print(io, "x             ")
      for label in variables
        @printf(io, "  %-14s", label)
      end
      println(io)

      # Data
      for idx in 1:length(xs)
        @printf(io, "%+10.8e", xs[idx])
        for variable_id in 1:length(variables)
          @printf(io, " %+10.8e", node_centered_data[idx, variable_id])
        end
        println(io)
      end
    end
  end
end

# Check if file is a data file
function is_solution_restart_file(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # If attribute "mesh_file" exists, this must be a data file
    return exists(attrs(file), "mesh_file")
  end
end

# Use data file to extract mesh filename from attributes
function extract_mesh_filename(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract filename relative to data file
    mesh_file = read(attrs(file)["mesh_file"])

    return joinpath(dirname(filename), mesh_file)
  end
end

# Read in mesh file and return relevant data
function read_meshfile(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Check dimension - only 1D supported
    if exists(attrs(file), "ndims")
      ndims_ = read(attrs(file)["ndims"])
    else
      ndims_ = read(attrs(file)["ndim"]) # FIXME once Trixi's 3D branch is merged & released
    end
    if ndims_ != 1
      error("currently only 1D files can be processed, but '$filename' is $(ndims_)D")
    end

    # Extract basic information
    n_cells = read(attrs(file)["n_cells"])
    n_leaf_cells = read(attrs(file)["n_leaf_cells"])
    center_level_0 = read(attrs(file)["center_level_0"])
    length_level_0 = read(attrs(file)["length_level_0"])

    # Extract coordinates, levels, child cells
    coordinates = Array{Float64}(undef, ndims_, n_cells)
    coordinates .= read(file["coordinates"])
    levels = Array{Int}(undef, n_cells)
    levels .= read(file["levels"])
    child_ids = Array{Int}(undef, 2^ndims_, n_cells)
    child_ids .= read(file["child_ids"])

    # Extract leaf cells (= cells to be plotted) and contract all other arrays accordingly
    leaf_cells = similar(levels)
    n_cells = 0
    for cell_id in 1:length(levels)
      if sum(child_ids[:, cell_id]) > 0
        continue
      end

      n_cells += 1
      leaf_cells[n_cells] = cell_id
    end
    leaf_cells = leaf_cells[1:n_cells]

    coordinates = coordinates[:, leaf_cells]
    levels = levels[leaf_cells]

    return center_level_0, length_level_0, leaf_cells, coordinates, levels
  end
end

# Read in data file and return all relevant information
function read_datafile(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract basic information
    if exists(attrs(file), "ndims")
      ndims_ = read(attrs(file)["ndims"])
    else
      ndims_ = read(attrs(file)["ndim"])
    end
    if exists(attrs(file), "polydeg")
      polydeg = read(attrs(file)["polydeg"])
    else
      polydeg = read(attrs(file)["N"])
    end
    n_elements = read(attrs(file)["n_elements"])
    n_variables = read(attrs(file)["n_vars"])
    time = read(attrs(file)["time"])

    # Extract labels for legend
    labels = Array{String}(undef, 1, n_variables)
    for v = 1:n_variables
      labels[1, v] = read(attrs(file["variables_$v"])["name"])
    end

    # Extract data arrays
    n_nodes = polydeg + 1

    if ndims_ == 1
      data = Array{Float64}(undef, n_nodes, n_elements, n_variables)
      for v = 1:n_variables
        vardata = read(file["variables_$v"])
        @views data[:, :, v][:] .= vardata
      end
    else
      error("Unsupported number of spatial dimensions: ", ndims_)
    end

    # Extract element variable arrays
    element_variables = Dict{String, Union{Vector{Float64}, Vector{Int}}}()
    index = 1
    while exists(file, "element_variables_$index")
      varname = read(attrs(file["element_variables_$index"])["name"])
      element_variables[varname] = read(file["element_variables_$index"])
      index +=1
    end

    return labels, data, n_elements, n_nodes, element_variables, time
  end
end

# Interpolate unstructured DG data to structured data (cell-centered)
function unstructured2structured(unstructured_data::AbstractArray{Float64},
                                 levels::AbstractArray{Int}, resolution::Int,
                                 nvisnodes_per_level::AbstractArray{Int})
  # Extract data shape information
  n_nodes_in, n_elements, n_variables = size(unstructured_data)

  # Get node coordinates for DG locations on reference element
  nodes_in, _ = gauss_lobatto_nodes_weights(n_nodes_in)

  # Calculate interpolation vandermonde matrices for each level
  max_level = length(nvisnodes_per_level) - 1
  vandermonde_per_level = []
  for l in 0:max_level
    n_nodes_out = nvisnodes_per_level[l + 1]
    dx = 2 / n_nodes_out
    nodes_out = collect(range(-1 + dx/2, 1 - dx/2, length=n_nodes_out))
    push!(vandermonde_per_level, polynomial_interpolation_matrix(nodes_in, nodes_out))
  end

  # Create output data structure
  structured = Array{Float64}(undef, resolution, n_variables)

  # For each variable, interpolate element data and store to global data structure
  for v in 1:n_variables
    first = 1

    # Reshape data array for use in interpolate_nodes function
    @views reshaped_data = reshape(unstructured_data[:, :, v], 1, n_nodes_in, n_elements)

    for element_id in 1:n_elements
      # Extract level for convenience
      level = levels[element_id]

      # Determine target indices
      n_nodes_out = nvisnodes_per_level[level + 1]
      last = first + (n_nodes_out - 1)

      # Interpolate data
      vandermonde = vandermonde_per_level[level + 1]
      @views structured[first:last, v] .= (
           reshape(multiply_dimensionwise(reshaped_data[:, :, element_id], vandermonde),
                   n_nodes_out))

      # Update first index for next iteration
      first += n_nodes_out
    end
  end

  return structured
end

# Convert cell-centered values to node-centered values by averaging over all
# four neighbors and making use of the periodicity of the solution
function cell2node(cell_centered_data::AbstractArray{Float64})
  # Create temporary data structure to make the averaging algorithm as simple
  # as possible (by using a ghost layer)
  tmp = similar(cell_centered_data, size(cell_centered_data) .+ (2, 0))

  # Fill center with original data
  tmp[2:end-1, :] .= cell_centered_data

  # Fill sides with opposite data (periodic domain)
  # x-direction
  tmp[1,   :] .= cell_centered_data[end, :]
  tmp[end, :] .= cell_centered_data[1,   :]

  # Create output data structure
  resolution_in, n_variables = size(cell_centered_data)
  resolution_out = resolution_in + 1
  node_centered_data = Array{Float64}(undef, resolution_out, n_variables)

  # Obtain node-centered value by averaging over neighboring cell-centered values
  for i in 1:resolution_out
    node_centered_data[i, :] = (tmp[i, :] + tmp[i+1, :]) / 2
  end

  return node_centered_data
end

end


text:
module DiscgolfRecord

using CSV, DataFrames, Distances, DataStructures, JSON, TOML, LightXML, TimeZones, Dates
using Interpolations: LinearInterpolation

export guess_and_check, parse_round_raw_csv, id

STATS_DIR = joinpath(homedir(),"dg_stats")

COURSES_DIR = joinpath(@__DIR__,"..","data","courses")
DISCS_CSV = joinpath(@__DIR__,"..","data","discs.csv")
SCORE_NAMES = Dict(
    -3 => "ALBATROSS",
    -2 => "EAGLE",
    -1 => "BIRDIE",
    0 => "PAR",
    1 => "BOGEY",
    2 => "DOUBLE BOGEY",
    3 => "TRIPLE BOGEY",
    4 => "QUADRUPLE BOGEY")

struct Location
    lat::Float64
    lon::Float64
end
lat(l::Location) = l.lat
lon(l::Location) = l.lon
dist(l1::Location, l2::Location) = Distances.haversine((lon(l1),lat(l1)),(lon(l2), lat(l2)))
Location(v) = Location(v[1], v[2])

struct POI
    id::String
    loc::Location
end
empty_poi() = POI("",Location(0.0,0.0))
loc(p::POI) = p.loc
id(p::POI) = p.id
lat(p::POI) = lat(loc(p))
lon(p::POI) = lon(loc(p))
dist(p1::POI, p2::POI) = dist(loc(p1), loc(p2))
function get_elem(elem_id, V::Vector{POI})
    for v in V
        if id(v) == elem_id
            return v
        end
    end
    error("No element of name $elem_id")
end

struct Hole
    id::String
    tees::Vector{POI}
    pins::Vector{POI}
    pars::Dict{String,Dict{String, Int}}
end
empty_hole() = Hole("",[],[],Dict(""=>Dict(""=>0)))
tees(h::Hole) = h.tees
pins(h::Hole) = h.pins
id(h::Hole) = h.id
pars(h::Hole) = h.pars
par(h::Hole,tee::String,pin::String) = pars(h)[tee][pin]
par(h::Hole,tee::POI,pin::POI) = pars(h)[id(tee)][id(pin)]
get_pin(h::Hole, id::String) = get_elem(id,pins(h))
get_tee(h::Hole, id::String) = get_elem(id,tees(h))
get_pin(h::Hole, p::POI) = get_elem(id(p),pins(h))
get_tee(h::Hole, t::POI) = get_elem(id(t),tees(h))

struct Course
    id::String
    name::String
    loc::Location
    holes::Vector{Hole}
end
loc(c::Course) = c.loc
id(c::Course) = c.id
holes(c::Course) = c.holes
name(c::Course) = c.name
lat(c::Course) = lat(loc(c))
lon(c::Course) = lon(loc(c))
tees(c::Course) = tees.(holes(c))
get_hole(c::Course, id::String) = get_elem(id,holes(c))

struct Throw
    loc1::Location
    loc2::Location
    type::String
    res::String
    pin::POI
    disc::String
end
loc1(t::Throw) = t.loc1
loc2(t::Throw) = t.loc2
lat(t::Throw) = lat(loc1(t))
lon(t::Throw) = lon(loc1(t))
disc(t::Throw) = t.disc

struct PlayedHole
    hole::Hole
    tee::POI
    pin::POI
    throws::Vector{Throw}
    res::Int
end
throws(ph::PlayedHole) = ph.throws
pin(ph::PlayedHole) = ph.pin
tee(ph::PlayedHole) = ph.tee
id(ph::PlayedHole) = id(ph.hole)
score(ph::PlayedHole) = length(throws(ph))
res(ph::PlayedHole) = ph.res
par(ph::PlayedHole) = par(ph.hole, ph.tee, ph.pin)
function res_name(ph::PlayedHole)
    if haskey(SCORE_NAMES, res(ph))
        return SCORE_NAMES[res(ph)]
    else
        return "BAD"
    end
end

struct Round
    id::String
    name::String
    notes::String
    start_time::DateTime
    end_time::DateTime
    course::Course
    holes::Vector{PlayedHole}
end
course(r::Round) = r.course
id(r::Round) = r.id
name(r::Round) = r.name
notes(r::Round) = r.notes
holes(r::Round) = r.holes
date(r::Round) = Dates.format(r.start_time,"yyyy-mm-dd")
score(r::Round) = sum([res(h) for h in holes(r)])
num_holes(r::Round) = length(holes(r))
par(r::Round) = sum([par(h) for h in holes(r)])
total(r::Round) = par(r) + score(r)

#get_scorecard(r::Round)
#make_vis(r::Round)

# function read_round(input_csv, label, notes)

#     Round(id, label, time, course, holes, notes)
# end

# function make_preview(r::Round, preview_file)
#     return -1
# end

# function dg_preview(input_csv; preview_file = "preview.kml")

#     round = make_round(input_csv)
#     make_preview(round, preview_file)

# end

function read_courses_dir(courses_dir = COURSES_DIR)
    files = readdir(courses_dir)

    C = Course[]
    for f in files
        d = JSON.parsefile(joinpath(courses_dir, f), dicttype=OrderedDict)
        c = process_course_dict(d)
        push!(C,c)
    end

    C
end

function process_course_dict(d)
    id = d["id"]
    name = d["name"]
    loc = Location(d["loc"])

    holes = Hole[]
    for (k,v) in d["holes"]
        h_id = String(k)

        pins = POI[]
        for (kk,vv) in v["pins"]
            push!(pins, POI(String(kk), Location(vv)))
        end

        tees = POI[]
        for (kk,vv) in v["tees"]
            push!(tees, POI(String(kk), Location(vv)))
        end

        ###make pars
        pars = Dict{String,typeof(Dict{String, Int}())}()
        for (kk,vv) in v["pars"]
            ppars = Dict{String, Int}()
            for (kkk,vvv) in vv
                ppars[kkk] = vvv
            end
            pars[kk] = ppars
        end

        push!(holes, Hole(h_id, tees, pins, pars))

    end

    Course(id, name, loc, holes)
end

function infer_course(l::Location, courses::Vector{Course})
    best_dist = 1000
    best_course = courses[1]

    for c in courses
        this_dist = dist(l, loc(c))
        if this_dist < best_dist
            best_dist = this_dist
            best_course = c
        end
    end

    if best_dist > 999
        error("No courses found within 1 km.")
    end

    return best_course

end

function infer_tee(l::Location, c::Course)
    best_dist = 100
    best_hole = empty_hole()    #placeholder
    best_tee = empty_poi()      #placeholder
    for h in holes(c)
        for t in tees(h)
            this_dist = dist(l,loc(t))
            if this_dist < best_dist
                best_dist = this_dist
                best_hole = h
                best_tee = t
            end
        end
    end
    best_hole, best_tee, best_dist
end

function parse_round_raw_csv(round_raw_csv)

    COURSES = read_courses_dir()

    df = CSV.read(round_raw_csv, DataFrame)
    nrows = size(df,1)

    t_start = df[1,:time]
    t_end = df[end,:time]

    loc_i = Location(df[1,:lat], df[1,:lon])
    c = infer_course(loc_i, COURSES)

    PH = PlayedHole[]

    cutoff_dists = [10,10,20] #this is placeholder - should be a better way to do this

    #When we are just starting up, let's snag the teebox and the hole we are playing
    cur_i = 1
    h, t, _ = infer_tee(loc_i, c)

    completed_read = false

    # while !eof
    while ~completed_read

        #Now we need to find where the hole ends, this will happen when we mark near a pin of the current hole, and then near a tee of another hole, and then not near a tee of that hole
        p = empty_poi() #to overwrite
        next_t = empty_poi()
        next_h = empty_hole()
        idx = 0

        for i in cur_i:nrows-1

            r = df[i,:]
            this_loc = Location(r.lat, r.lon)
            pds = [dist(loc(p), this_loc) for p in pins(h)]
            pd, p_num = findmin(pds)
            p = pins(h)[p_num]

            if i==nrows-1
                idx = i+1
                completed_read = true
                break
            end

            n_r = df[i+1,:]
            n_loc = Location(n_r.lat, n_r.lon)
            next_h, next_t, bd = infer_tee(n_loc, c)

            nn_r = df[i+2,:]
            nn_loc = Location(nn_r.lat, nn_r.lon)
            ndd = dist(nn_loc, loc(next_t))

            if pd<cutoff_dists[1] && bd<cutoff_dists[2] && ndd>cutoff_dists[3] #then we have found the last putt of the hole!
                idx = i
                break
            end
        end

        #Now we know what throws were made on the hole, which teebox, and which pin. We can make some Throws and a PlayedHole!
        i1 = cur_i
        i2 = idx-1
        T = Throw[]
        for i in i1:i2
            loc1 = Location(df[i,:lat], df[i,:lon])
            loc2 = Location(df[i+1,:lat], df[i+1,:lon])
            if i==i1
                type = "DRIVE"
                loc1 = loc(t)
            else
                type = "THROW"
            end

            if i==i2
                res = "BASKET"
                loc2 = loc(p)
            else
                res = ""
            end

            pin = p
            disc = df[i,:disc]

            thrw = Throw(loc1, loc2, type, res, pin, disc)

            push!(T, thrw)
        end

        s = length(T) - par(h, t, p)
        ph = PlayedHole(h,t,p,T,s)
        push!(PH, ph)

        t = next_t
        h = next_h
        cur_i = idx+1

    end

    start_minute = Dates.format(t_start, "yyyy-mm-dd-HH-MM")
    round_id = start_minute * "_" * id(c)
    name = "Sample Name"
    notes = "Sample Notes"
    start_time = t_start
    end_time = t_end

    Round(round_id, name, notes, start_time, end_time, c, PH)
end

function guess_and_check(round_raw_csv, out="check_round.json")

    r = parse_round_raw_csv(round_raw_csv)

    df_raw = CSV.read(round_raw_csv, DataFrame)

    json_both = draw_both(r,df_raw)

    open(out, "w") do f
        JSON.print(f, json_both, 2)
    end

end

function draw_both(r::Round, df)
    D = OrderedDict()
    D["type"] = "FeatureCollection"

    #Some metadata up front
    my_discs = CSV.read(joinpath(STATS_DIR, "discs", "discs.csv"), DataFrame)
    discs = OrderedDict[]
    for r in eachrow(my_discs)
        push!(discs, OrderedDict(
            "id" => r.my_id,
            "image" => r.image
        ))
    end
    D["discs"] = discs

    summary = []
    running = 0
    for h in holes(r)
        running += res(h)
        h_summary = OrderedDict(
            "hole" => id(h),
            "tee" => id(tee(h)),
            "pin" => id(pin(h)),
            "par" => par(h),
            "score" => score(h),
            "running" => running
        )
        push!(summary, h_summary)
    end
    D["summary"] = summary

    features = OrderedDict[]

    #Now let's put all the raw marks on the map
    for row in eachrow(df)
        d = OrderedDict()
        d["type"] = "Feature"
        d["properties"] = OrderedDict(
            "thing" => "raw_mark",
            "disc_name" => row.disc,
            "name"  => "$(rownumber(row))"
            )
        d["geometry"] = OrderedDict(
                "type" => "Point",
                "coordinates" => [row.lon, row.lat]
            )
            push!(features, d)
    end

    #Now draw the course
    c = course(r)
    for h in holes(c)
        for t in tees(h)
            d = OrderedDict()
            d["type"] = "Feature"
            d["properties"] = OrderedDict(
                "thing" => "tee",
                "name"  => id(h) * id(t)
                )
            d["geometry"] = OrderedDict(
                "type" => "Point",
                "coordinates" => [lon(t), lat(t)]
            )
            push!(features, d)
        end

        for p in pins(h)
            d = OrderedDict()
            d["type"] = "Feature"
            d["properties"] = OrderedDict(
                "thing" => "pin",
                "name" => id(h) * id(p))
            d["geometry"] = OrderedDict(
                "type" => "Point",
                "coordinates" => [lon(p), lat(p)]
            )
            push!(features, d)
        end
    end

    #Now draw the throws as they were interpreted
    for h in holes(r)
        for t in h.throws
            d = OrderedDict()
            d["type"] = "Feature"
            d["properties"] = OrderedDict(
                "thing" => "throw",
                "hole_res" => h.res)
            d["geometry"] = OrderedDict(
                "type" => "LineString",
                "coordinates" => [[lon(loc1(t)), lat(loc1(t))],[lon(loc2(t)), lat(loc2(t))]]
            )
            push!(features, d)
        end
    end

    D["features"] = features

    D
end

feet(x) = round(x*3.2808)
llstring(x) = string(round(x, digits=5))

function append_all_throws!(df,r::Round)
    for h in holes(r)
        shot=0

        for t in throws(h)
            shot += 1

            if score(h) == 1
                t_res = "ACE"
            elseif shot == 1
                t_res = "DRIVE"
            elseif shot == score(h)
                t_res = "MAKE"
            else
                t_res = "THROW"
            end

            push!(df, (
                        date(r) |> string,
                        id(course(r)),
                        id(r),
                        id(h),
                        shot |> string,
                        disc(t),
                        lat(loc1(t)) |> llstring,
                        lon(loc1(t)) |> llstring,
                        lat(loc2(t)) |> llstring,
                        lon(loc2(t)) |> llstring,
                        lat(pin(h)) |> llstring,
                        lon(pin(h)) |> llstring,
                        dist(loc1(t), loc2(t)) |> feet |> string,
                        dist(loc1(t), loc(pin(h))) |> feet |> string,
                        res_name(h),
                        t_res)
                )
        end
    end
end

function append_all_holes!(df, r::Round)
    for h in holes(r)
        push!(df, (
                    date(r) |> string,
                    id(course(r)),
                    id(r),
                    id(h),
                    par(h) |> string,
                    id(tee(h)),
                    id(pin(h)),
                    dist(tee(h), pin(h)) |> feet |> string,
                    dist(loc1(throws(h)[1]), loc2(throws(h)[1])) |> feet |> string,
                    dist(loc1(throws(h)[end]), loc2(throws(h)[end])) |> feet |> string,
                    score(h) |> string,
                    res_name(h))
            )
    end
end

function append_round!(df, r::Round)
    push!(df, (
        date(r) |> string,
        id(course(r)),
        id(r),
        num_holes(r) |> string,
        par(r) |> string,
        total(r) |> string,
        score(r) |> string
    ))
end

function get_dash_json(throws_df, holes_df, rounds_df)
    D = OrderedDict()

    #Recent Scores
    recent_scores = OrderedDict()
    recent_scores["labels"] = rounds_df.date
    recent_scores["datasets"] = [OrderedDict(
        "label" => "Recent Scores",
        "data" => rounds_df.result,
        "backgroundColor" => "rgba(60,141,188,0.9)",
        "borderColor" => "rgba(60,141,188,0.8)",
        "pointRadius" => 5,
        "pointColor" => "#3b8bba",
        "pointStrokeColor" => "rgba(60,141,188,1)",
        "pointHighlightFill" => "#fff",
        "pointHighlightStroke" => "rgba(60,141,188,1)",
        "showLine" => true,
        "fill" => false
    )]
    D["recent_scores"] = recent_scores

    #Disc Bag
    my_discs = CSV.read(joinpath(STATS_DIR, "discs", "discs.csv"), DataFrame)
    discs_db = CSV.read(DISCS_CSV, DataFrame, type=String)
    all_discs = innerjoin(my_discs, discs_db, on = :disc_id => :id)
    discs = OrderedDict[]
    for r in eachrow(all_discs)
        push!(discs, OrderedDict(
            "id" => r.my_id,
            "image" => r.image,
            "mold" => "$(r.brand) $(r.mold)",
            "plastic" => r.plastic,
            "weight" => r.weight,
            "numbers" => "$(r.speed), $(r.glide), $(r.turn), $(r.fade)"
        ))
    end
    D["discs"] = discs

    #All Rounds
    rounds = OrderedDict[]
    COURSES = read_courses_dir()
    for r in eachrow(rounds_df)
        c_name = "Unknown"
        for c in COURSES
            if r.course == id(c)
                c_name = name(c)
                break
            end
        end
        push!(rounds, OrderedDict(
            "id" => r.round,
            "date" => r.date,
            "course" => c_name,
            "result" => r.result
        ))
    end
    D["rounds"] = rounds

    D
end

function save_round(round_raw_csv)

    r = parse_round_raw_csv(round_raw_csv)

    #write out round.json

    #append all throws to all_throws.csv
    throws_csv = joinpath(STATS_DIR,"stats","all_throws.csv")
    throws_df = CSV.read(throws_csv, DataFrame, type=String)
    append_all_throws!(throws_df,r)
    CSV.write(throws_csv, throws_df)

    #append all holes to all_holes.csv
    holes_csv = joinpath(STATS_DIR,"stats","all_holes.csv")
    holes_df = CSV.read(holes_csv, DataFrame, type=String)
    append_all_holes!(holes_df,r)
    CSV.write(holes_csv, holes_df)

    #append round to all_rounds.csv
    rounds_csv = joinpath(STATS_DIR,"stats","all_rounds.csv")
    rounds_df = CSV.read(rounds_csv, DataFrame, type=String)
    append_round!(rounds_df,r)
    CSV.write(rounds_csv, rounds_df)

    #write dash.json with updated data
    dash_json = joinpath(STATS_DIR,"stats","dash.json")
    dash = get_dash_json(throws_df, holes_df, rounds_df)
    open(dash_json, "w") do f
        JSON.print(f, dash, 2)
    end

end

end


text:
@testset "ResidueField" begin

  Qx, x = FlintQQ["x"]
  f = x^6 + x^5 + 41*x^4 - 34*x^3 + 355*x^2 - 100*x + 125
  K, a = number_field(f, cached = false);
  OK = maximal_order(K);
  lp = prime_decomposition(OK, 5)
  B = basis(OK)
  for i = 1:length(lp)
    P = lp[i][1]
    F, mF = ResidueField(OK, P)
    F1, mF1 = Hecke.ResidueFieldSmall(OK, P)
    F2, mF2 = Hecke.ResidueFieldSmallDegree1(OK, P)
    F3, mF3 = Hecke.ResidueFieldDegree1(OK, P)
    @test order(F1) == order(F2)
    @test order(F1) == order(F3)
    @test order(F1) == order(F)

    for b in B
      a1 = mF(b)
      a2 = mF1(b)
      a3 = mF2(b)
      a4 = mF3(b)
      @test (mF\(a1) - b) in P
      @test (mF1\(a2) - b) in P
      @test (mF2\(a3) - b) in P
      @test (mF3\(a4) - b) in P
    end
  end

  lp = prime_decomposition(OK, 47)
  for i = 1:length(lp)
    P = lp[i][1]
    F, mF = ResidueField(OK, P)
    F1, mF1 = Hecke.ResidueFieldSmall(OK, P)
    F2, mF2 = Hecke.ResidueFieldSmallDegree1(OK, P)
    F3, mF3 = Hecke.ResidueFieldDegree1(OK, P)
    @test order(F1) == order(F2)
    @test order(F1) == order(F3)
    @test order(F1) == order(F)

    for b in B
      a1 = mF(b)
      a2 = mF1(b)
      a3 = mF2(b)
      a4 = mF3(b)
      @test (mF\(a1) - b) in P
      @test (mF1\(a2) - b) in P
      @test (mF2\(a3) - b) in P
      @test (mF3\(a4) - b) in P
    end
  end
end


text:
module HistogramThresholding

using Images
using LinearAlgebra

abstract type HistogramDistribution end
struct Unimodal <: HistogramDistribution end
struct Bimodal <: HistogramDistribution end

abstract type ThresholdAlgorithm end
struct Otsu <: ThresholdAlgorithm end
struct Moments <: ThresholdAlgorithm end

include("otsu.jl")
include("moments.jl")

export
	# main functions
    find_threshold,
	Otsu,
	Moments
end # module

max_stars_repo_name: AbbeyMac/HistogramThresholding

max_issues_repo_name: AbbeyMac/HistogramThresholding

max_forks_repo_name: AbbeyMac/HistogramThresholding


text:
import PackageCompiler

_default_docker_parent_image() = "ubuntu:latest"

_default_julia_version() = VersionNumber(VERSION.major, VERSION.minor, VERSION.patch)

_default_packagecompiler_installation_command() = "import Pkg; Pkg.add(Pkg.PackageSpec(name = \"PackageCompiler\", version = \"1.2.2 - 1\"));"

_default_simplecontainergenerator_installation_command() = "import Pkg; Pkg.add(Pkg.PackageSpec(name = \"SimpleContainerGenerator\"));"

_default_julia_cpu_target() = "$(PackageCompiler.default_app_cpu_target());$(PackageCompiler.NATIVE_CPU_TARGET)"
_recommended_julia_cpu_target() = "$(PackageCompiler.default_app_cpu_target());$(PackageCompiler.NATIVE_CPU_TARGET)"
_fastest_nonportable_julia_cpu_target() = "$(PackageCompiler.NATIVE_CPU_TARGET)"

 _default_precompile_execution_env_vars() = Dict{String, String}(
    "PREDICTMD_TEST_GROUP" => "all",
    "PREDICTMD_TEST_PLOTS" => "true"
)

_default_wrapper_script_env_vars() = Dict{String, String}(
    "PREDICTMD_TEST_GROUP" => "all",
    "PREDICTMD_TEST_PLOTS" => "true",
    "GKSwstype" => "100"
)

_predictmd_apt() = String[]

_default_apt() = String[
    "apt-utils",
    "build-essential",
    "bzip2",
    "ca-certificates",
    "cmake",
    "coreutils",
    "curl",
    "emacs ",
    "emacs nano vim",
    "fonts-liberation",
    "gettext",
    "gfortran",
    "git",
    "git-all",
    "git-flow",
    "git-lfs",
    "gnupg",
    "gpg",
    "gpg-agent",
    "gzip",
    "hdf5-tools",
    "libcurl4-openssl-dev",
    "libgconf-2-4",
    "libgtk2.0-0",
    "libnss3",
    "libpgf-dev",
    "libpgf6",
    "libpgf6-dbg",
    "libpng-dev",
    "libssl-dev",
    "libxss1",
    "libxtst6",
    "locales",
    "lsb-release",
    "m4",
    "nano",
    "openssh-client",
    "openssl",
    "pdf2svg",
    "poppler-utils",
    "qt5-default",
    "screen",
    "sudo",
    "texlive-binaries",
    "texlive-latex-base",
    "texlive-latex-extra",
    "texlive-luatex",
    "texlive-pictures",
    "tmux ",
    "tree",
    "unzip",
    "vim",
    "wget",
    "xdg-utils",
    "zip",
    "zlib1g-dev",
]


text:
if VERSION < v"1.7.0-DEV.849"
    @error "" VERSION < v"1.7.0-DEV.849"
    error()
end

include("generator_utils.jl")

using REPL
write_dict_to_file("repl_emoji_symbols_170.jl",
                   "repl_emoji_symbols_170",
                   REPL.REPLCompletions.emoji_symbols)


text:
using Test, OMEinsum
using ForwardDiff
using FiniteDifferences
using TensorNetworkEvolve.TensorAD

include("ADTest.jl")
using .ADTest: match_jacobian, match_random, match_hessian
include("rules.jl")

@testset "diff tensor" begin
    function f(x, y)
        z = ein"ij,jk->ik"(x, y)
        return ein"ii->"(sin.(cos.(cos.(z))) .* x .* y)
    end
    x = randn(10, 10)
    y = randn(10, 10)
    gs = TensorAD.gradient(f, DiffTensor(x), DiffTensor(y))
    hs = ForwardDiff.gradient(x->f(reshape(x[1:100], 10, 10), reshape(x[101:200], 10, 10))[], vcat(vec(x), vec(y)))
    @test vcat(vec(gs[1].data), vec(gs[2].data)) ≈ hs
end

@testset "hessian" begin
    function f1(x, y)
        z = ein"ij,jk->ik"(x, y)
        return ein"ii->"(z)
    end
    @test match_hessian(f1, randn(5, 5), randn(5, 5))
    function f2(xy)
        sin.(sin.(sin.(xy)))
    end
    @test match_hessian(f2, fill(0.5))
    function f3(xy)
        x = xy[1:1]
        y = xy[1:1]
        z = xy[1:1]
        ein"i->"((x .* y) .* z)
    end
    @test match_hessian(f3, [0.2, 0.4, 0.8])
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# -------------
# 2D ROTATIONS
# -------------

# temporary fix for https://github.com/JuliaSpace/ReferenceFrameRotations.jl/issues/22

"""
    ClockwiseAngle(θ)

Clockwise rotation in 2D space by angle `θ`.
"""
struct ClockwiseAngle{T}
  θ::T
end

Base.inv(cw::ClockwiseAngle) = CounterClockwiseAngle(cw.θ)

function Base.convert(::Type{DCM{T}}, cw::ClockwiseAngle) where {T}
  s, c = sincos(cw.θ)
  SMatrix{2,2,T}([c s; -s c])
end

"""
    ClockwiseAngle(θ)

Counter-clockwise rotation in 2D space by angle `θ`.
"""
struct CounterClockwiseAngle{T}
  θ::T
end

Base.inv(ccw::CounterClockwiseAngle) = ClockwiseAngle(ccw.θ)

Base.convert(::Type{<:DCM}, ccw::CounterClockwiseAngle) =
  convert(DCM, ClockwiseAngle(-ccw.θ))

# -------------
# 3D ROTATIONS
# -------------

TaitBryanAngles(θs...) = EulerAngles((-1 .* θs)..., :ZXY) |> inv


text:
const opencl_builtins = String["printf"]

# OpenCL functions need to be mangled according to the C++ Itanium spec. We implement a very
# limited version of that spec here, just enough to support OpenCL built-ins.
#
# This macro also keeps track of called builtins, generating `ccall("extern...", llvmcall)`
# expressions for them (so that we can exclude them during IR verification).
macro builtin_ccall(name, ret, argtypes, args...)
    @assert Meta.isexpr(argtypes, :tuple)
    argtypes = argtypes.args

    function mangle(T::Type)
        if T == Cint
            "i"
        elseif T == Cuint
            "j"
        elseif T == Clong
            "l"
        elseif T == Culong
            "m"
        elseif T == Cshort
            "s"
        elseif T == Cushort
            "t"
        elseif T == Cchar
            "c"
        elseif T == Cuchar
            "h"
        elseif T == Cfloat
            "f"
        elseif T == Cdouble
            "d"
        elseif T <: LLVMPtr
            elt, as = T.parameters

            # mangle address space
            ASstr = if as == AS.Global
                "CLglobal"
            #elseif as == AS.Global_device
            #    "CLdevice"
            #elseif as == AS.Global_host
            #    "CLhost"
            elseif as == AS.Local
                "CLlocal"
            elseif as == AS.Constant
                "CLconstant"
            elseif as == AS.Private
                "CLprivate"
            elseif as == AS.Generic
                "CLgeneric"
            else
                error("Unknown address space $AS")
            end

            # encode as vendor qualifier
            ASstr = "U" * string(length(ASstr)) * ASstr

            # XXX: where does the V come from?
            "P" * ASstr * "V" * mangle(elt)
        else
            error("Unknown type $t")
        end
    end

    # C++-style mangling; very limited to just support these intrinsics
    # TODO: generalize for use with other intrinsics? do we need to mangle those?
    mangled = "_Z$(length(name))$name"
    for t in argtypes
        # with `@eval @builtin_ccall`, we get actual types in the ast, otherwise symbols
        t = (isa(t, Symbol) || isa(t, Expr)) ? eval(t) : t
        mangled *= mangle(t)
    end

    push!(opencl_builtins, mangled)
    esc(quote
        @typed_ccall($mangled, llvmcall, $ret, ($(argtypes...),), $(args...))
    end)
end

## device overrides

# local method table for device functions
@static if isdefined(Base.Experimental, Symbol("@overlay"))
Base.Experimental.@MethodTable(method_table)
else
const method_table = nothing
end

# list of overrides (only for Julia 1.6)
const overrides = quote end

macro device_override(ex)
    code = quote
        $GPUCompiler.@override($method_table, $ex)
    end
    if isdefined(Base.Experimental, Symbol("@overlay"))
        return esc(code)
    else
        push!(overrides.args, code)
        return
    end
end

macro device_function(ex)
    ex = macroexpand(__module__, ex)
    def = ExprTools.splitdef(ex)

    # generate a function that errors
    def[:body] = quote
        error("This function is not intended for use on the CPU")
    end

    esc(quote
        $(ExprTools.combinedef(def))
        @device_override $ex
    end)
end


text:
baremodule Config
    julia_path = "/Applications/Julia-0.7.app/Contents/Resources/julia/bin/julia"

end


text:
"Convert Cartesian to Spherical coordinates"
function cart2pol(x::Real, y::Real, z::Real)
    ρ = sqrt(x^2 + y^2 + z^2)  # distance
    θ = atand(y, x)  # horizontal angle
    ϕ = acosd(z / ρ)  # incidence angle
    (ρ, θ, ϕ)
end

"Convert Spherical to Cartesian coordinates"
function pol2cart(ρ::Real, θ::Real, ϕ::Real)
    sθ, cθ = sincosd(θ)
    sϕ, cϕ = sincosd(ϕ)
    x = ρ * sϕ * cθ
    y = ρ * sϕ * sθ
    z = ρ * cϕ
    (x, y, z)
end


text:
include("poisson.jl")
include("negativebinomial.jl")

const EventLikelihood = Union{PoissonLikelihood,NegBinomialLikelihood}

""" Return the labels in a vector of vectors for multiple outputs"""
function treat_labels!(y::AbstractVector{<:Int}, ::EventLikelihood)
    return y
end

function treat_labels!(::AbstractVector{<:Real}, ::EventLikelihood)
    return error("For event count target(s) should be integers")
end


text:
using EngEconomics, Roots

# Given
aInit = -1500000
aAnnual = 900000 - 600000
aScrap = 100000
aLife = 5

bInit = -2000000
bAnnual = 1100000 - 800000
bScrap = 200000
bLife = 10

MARR = 0.08

# Find Annual Worth
aAW = aInit * capitalRecoveryFactor(MARR, aLife) + aAnnual + aScrap * sinkingFundFactor(MARR, aLife)
bAW = bInit * capitalRecoveryFactor(MARR, bLife) + bAnnual + bScrap * sinkingFundFactor(MARR, bLife)

if aAW >= 0 && bAW >= 0
	println("Both")
elseif aAW >= 0
	println("A")
elseif bAW >= 0
	println("B")
else
	println("Neither")
end

# Find Present Worth
aPW = aInit + aAnnual * seriesPresentAmountFactor(MARR, aLife) + aScrap * presentWorthFactor(MARR, aLife)
bPW = bInit + bAnnual * seriesPresentAmountFactor(MARR, bLife) + bScrap * presentWorthFactor(MARR, bLife)

if aPW >= 0 && bPW >= 0
	println("Both")
elseif aPW >= 0
	println("A")
elseif bPW >= 0
	println("B")
else
	println("Neither")
end

# Payback Period
# - FirstCost / AnnualSavings
aPP = -aInit / aAnnual
bPP = -bInit / bAnnual
minPP = 3

if aPP <= minPP && bPP <= minPP
	println("Both")
elseif aPP <= minPP
	println("A")
elseif bPP <= minPP
	println("B")
else
	println("Neither")
end


text:
#!/usr/bin/julia
print("List Iteration\n")
l = ["geeks", "for", "geeks"]
for i in l
    println(i)
end

# Iterating over a tuple (immutable)
print("\nTuple Iteration\n")
t = ("geeks", "for", "geeks")
for i in t
    println(i)
end

# Iterating over a String
print("\nString Iteration\n")
s = "Geeks"
for i in s
    println(i)
end


text:
#Define AbstractTransformable type for transformation forces/moments/etc.
abstract type AbstractTransformable end


text:
@system Plant(
    Mass,
    Area,
    Count,
    Ratio,
    Carbon,
    #Nitrogen,
    Water,
    Weight,
    Density,
    Photosynthesis
) begin
    calendar(context) ~ ::Calendar
    weather(context, calendar) ~ ::Weather
    sun(context, calendar, weather) ~ ::Sun
    soil(context) ~ ::Soil
    pheno(context, calendar, weather, sun, soil): phenology ~ ::Phenology

    primordia => 5 ~ preserve::int(parameter)

    #bulb => begin end ~ produce::Bulb

    #scape => begin end ~ produce::Scape

    root(pheno, emerging=pheno.emerging) => begin
        if emerging
            #TODO import_carbohydrate(soil.total_root_weight)
            produce(Root, phenology=pheno)
        end
    end ~ produce::Root

    #TODO pass PRIMORDIA as initial_leaves
    NU(NU, pheno, primordia, germinated=pheno.germinated, dead=pheno.dead, l=pheno.leaves_initiated): nodal_units => begin
        if germinated && !dead
            n = length(NU)
            if n == 0
                [produce(NodalUnit, phenology=pheno, rank=i) for i in 1:primordia]
            elseif n < l
                [produce(NodalUnit, phenology=pheno, rank=i) for i in (n+1):l]
            end
        end
    end ~ produce::NodalUnit[]

    # dry_yield(bulb_mass, PD) => bulb_mass * PD ~ track(u"g/m^2")
    # fresh_yield(dry_yield, BMC) => begin
    #     dry_yield * 1 / (1 - BMC)
    # end ~ track(u"g/m^2")
    # BMC: bulb_moisture_content => 0.85 ~ preserve

    SLA(leaf_area, leaf_mass): specific_lear_area => begin
        leaf_area / leaf_mass
    end ~ track(u"cm^2/g")
    LMA(SLA): leaf_mass_per_area => 1 / SLA ~ track(u"g/cm^2")

    DAP(pheno.DAP): day_after_planting ~ track::int(u"d")
    time(calendar.time) ~ track::datetime
    date(calendar.date) ~ track::date
end

@system Model(Plant, Controller)


text:
module ElementSymbols

export element_symbol, zahl

const element_symbols = [
    "H", "He", "Li", "Be",  "B",  "C",  "N",  "O",  "F", "Ne",
   "Na", "Mg", "Al", "Si",  "P",  "S", "Cl", "Ar",  "K", "Ca",
   "Sc", "Ti",  "V", "Cr", "Mn", "Fe", "Co", "Ni", "Cu", "Zn",
   "Ga", "Ge", "As", "Se", "Br", "Kr", "Rb", "Sr",  "Y", "Zr",
   "Nb", "Mo", "Tc", "Ru", "Rh", "Pd", "Ag", "Cd", "In", "Sn",
   "Sb", "Te",  "I", "Xe", "Cs", "Ba", "La", "Ce", "Pr", "Nd",
   "Pm", "Sm", "Eu", "Gd", "Tb", "Dy", "Ho", "Er", "Tm", "Yb",
   "Lu", "Hf", "Ta",  "W", "Re", "Os", "Ir", "Pt", "Au", "Hg",
   "Tl", "Pb", "Bi", "Po", "At", "Rn", "Fr", "Ra", "Ac", "Th",
   "Pa",  "U", "Np", "Pu", "Am", "Cm", "Bk", "Cf", "Es", "Fm",
   "Md", "No", "Lr", "Rf", "Db", "Sg", "Bh", "Hs", "Mt", "Ds",
   "Rg", "Cn", "Nh", "Fl", "Mc", "Lv", "Ts", "Og"
   ]

function element_symbol(zahl::Integer)
    if zahl == 0
        return "n"
    elseif zahl > 0 && zahl <= length(element_symbols)
        return element_symbols[zahl]
    else
        return "X"
    end
end

zahl(symbol::AbstractString) = indexin([symbol], element_symbols)[begin]

end # module


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

## Create an extractor that extracts the modified original matrix, e.g.
## LD for BunchKaufman, UL for CholeskyDense, LU for LUDense and
## define size methods for Factorization types using it.

struct BunchKaufman{T,S<:AbstractMatrix} <: Factorization{T}
    LD::S
    ipiv::Vector{BlasInt}
    uplo::Char
    symmetric::Bool
    rook::Bool
    info::BlasInt
end
BunchKaufman(A::AbstractMatrix{T}, ipiv::Vector{BlasInt}, uplo::Char, symmetric::Bool,
             rook::Bool, info::BlasInt) where {T} =
        BunchKaufman{T,typeof(A)}(A, ipiv, uplo, symmetric, rook, info)

"""
    bkfact!(A, uplo::Symbol=:U, symmetric::Bool=issymmetric(A), rook::Bool=false) -> BunchKaufman

`bkfact!` is the same as [`bkfact`](@ref), but saves space by overwriting the
input `A`, instead of creating a copy.
"""
function bkfact!(A::StridedMatrix{<:BlasReal}, uplo::Symbol = :U,
    symmetric::Bool = issymmetric(A), rook::Bool = false)

    if !symmetric
        throw(ArgumentError("Bunch-Kaufman decomposition is only valid for symmetric matrices"))
    end
    if rook
        LD, ipiv, info = LAPACK.sytrf_rook!(char_uplo(uplo), A)
    else
        LD, ipiv, info = LAPACK.sytrf!(char_uplo(uplo), A)
    end
    BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric, rook, info)
end
function bkfact!(A::StridedMatrix{<:BlasComplex}, uplo::Symbol=:U,
    symmetric::Bool=issymmetric(A), rook::Bool=false)

    if rook
        if symmetric
            LD, ipiv, info = LAPACK.sytrf_rook!(char_uplo(uplo), A)
        else
            LD, ipiv, info = LAPACK.hetrf_rook!(char_uplo(uplo), A)
        end
    else
        if symmetric
            LD, ipiv, info = LAPACK.sytrf!(char_uplo(uplo),  A)
        else
            LD, ipiv, info = LAPACK.hetrf!(char_uplo(uplo), A)
        end
    end
    BunchKaufman(LD, ipiv, char_uplo(uplo), symmetric, rook, info)
end

"""
    bkfact(A, uplo::Symbol=:U, symmetric::Bool=issymmetric(A), rook::Bool=false) -> BunchKaufman

Compute the Bunch-Kaufman [^Bunch1977] factorization of a symmetric or Hermitian
matrix `A` and return a `BunchKaufman` object.
`uplo` indicates which triangle of matrix `A` to reference.
If `symmetric` is `true`, `A` is assumed to be symmetric. If `symmetric` is `false`,
`A` is assumed to be Hermitian. If `rook` is `true`, rook pivoting is used. If
`rook` is false, rook pivoting is not used.
The following functions are available for
`BunchKaufman` objects: [`size`](@ref), `\\`, [`inv`](@ref), [`issymmetric`](@ref), [`ishermitian`](@ref).

[^Bunch1977]: J R Bunch and L Kaufman, Some stable methods for calculating inertia and solving symmetric linear systems, Mathematics of Computation 31:137 (1977), 163-179. [url](http://www.ams.org/journals/mcom/1977-31-137/S0025-5718-1977-0428694-0/).

"""
bkfact(A::StridedMatrix{<:BlasFloat}, uplo::Symbol=:U, symmetric::Bool=issymmetric(A),
    rook::Bool=false) =
        bkfact!(copy(A), uplo, symmetric, rook)
bkfact(A::StridedMatrix{T}, uplo::Symbol=:U, symmetric::Bool=issymmetric(A),
    rook::Bool=false) where {T} =
        bkfact!(convert(Matrix{promote_type(Float32, typeof(sqrt(one(T))))}, A),
                uplo, symmetric, rook)

convert(::Type{BunchKaufman{T}}, B::BunchKaufman{T}) where {T} = B
convert(::Type{BunchKaufman{T}}, B::BunchKaufman) where {T} =
    BunchKaufman(convert(Matrix{T}, B.LD), B.ipiv, B.uplo, B.symmetric, B.rook, B.info)
convert(::Type{Factorization{T}}, B::BunchKaufman{T}) where {T} = B
convert(::Type{Factorization{T}}, B::BunchKaufman) where {T} = convert(BunchKaufman{T}, B)

size(B::BunchKaufman) = size(B.LD)
size(B::BunchKaufman, d::Integer) = size(B.LD, d)
issymmetric(B::BunchKaufman) = B.symmetric
ishermitian(B::BunchKaufman) = !B.symmetric

issuccess(B::BunchKaufman) = B.info == 0

function Base.show(io::IO, F::BunchKaufman)
    println(io, summary(F))
    print(io, "successful: $(issuccess(F))")
end

function inv(B::BunchKaufman{<:BlasReal})
    if !issuccess(B)
        throw(SingularException(B.info))
    end

    if B.rook
        copytri!(LAPACK.sytri_rook!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
    else
        copytri!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
    end
end

function inv(B::BunchKaufman{<:BlasComplex})
    if !issuccess(B)
        throw(SingularException(B.info))
    end

    if issymmetric(B)
        if B.rook
            copytri!(LAPACK.sytri_rook!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
        else
            copytri!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
        end
    else
        if B.rook
            copytri!(LAPACK.hetri_rook!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
        else
            copytri!(LAPACK.hetri!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
        end
    end
end

function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasReal
    if !issuccess(B)
        throw(SingularException(B.info))
    end

    if B.rook
        LAPACK.sytrs_rook!(B.uplo, B.LD, B.ipiv, R)
    else
        LAPACK.sytrs!(B.uplo, B.LD, B.ipiv, R)
    end
end
function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{T}) where T<:BlasComplex
    if !issuccess(B)
        throw(SingularException(B.info))
    end

    if B.rook
        if issymmetric(B)
            LAPACK.sytrs_rook!(B.uplo, B.LD, B.ipiv, R)
        else
            LAPACK.hetrs_rook!(B.uplo, B.LD, B.ipiv, R)
        end
    else
        if issymmetric(B)
            LAPACK.sytrs!(B.uplo, B.LD, B.ipiv, R)
        else
            LAPACK.hetrs!(B.uplo, B.LD, B.ipiv, R)
        end
    end
end
# There is no fallback solver for Bunch-Kaufman so we'll have to promote to same element type
function A_ldiv_B!(B::BunchKaufman{T}, R::StridedVecOrMat{S}) where {T,S}
    TS = promote_type(T,S)
    return A_ldiv_B!(convert(BunchKaufman{TS}, B), convert(AbstractArray{TS}, R))
end

function logabsdet(F::BunchKaufman)
    M = F.LD
    p = F.ipiv
    n = size(F.LD, 1)

    if !issuccess(F)
        return eltype(F)(-Inf), zero(eltype(F))
    end
    s = one(real(eltype(F)))
    i = 1
    abs_det = zero(real(eltype(F)))
    while i <= n
        if p[i] > 0
            elm = M[i,i]
            s *= sign(elm)
            abs_det += log(abs(elm))
            i += 1
        else
            # 2x2 pivot case. Make sure not to square before the subtraction by scaling
            # with the off-diagonal element. This is safe because the off diagonal is
            # always large for 2x2 pivots.
            if F.uplo == 'U'
                elm = M[i, i + 1]*(M[i,i]/M[i, i + 1]*M[i + 1, i + 1] -
                    (issymmetric(F) ? M[i, i + 1] : conj(M[i, i + 1])))
                s *= sign(elm)
                abs_det += log(abs(elm))
            else
                elm = M[i + 1,i]*(M[i, i]/M[i + 1, i]*M[i + 1, i + 1] -
                    (issymmetric(F) ? M[i + 1, i] : conj(M[i + 1, i])))
                s *= sign(elm)
                abs_det += log(abs(elm))
            end
            i += 2
        end
    end
    return abs_det, s
end

## reconstruct the original matrix
## TODO: understand the procedure described at
## http://www.nag.com/numeric/FL/nagdoc_fl22/pdf/F07/f07mdf.pdf

max_stars_repo_name: omus/julia

max_issues_repo_name: omus/julia

max_forks_repo_name: omus/julia


text:
module TailwindCSS

using Pkg.Artifacts
using Scratch
using JSON3

function __init__()
    global DIR = get_scratch!("TailwindCSS")
end

#-----------------------------------------------------------------------------# artifacts
artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
prefix = "https://github.com/tailwindlabs/tailwindcss/releases/download/v3.0.18"
tailwindcli_hash = artifact_hash("tailwindcli", artifacts_toml)

if tailwindcli_hash == nothing || !artifact_exists(tailwindcli_hash)
    tailwindcli_hash = create_artifact() do dir
        url = if Sys.isapple() && Sys.ARCH === :aarch64
            "$prefix/tailwindcss-macos-arm64"
        elseif Sys.isapple() && Sys.ARCH === :x86_64
            "$prefix/tailwindcss-macos-x64"
        elseif Sys.islinux() && Sys.ARCH === :aarch64
            "$prefix/tailwindcss-linux-arm64"
        elseif Sys.islinux() && Sys.ARCH === :x86_64
            "$prefix/tailwindcss-linux-x64"
        elseif Sys.iswindows()
            "$prefix/tailwindcss-windows-x64.exe"
        else
            error("TailwindCSS does not provide a CLI for your machine: $(Sys.MACHINE)")
        end
        download(url, joinpath(dir, "tailwindcli"))
        chmod(joinpath(dir, "tailwindcli"), 0o770, recursive=true)
    end
    bind_artifact!(artifacts_toml, "tailwindcli", tailwindcli_hash)
end
cli = joinpath(artifact_path(tailwindcli_hash), "tailwindcli")

#-----------------------------------------------------------------------------# function
help() = run(`$cli -h`)

#-----------------------------------------------------------------------------# init
"""
    init(file::String)
    init(file::String, config::AbstractDict)

Create a `tailwind.config.js` file, optionally using entries in `config`.

# Example

    using TailwindCSS: init

    init("tailwind.config.js", Dict("content" => ["*.html", "*.js"]))
"""
function init(file::String = "")
    args = [cli, "init"]
    !isempty(file) && push!(args, file)
    run(`$args`)
end
function init(file::String, config::AbstractDict)
    open(touch(file), "w") do io
        print(io, "module.exports = ")
        JSON3.write(io, config)
    end
end

#-----------------------------------------------------------------------------# minify
"""
    minify(inputfile, outputfile, config)

Write a minified CSS file `outputfile` based on `inputfile` and `config`.

- `config` can be a `String` (path to tailwind.config.js) or an `AbstractDict` (see `TailwindCSS.init`).
"""
function minify(input::String, output::String, config::String)
    args = [cli, "--input", input, "--output", output, "--config", config, "--minify"]
    run(`$args`)
end

function minify(input::String, output::String, config::AbstractDict)
    file = joinpath(DIR, "generated_tailwind.config.js")
    config = init(file, config)
    minify(input, output, file)
end

end


text:
function mdl_Zload!(voltage_r, voltage_i, current_r, current_i, device::PSY.PowerLoad)

    #Load squared voltage magnitude at steady state
    bus = PSY.get_bus(device)
    Vmag_sq = PSY.get_magnitude(bus)^2
    #Vmag_sq = 1.0

    #Load device parameters
    P = PSY.get_active_power(device)
    Q = PSY.get_reactive_power(device)

    #Compute impedance and load: Z = |V|^2/conj(S)
    #Z_load =  Vmag_sq/(P-Q*1im) #in pu
    #I = -(voltage_r[1] + voltage_i[1]*1im)/Z_load #in pu flowing out
    #current_r[1] += real(I)
    #current_i[1] += imag(I)

    #Update current
    #This model creates an equivalent RL/RC circuit based on steady state voltage
    current_r[1] += -(1.0 / Vmag_sq) * (voltage_r[1] * P + voltage_i[1] * Q) #in system pu flowing out
    current_i[1] += -(1.0 / Vmag_sq) * (voltage_i[1] * P - voltage_r[1] * Q) #in system pu flowing out

    return
end


text:
using PackageCompiler

create_sysimage(:PartitionedPoisson,
  sysimage_path=joinpath(@__DIR__,"..","PartitionedPoisson.so"),
  precompile_execution_file=joinpath(@__DIR__,"warmup.jl"))


text:
function sector_performance(datatype::String = "json")
    @argcheck in(datatype, ["json", "csv"])
    uri = "$(alphavantage_api)query?function=SECTOR&apikey=" * ENV["ALPHA_VANTAGE_API_KEY"]
    data = _get_request(uri)
    return _parse_response(data, datatype)
end
export sector_performance