text:
# Data analysis and plotting related to decoupling of transmission and distribution

using Printf
using DataFrames
using StatsPlots

function report_dist_candidates_pcc_power(
        dist_candidates::Dict{String,Any},
        out_dir::String;
        filename::String = "pcc_power",
        candidate_ids::Vector{String} = sort(collect(keys(dist_candidates))),
        report_deviation::Bool = true,
        report_flex_band::Bool = true,
        plot::Bool = true,
        plot_ext::String = "pdf",
        plot_kwargs...
    )
    mkpath(out_dir)
    periods = parse.(Int,first(values(dist_candidates))["ids"]["nw"])
    pcc_power = DataFrame()
    pcc_power.period = periods
    for cand_id in candidate_ids
        res = dist_candidates[cand_id]["result"]
        pcc_power[!,"$(cand_id)_import"] = [res["import"]["solution"]["nw"]["$t"]["td_coupling"]["p"] for t in periods]
        if haskey(res, "base")
            pcc_power[!,"$(cand_id)_base"] = [res["base"]["solution"]["nw"]["$t"]["td_coupling"]["p"] for t in periods]
        end
        pcc_power[!,"$(cand_id)_export"] = [res["export"]["solution"]["nw"]["$t"]["td_coupling"]["p"] for t in periods]
    end
    CSV.write(normpath(out_dir,"$filename.csv"), pcc_power)
    if plot
        plt = @df pcc_power Plots.plot(:period, cols(2:ncol(pcc_power));
            title  = "Power exchange at PCC",
            ylabel = "Imported power [MW]",
            xlabel = "Period",
            legend = :outertopright,
            plot_kwargs...
        )
        savefig(plt, normpath(out_dir,"$filename.$plot_ext"))
    end
    if report_deviation
        pcc_deviation = DataFrame()
        pcc_deviation.period = pcc_power.period
        for cand_id in candidate_ids
            res = dist_candidates[cand_id]["result"]
            if haskey(res, "base")
                pcc_deviation[!,"$(cand_id)_import"] = pcc_power[!,"$(cand_id)_import"] - pcc_power[!,"$(cand_id)_base"]
                pcc_deviation[!,"$(cand_id)_export"] = pcc_power[!,"$(cand_id)_export"] - pcc_power[!,"$(cand_id)_base"]
            end
        end
        CSV.write(normpath(out_dir,"$(filename)_deviation.csv"), pcc_deviation)
        if plot
            plt = @df pcc_deviation Plots.plot(:period, cols(2:ncol(pcc_deviation));
                title  = "Difference in power exchange at PCC w.r.t. base case",
                ylabel = "Power [MW]",
                xlabel = "Period",
                legend = :outertopright,
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_deviation.$plot_ext"))
        end
    end
    if report_flex_band
        pcc_flex_band = DataFrame()
        pcc_flex_band.period = pcc_power.period
        for cand_id in candidate_ids
            pcc_flex_band[!,"$(cand_id)"] = pcc_power[!,"$(cand_id)_import"] - pcc_power[!,"$(cand_id)_export"]
        end
        CSV.write(normpath(out_dir,"$(filename)_flex_band.csv"), pcc_flex_band)
        if plot
            plt = @df pcc_flex_band Plots.plot(:period, cols(2:ncol(pcc_flex_band));
                title  = "Flexibility band at PCC",
                ylabel = "Power [MW]",
                xlabel = "Period",
                legend = :outertopright,
                ylims  = (0.0, Inf),
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_flex_band.$plot_ext"))
        end
    end
end

function report_dist_candidates_branch(
        dist_candidates::Dict{String,Any},
        out_dir::String,
        mn_data::Dict{String,Any};
        filename::String = "branch",
        candidate_ids::Vector{String} = sort(collect(keys(dist_candidates))),
        plot::Bool = true,
        plot_ext::String = "pdf",
        plot_kwargs...
    )
    sn_data = mn_data["nw"]["$(first(_FP.nw_ids(mn_data; hour=1, scenario=1)))"]
    function _calc_p_rel(p, q, rated)
        p = abs(p)
        q = abs(q)
        if q <= sin(π/8) * rated
            return p / (cos(π/8)*rated)
        else
            return p / ((sin(π/8)+cos(π/8))*rated - q)
        end
    end
    mkpath(out_dir)
    branch = DataFrame(candidate = String[], result = String[], period = Int[], type = String[], id = Int[], p = Float64[], q = Float64[], p_rel = Float64[])
    for cand_id in candidate_ids
        for (res_id, res) in dist_candidates[cand_id]["result"]
            for (n, nw) in res["solution"]["nw"]
                for comp in ("branch", "ne_branch")
                    for (b, br) in get(nw, comp, Dict{String,Any}())
                        rated = sn_data[comp][b]["rate_a"]
                        p = br["pf"]
                        q = br["qf"]
                        p_rel = _calc_p_rel(p, q, rated)
                        push!(branch, (cand_id, res_id, parse(Int,n), comp, parse(Int, b), p, q, p_rel))
                    end
                end
            end
        end
    end
    sort!(branch, 1:5)
    CSV.write(normpath(out_dir,"$(filename).csv"), branch)

    if plot
        branch.cand_res = branch.candidate .* "_" .* branch.result
        branch.b = branch.type .* "_" .* string.(branch.id)
        branch_prel = DataFrames.unstack(branch, [:cand_res, :period], :b, :p_rel)
        for cand_res_branch in DataFrames.groupby(branch_prel, :cand_res)
            cand_res = cand_res_branch[1, :cand_res]
            plt = @df cand_res_branch Plots.plot(:period, cols(3:ncol(branch_prel));
                title  = "Relative active power of branches, $cand_res",
                xlabel = "Period",
                legend = :outertopright,
                legendfontsize = 5,
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_$(cand_res).$plot_ext"))
        end
    end
end

function report_dist_candidates_storage(
        dist_candidates::Dict{String,Any},
        out_dir::String;
        filename::String = "storage",
        candidate_ids::Vector{String} = sort(collect(keys(dist_candidates))),
        plot::Bool = true,
        plot_ext::String = "pdf",
        plot_kwargs...
    )
    mkpath(out_dir)
    storage_power = DataFrame(candidate = String[], result = String[], period = Int[], type = String[], id = Int[], p = Float64[])
    for cand_id in candidate_ids
        for (res_id, res) in dist_candidates[cand_id]["result"]
            for (n, nw) in res["solution"]["nw"]
                for (comp, p_key) in ("storage" => "ps", "ne_storage" => "ps_ne")
                    for (s, st) in get(nw, comp, Dict{String,Any}())
                        push!(storage_power, (cand_id, res_id, parse(Int,n), comp, parse(Int, s), st[p_key]))
                    end
                end
            end
        end
    end
    sort!(storage_power, 1:5)
    if nrow(storage_power) > 0
        CSV.write(normpath(out_dir,"$(filename)_power.csv"), storage_power)
        storage_power.cand_res = storage_power.candidate .* "_" .* storage_power.result
        storage_power.s = storage_power.type .* "_" .* string.(storage_power.id)
        storage_p = DataFrames.unstack(storage_power, [:s, :period], :cand_res, :p)
        for a_storage in DataFrames.groupby(storage_p, :s)
            s = a_storage[1, :s]
            plt = @df a_storage Plots.plot(:period, cols(3:ncol(storage_p));
                title  = "$s",
                ylabel = "Absorbed power [MW]",
                xlabel = "Period",
                legend = :outertopright,
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_power_$(s).$plot_ext"))
        end

        storage_power_total = combine(DataFrames.groupby(storage_p, :period), names(storage_p)[3:end].=>sum.=>names(storage_p)[3:end])
        CSV.write(normpath(out_dir,"$(filename)_power_total.csv"), storage_power_total)
        if plot
            plt = @df storage_power_total Plots.plot(:period, cols(2:ncol(storage_power_total));
                title  = "Total storage",
                ylabel = "Absorbed power [MW]",
                xlabel = "Period",
                legend = :outertopright,
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_power_total.$plot_ext"))
        end
    end

    storage_energy = DataFrame(candidate = String[], result = String[], period = Int[], type = String[], id = Int[], e = Float64[])
    for cand_id in candidate_ids
        for (res_id, res) in dist_candidates[cand_id]["result"]
            if res_id == "base"
                for (n, nw) in res["solution"]["nw"]
                    for (comp, e_key) in ("storage" => "se", "ne_storage" => "se_ne")
                        for (s, st) in get(nw, comp, Dict{String,Any}())
                            push!(storage_energy, (cand_id, res_id, parse(Int,n), comp, parse(Int, s), st[e_key]))
                        end
                    end
                end
            end
        end
    end
    sort!(storage_energy, 1:5)
    if nrow(storage_energy) > 0
        CSV.write(normpath(out_dir,"$(filename)_energy.csv"), storage_energy)
        storage_energy.cand_res = storage_energy.candidate .* "_" .* storage_energy.result
        storage_energy.s = storage_energy.type .* "_" .* string.(storage_energy.id)
        storage_e = DataFrames.unstack(storage_energy, [:s, :period], :cand_res, :e)
        for a_storage in DataFrames.groupby(storage_e, :s)
            s = a_storage[1, :s]
            plt = @df a_storage Plots.plot(:period, cols(3:ncol(storage_e));
                title  = "$s",
                ylabel = "Stored energy [MWh]",
                xlabel = "Period",
                legend = :outertopright,
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_energy_$(s).$plot_ext"))
        end

        storage_energy_total = combine(DataFrames.groupby(storage_e, :period), names(storage_e)[3:end].=>sum.=>names(storage_e)[3:end])
        CSV.write(normpath(out_dir,"$(filename)_energy_total.csv"), storage_energy_total)
        if plot
            plt = @df storage_energy_total Plots.plot(:period, cols(2:ncol(storage_energy_total));
                title  = "Total storage",
                ylabel = "Stored energy [MWh]",
                xlabel = "Period",
                legend = :outertopright,
                plot_kwargs...
            )
            savefig(plt, normpath(out_dir,"$(filename)_energy_total.$plot_ext"))
        end
    end
end

function report_dist_candidates_investment(
        dist_candidates::Dict{String,Any},
        out_dir::String;
        filename::String = "investment",
        candidate_ids::Vector{String} = sort(collect(keys(dist_candidates)))
    )
    mkpath(out_dir)
    investment = DataFrame(candidate=String[], component=String[], id=Int[], built=Bool[])
    for (cand_id, candidate) in dist_candidates
        for (comp_name, component) in candidate["investment"]
            for (comp_id, comp_activation) in component
                push!(investment, (cand_id, comp_name, parse(Int,comp_id), comp_activation))
            end
        end
    end
    investments = DataFrames.unstack(investment, :candidate, :built)
    select!(investments, ["component", "id", candidate_ids...]) # To change the order of columns
    CSV.write(normpath(out_dir,"$filename.csv"), investments)
end

function report_dist_candidates_cost(
        dist_candidates::Dict{String,Any},
        out_dir::String;
        filename::String = "cost",
        candidate_ids::Vector{String} = sort(collect(keys(dist_candidates)))
    )
    mkpath(out_dir)
    cost = DataFrame(candidate=String[], cost=Float64[])
    for cand_id in candidate_ids
        push!(cost, (cand_id, dist_candidates[cand_id]["cost"]))
    end
    CSV.write(normpath(out_dir,"$filename.csv"), cost)
end

function report_dist_candidates_nw_summary(
        dist_candidates::Dict{String,Any},
        out_dir::String;
        subdir::String = "nw_summary",
        candidate_ids::Vector{String} = sort(collect(keys(dist_candidates)))
    )
    for cand_id in candidate_ids
        for (res_id, res) in dist_candidates[cand_id]["result"]
            summary_dir = mkpath(normpath(out_dir, subdir, cand_id, res_id))
            for (n, nw) in res["solution"]["nw"]
                file = normpath(summary_dir,"$(@sprintf("%04i",parse(Int,n))).txt")
                open(file, "w") do io
                    _PM.summary(io, nw)
                end
            end
        end
    end
end


text:
#---------------------------------------------------------------------
#*	liquid and vapour phases
#*--------------------------------------------------------------------
type tank_liqvap
	tank_liqvap()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External physical properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of components",
				:Default=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0,
				:Symbol=>"_{in}"
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Intermediary liquid outlet stream",
				:Symbol=>"_{outmL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Outlet vapour stream",
				:Symbol=>"_{outV}"
			)),
			vol_tank(Dict{Symbol,Any}(
				:Brief=>"Routine to volume tank calculation",
				:Symbol=>"_{tank}"
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Component molar holdup"
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup"
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup"
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Internal energy"
			)),
			heat_rate(Dict{Symbol,Any}(
				:Brief=>"Reactor duty",
				:Default=>0
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume"
			)),
			[
				:(diff(M) = Inlet.F*Inlet.z - (OutletmL.F*OutletmL.z + OutletV.F*OutletV.z)),
				:(M = ML*OutletmL.z + MV*OutletV.z),
				:(sum(OutletmL.z) = 1),
				:(sum(OutletmL.z) = sum(OutletV.z)),
				:(OutletV.v = 1),
				:(OutletmL.v = 0),
				:(diff(E) = Inlet.F*Inlet.h - (OutletmL.F*OutletmL.h + OutletV.F*OutletV.h) + Q),
				:(E = ML*OutletmL.h + MV*OutletV.h),
				:(Tank.V = ML*vL + MV*PP.VapourVolume(OutletV.T,OutletV.P,OutletV.z)),
				:(PP.LiquidFugacityCoefficient(OutletmL.T,OutletmL.P,OutletmL.z)*OutletmL.z = PP.VapourFugacityCoefficient(OutletV.T,OutletV.P,OutletV.z)*OutletV.z),
				:(OutletmL.P = OutletV.P),
				:(OutletmL.T = OutletV.T),
				:(vL = PP.LiquidVolume(OutletmL.T,OutletmL.P,OutletmL.z)),
				:(ML*vL = Tank.V),
			],
			[
				"Component molar balance","Molar holdup","Mole fraction normalisation","Mole fraction normalisation","Vapourisation fraction","Vapourisation fraction","Energy balance","Total internal energy","Geometry constraint","Chemical Equilibrium","Mechanical Equilibrium","Thermal Equilibrium","Liquid Volume","Tank Level",
			],
			[:PP,:NComp,],
			[:Inlet,:OutletmL,:OutletV,:Tank,:M,:ML,:MV,:E,:Q,:vL,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Inlet::stream
	OutletmL::liquid_stream
	OutletV::vapour_stream
	Tank::vol_tank
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	Q::heat_rate
	vL::volume_mol
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export tank_liqvap
function setEquationFlow(in::tank_liqvap)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
end
function atributes(in::tank_liqvap,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Brief]="Model of a generic two-phase tank"
	drive!(fields,_)
	return fields
end
tank_liqvap(_::Dict{Symbol,Any})=begin
	newModel=tank_liqvap()
	newModel.attributes=atributes(newModel,_)
	newModel
end

ext: jl


text:
"""
    random_init(
        ncellx,ncelly,
        mean_ρ,range_ρ,
        mean_θ,range_θ;
        bcond_x="periodic",bcond_y="periodic"
    )

Sample a random data `(ρ,u,v)` of size `(ncellx+2)*(ncelly+2)`. The density `ρ` is sampled
uniformly around `mean_ρ` in an interval of size `range_ρ`. The velocity component `u` and
`v` are respectively given by `cos(θ)` and `sin(θ)` where `θ` is sampled uniformly around
`mean_θ` in an interval of size `range_θ`. Boundary conditions are then applied.
"""
function random_init(
    ncellx,ncelly,
    mean_ρ,range_ρ,
    mean_θ,range_θ;
    bcond_x="periodic",bcond_y="periodic"
)
    ρ = mean_ρ .+ range_ρ .* (rand(ncellx+2,ncelly+2) .- 0.5)
    θ = mean_θ .+ range_θ .* (rand(ncellx+2,ncelly+2) .- 0.5)
    u = cos.(θ)
    v = sin.(θ)
    boundary_conditions_ρuv!(ρ,u,v,bcond_x,bcond_y)
    return ρ,u,v
end

"""
    quadratic_potential_force(ncellx,ncelly,Δx,Δy,C)

Return the components `Fx` and `Fy` of the force exerted by the quadratic radial
potential V(r) = Cr²/2 centered around the midpoint of the box `[0,Lx]*[0,Ly]`
where `Lx = ncellx * Δx` and `Ly = ncelly * Δy`. For convenience, the outputs `Fx` and `Fy` are
matrices of size `(ncellx+2)*(ncelly+2)` with zeros on the boundary.
"""
function quadratic_potential_force(ncellx,ncelly,Δx,Δy,C)
    Fx = zeros(ncellx+2,ncelly+2)
    Fy = zeros(ncellx+2,ncelly+2)
    Lx = ncellx * Δx
    Ly = ncelly * Δy
    for i in 2:ncellx+1
        for j in 2:ncelly+1
            mx = Lx/2
            my = Ly/2
            xi = (i-2)*Δx + Δx/2
            yj = (j-2)*Δy + Δy/2
            Fx[i,j] = -C * (xi - mx)
            Fy[i,j] = -C * (yj - my)
        end
    end
    return Fx, Fy
end
"""
    flat_quadratic_potential_force(ncellx,ncelly,Δx,Δy,r0,C)

Return the components `Fx` and `Fy` of the force exerted by the quadratic radial
potential defined by V(r) = Cr²/2 for r>r0 and V(r)=0 otherwise,
centered around the midpoint of the box `[0,Lx]*[0,Ly]` where `Lx = ncellx * Δx`
and `Ly = ncelly * Δy`. For convenience, the outputs `Fx` and `Fy` are
matrices of size `(ncellx+2)*(ncelly+2)` with zeros on the boundary.
"""
function flat_quadratic_potential_force(ncellx,ncelly,Δx,Δy,r0,C)
    Fx = zeros(ncellx+2,ncelly+2)
    Fy = zeros(ncellx+2,ncelly+2)
    Lx = ncellx * Δx
    Ly = ncelly * Δy
    for i in 2:ncellx+1
        for j in 2:ncelly+1
            mx = Lx/2
            my = Ly/2
            xi = (i-2)*Δx + Δx/2
            yj = (j-2)*Δy + Δy/2
            rij = sqrt((xi-mx)^2+(yj-my)^2)
            if rij > r0
                Fx[i,j] = -C * (rij - r0) * (xi - mx)/rij
                Fy[i,j] = -C * (rij - r0) * (yj - my)/rij
            end
        end
    end
    return Fx, Fy
end


text:
module ProgressLogging

# Use README as the docstring of the module:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String) ProgressLogging

export @progress, @progressid, @withprogress, @logprogress

using Base.Meta: isexpr
using UUIDs: UUID, uuid4
using Logging: Logging, @logmsg, LogLevel

if VERSION >= v"1.1-"
    using UUIDs: uuid5
else
    import SHA
    function uuid5(ns::UUID, name::String)
        nsbytes = zeros(UInt8, 16)
        nsv = ns.value
        for idx in Base.OneTo(16)
            nsbytes[idx] = nsv >> 120
            nsv = nsv << 8
        end
        hash_result = SHA.sha1(append!(nsbytes, convert(Vector{UInt8}, codeunits(unescape_string(name)))))
        # set version number to 5
        hash_result[7] = (hash_result[7] & 0x0F) | (0x50)
        hash_result[9] = (hash_result[9] & 0x3F) | (0x80)
        v = zero(UInt128)
        #use only the first 16 bytes of the SHA1 hash
        for idx in Base.OneTo(16)
            v = (v << 0x08) | hash_result[idx]
        end
        return UUID(v)
    end
end

const ProgressLevel = LogLevel(-1)

"""
    ProgressLogging.ROOTID

This is used as `parentid` of root [`Progress`](@ref)es.
"""
const ROOTID = UUID(0)

"""
    ProgressLogging.Progress(id, [fraction]; [parentid, name, done])

# Usage: Progress log record provider

Progress log record can be created by using the following pattern

```julia
id = uuid4()
try
    @info Progress(id)  # create a progress bar
    # some time consuming job
    # ...
    @info Progress(id, 0.1)  # update progress to 10%
    # ...
finally
    @info Progress(id, done = true)  # close the progress bar
end
```

It is recommended to use [`@withprogress`](@ref),
[`@logprogress`](@ref), and optionally [`@progressid`](@ref) to create
log records.

# Usage: Progress log record consumer (aka progress monitor)

It is recommended to use [`ProgressLogging.asprogress`](@ref) instead
of checking `message isa Progress`.  Progress monitors can retrieve
progress-related information from the following properties.

# Properties
- `fraction::Union{Float64,Nothing}`: it can take following values:
  - `0 <= fraction < 1`
  - `fraction >= 1`: completed
  - `fraction = nothing`: indeterminate progress
- `id::UUID`: Identifier of the job whose progress is at `fraction`.
- `parentid::UUID`: The ID of the parent progress.  It is set to
  [`ProgressLogging.ROOTID`](@ref) when there is no parent progress.
  This is used for representing progresses of nested jobs.  Note that
  sub-jobs may be executed concurrently; i.e., there can be multiple
  child jobs for one parent job.
- `name::String`: Name of the progress bar.
- `done::Bool`: `true` if the job is done.
"""
struct Progress
    id::UUID
    parentid::UUID
    fraction::Union{Float64,Nothing}
    name::String
    done::Bool

    function Progress(id, parentid, fraction, name, done)
        if fraction isa Real && isnan(fraction)
            fraction = nothing
        end
        return new(id, parentid, fraction, name, done)
    end
end

Progress(;
    id::UUID,
    parentid::UUID = ROOTID,  # not nested by default
    fraction::Union{Real,Nothing} = nothing,
    name::String = "",
    done::Bool = false,
) = Progress(id, parentid, fraction, name, done)

Progress(id::UUID, fraction::Union{Real,Nothing} = nothing; kwargs...) =
    Progress(; kwargs..., fraction = fraction, id = id)

# Define `string`/`print` so that progress log records are (somewhat)
# readable even without specific log monitors.
function Base.print(io::IO, progress::Progress)
    print(io, isempty(progress.name) ? "Progress" : progress.name)
    if progress.parentid !== ROOTID
        print(io, " (sub)")
    end
    print(io, ": ")
    if progress.fraction === nothing
        print(io, "??%")
    else
        print(io, floor(Int, progress.fraction * 100), '%')
    end
    return
end

const PROGRESS_LOGGING_UUID_NS = UUID("1e962757-ea70-431a-b9f6-aadf988dcb7f")

asuuid(id::UUID) = id
asuuid(id) = uuid5(PROGRESS_LOGGING_UUID_NS, repr(id))

"""
    ProgressLogging.asprogress(_, name, _, _, id, _, _; progress, ...) :: Union{Progress, Nothing}

Pre-process log record to obtain a [`Progress`](@ref) object if it is
one of the supported format.  This is mean to be used with the
`message` positional argument and _all_ keyword arguments passed to
`Logging.handle_message`.  Example:

```julia
function Logging.handle_message(logger::MyLogger, args...; kwargs...)
    progress = ProgressLogging.asprogress(args...; kwargs...)
    if progress !== nothing
        return # handle progress log record
    end
    # handle normal log record
end
```
"""
asprogress(_level, progress::Progress, _args...; _...) = progress
function asprogress(
    _level,
    name,
    _module,
    _group,
    id,
    _file,
    _line;
    progress = undef,  # `undef` is an arbitrary unsupported value
    kwargs...,
)
    if progress isa Union{Nothing,Real,AbstractString}
        return _asprogress(name, id; progress = progress, kwargs...)
    else
        return nothing
    end
end

# `parentid` is used from `@logprogress`.
function _asprogress(name, id, parentid = ROOTID; progress, _...)
    if progress isa Union{Nothing,Real}
        fraction = progress
    elseif progress == "done"
        fraction = nothing
    else
        return nothing
    end
    return Progress(
        fraction = fraction,
        name = name,
        id = asuuid(id),
        parentid = parentid,
        done = progress == "done",
    )
end

# To pass `Progress` value without breaking progress monitors with the
# previous `progress` key based specification, we create a custom
# string type that has `Progress` attached to it.  This is used as the
# third argument `message` of `Logging.handle_message`.
struct ProgressString <: AbstractString
    progress::Progress
end

asprogress(_level, str::ProgressString, _args...; _...) = str.progress

# Since `Base.string(::AbstractString)` is defined to be an `identity`
# function, we overload it to make sure that `string(message)`
# typically used in the loggers converts `ProgressString` to a vanilla
# `String` as soon as possible.  It may not be needed if we define
# `ProgressString` perfectly.  But let's play on the safe side.
Base.string(str::ProgressString) = str.progress.name

Base.print(io::IO, str::ProgressString) = print(io, string(str))
Base.convert(::Type{ProgressString}, str::ProgressString) = str
Base.convert(::Type{T}, str::ProgressString) where {T<:AbstractString} =
    convert(T, str.progress.name)

# Define `cmp` to make `==` etc. work
Base.cmp(a::AbstractString, b::ProgressString) = cmp(a, string(b))
Base.cmp(a::ProgressString, b::AbstractString) = cmp(string(a), b)
Base.cmp(a::ProgressString, b::ProgressString) = cmp(string(a), string(b))

# Avoid using `show(::IO, ::AbstractString)` which expects
# `Base.print_quoted` to work.
function Base.show(io::IO, str::ProgressString)
    if get(io, :typeinfo, Any) === ProgressString
        show(io, string(str))
        return
    end
    print(io, @__MODULE__, ".")
    print(io, "ProgressString(")
    show(io, str.progress)
    print(io, ")")
end

"""
    progress(f::Function; name = "")

Evaluates `f` with `id` as its argument and makes sure to destroy the progress
bar afterwards. To update the progress bar in `f` you can call a logging statement
like `@info` or even just `@logmsg` with `_id=id` and `progress` as arguments.

`progress` can take either of the following values:
  - `0 <= progress < 1`: create or update progress bar
  - `progress == nothing || progress = NaN`: set progress bar to indeterminate progress
  - `progress >= 1 || progress == "done"`: destroy progress bar

The logging message (e.g. `"foo"` in `@info "foo"`) will be used as the progress
bar's name.

Log level must be higher or equal to `$ProgressLevel`.

```julia
ProgressLogging.progress() do id
    for i = 1:10
        sleep(0.5)
        @info "iterating" progress=i/10 _id=id
    end
end
```
"""
function progress(f; name = "")
    _id = uuid4()
    @logmsg ProgressLevel name progress = nothing _id = _id
    try
        f(_id)
    finally
        @logmsg ProgressLevel name progress = "done" _id = _id
    end
end

const _id_var = gensym(:progress_id)
const _parentid_var = gensym(:progress_parentid)
const _name_var = gensym(:progress_name)

"""
    @withprogress [name=""] [parentid=uuid4()] ex

Create a lexical environment in which [`@logprogress`](@ref) can be used to
emit progress log events without manually specifying the log level, `_id`,
and name (log message).

```julia
@withprogress name="iterating" begin
    for i = 1:10
        sleep(0.5)
        @logprogress i/10
    end
end
```
"""
macro withprogress(exprs...)
    _withprogress(exprs...)
end

function _withprogress(exprs...)
    length(exprs) == 0 &&
        throw(ArgumentError("`@withprogress` requires at least one expression."))

    m = @__MODULE__

    kwargs = Dict(:name => "", :parentid => :($m.@progressid()))
    unsupported = []
    for kw in exprs[1:end-1]
        if isexpr(kw, :(=)) && length(kw.args) == 2 && haskey(kwargs, kw.args[1])
            kwargs[kw.args[1]] = kw.args[2]
        else
            push!(unsupported, kw)
        end
    end

    # Error on invalid input expressions:
    if !isempty(unsupported)
        msg = sprint() do io
            println(io, "Unsupported optional arguments:")
            for kw in unsupported
                println(io, kw)
            end
            print(io, "`@withprogress` supports only following keyword arguments: ")
            join(io, keys(kwargs), ", ")
        end
        throw(ArgumentError(msg))
    end

    ex = exprs[end]
    quote
        let $_parentid_var = $(kwargs[:parentid]),
            $_id_var = $uuid4(),
            $_name_var = $(kwargs[:name])
            $m.@logprogress nothing
            try
                $ex
            finally
                $m.@logprogress "done"
            end
        end
    end |> esc
end

"""
    @logprogress [name] progress [key1=val1 [key2=val2 ...]]

This macro must be used inside [`@withprogress`](@ref) macro.

Log a progress event with a value `progress`.  The expression
`progress` must be evaluated to be a real number between `0` and `1`
(inclusive), a `NaN`, or a string `"done"`.

Optional first argument `name` can be used to change the name of the
progress bar.  Additional keyword arguments are passed to `@logmsg`.
"""
macro logprogress(name, progress = nothing, args...)
    name_expr = :($Base.@isdefined($_name_var) ? $_name_var : "")
    if progress == nothing
        # Handle: @logprogress progress
        kwargs = (:(progress = $name), args...)
        progress = name
        name = name_expr
    elseif isexpr(progress, :(=)) && progress.args[1] isa Symbol
        # Handle: @logprogress progress key1=val1 ...
        kwargs = (:(progress = $name), progress, args...)
        progress = name
        name = name_expr
    else
        # Otherwise, it's: @logprogress name progress key1=val1 ...
        kwargs = (:(progress = $progress), args...)
    end

    id_err = "`@logprogress` must be used inside `@withprogress` or with `_id` keyword argument"
    id_expr = :($Base.@isdefined($_id_var) ? $_id_var : $error($id_err))
    for x in kwargs
        if isexpr(x, :(=)) && x.args[1] === :_id
            id_expr = :($asuuid($(x.args[2])))
            # last set wins; so not `break`ing
        end
    end

    @gensym id_tmp
    # Emitting progress log record as old/open API (i.e., using
    # `progress` key) and _also_ as new API based on `Progress` type.
    msgexpr = :($ProgressString($_asprogress(
        $name,
        $id_tmp,
        $_parentid_var;
        progress = $progress,
    )))
    quote
        $id_tmp = $id_expr
        $Logging.@logmsg($ProgressLevel, $msgexpr, $(kwargs...), _id = $id_tmp)
    end |> esc
end

"""
    @progressid

Get the progress ID of current lexical scope.
"""
macro progressid()
    quote
        $Base.@isdefined($_id_var) ? $_id_var : $ROOTID
    end |> esc
end

"""
    @progress [name="", threshold=0.005] for i = ..., j = ..., ...
    @progress [name="", threshold=0.005] x = [... for i = ..., j = ..., ...]

Show a progress meter named `name` for the given loop or array comprehension
if possible. Update frequency is limited by `threshold` (one update per 0.5% of
progress by default).
"""
macro progress(args...)
    esc(_progress(args...))
end

_progress(ex) = _progress("", 0.005, ex)
_progress(name::Union{AbstractString,Expr}, ex) = _progress(name, 0.005, ex)
_progress(thresh::Real, ex) = _progress("", thresh, ex)

function _progress(name, thresh, ex)
    if ex.head == Symbol("=") &&
       ex.args[2].head == :comprehension && ex.args[2].args[1].head == :generator
        # comprehension: <target> = [<body> for <iter_var> in <range>,...]
        loop = _comprehension
        target = ex.args[1]
        result = target
        gen_ex = ex.args[2].args[1]
        body = gen_ex.args[1]
        iter_exprs = gen_ex.args[2:end]
        iter_vars = [e.args[1] for e in iter_exprs]
        ranges = [e.args[2] for e in iter_exprs]
    elseif ex.head == :for && ex.args[1].head == Symbol("=") && ex.args[2].head == :block
        # single-variable for: for <iter_var> = <range>; <body> end
        loop = _for
        target = :_
        result = :nothing
        iter_vars = [ex.args[1].args[1]]
        ranges = [ex.args[1].args[2]]
        body = ex.args[2]
    elseif ex.head == :for && ex.args[1].head == :block && ex.args[2].head == :block
        # multi-variable for: for <iter_var> = <range>,...; <body> end
        loop = _for
        target = :_
        result = :nothing
        # iter_vars and ranges are ordered from inner loop to outer loop, for
        # consistent computation of progress between for loops and comprehensions
        iter_vars = reverse([e.args[1] for e in ex.args[1].args])
        ranges = reverse([e.args[2] for e in ex.args[1].args])
        body = ex.args[2]
    else
        error("@progress requires a for loop (for i in irange, j in jrange, ...; <body> end) " *
              "or array comprehension with assignment (x = [<body> for i in irange, j in jrange, ...])")
    end
    _progress(name, thresh, ex, target, result, loop, iter_vars, ranges, body)
end

function _progress(name, thresh, ex, target, result, loop, iter_vars, ranges, body)
    count_vars = [gensym(Symbol("i$k")) for k = 1:length(iter_vars)]
    iter_exprs = [:(($i, $v) = $enumerate($r)) for (i, v, r) in zip(
        count_vars,
        iter_vars,
        ranges,
    )]
    @gensym count_to_frac val frac lastfrac
    m = @__MODULE__
    quote
        $target = $m.@withprogress name = $name begin
            $count_to_frac = $make_count_to_frac($(ranges...))
            $lastfrac = 0.0

            $(loop(
                iter_exprs,
                quote
                    $val = $body
                    $frac = $count_to_frac($(count_vars...))
                    if $frac - $lastfrac > $thresh
                        $m.@logprogress $frac
                        $lastfrac = $frac
                    end
                    $val
                end,
            ))
        end
        $result
    end
end

_comprehension(iter_exprs, body) =
    Expr(:comprehension, Expr(:generator, body, iter_exprs...))
_for(iter_exprs, body) = Expr(:for, Expr(:block, reverse(iter_exprs)...), body)

taccumulate(op, ::Tuple{}) = ()
function taccumulate(op::F, xs::Tuple) where {F}
    ys, = foldl(Base.tail(xs); init=((xs[1],), xs[1])) do (ys, acc), x
        acc = op(acc, x)
        (ys..., acc), acc
    end
    return ys
end

function make_count_to_frac(iterators...)
    lens = map(length, iterators)
    n = prod(lens)
    strides = (1, taccumulate(*, Base.front(lens))...)
    function count_to_frac(idxs...)
        offsets = map(i -> i - 1, idxs)
        total = sum(map(*, offsets, strides)) + 1
        return total / n
    end
    return count_to_frac
end

end # module


text:
(*)(op::ProjectionSampling, object::SynthesisOperator) = apply(op, object)

apply(op::SamplingOperator, dict::Dictionary; options...) =
	apply(op, SynthesisOperator(dict); options...)

(*)(op::SamplingOperator, object::SynthesisOperator) = apply(op, object)

function apply(sampling::SamplingOperator, op::SynthesisOperator; options...)
	T = operatoreltype(dictionary(op), dest(sampling))
    evaluation(T, dictionary(op), grid(sampling); options...)
end

function apply(analysis::ProjectionSampling, synthesis::SynthesisOperator; options...)
	T = operatoreltype(dictionary(synthesis), dest(analysis))
	synthesisdict = dictionary(synthesis)
	analysismeasure = measure(analysis)
	analysisdict = dictionary(analysis)
	# Check whether we can just compute the gram matrix instead.
	if isequal(analysisdict, synthesisdict) && size(analysisdict)==size(synthesisdict)
		gram(T, synthesisdict, analysismeasure; options...)
	else
		mixedgram(T, analysisdict, synthesisdict, analysismeasure; options...)
	end
end

function adjoint(op::SynthesisOperator)
    if op.measure == nothing
        @warn("Adjoint of $(typeof(op)) does not exist if measure is not defined. Returning inverse.")
        inv(op)
    else
        AnalysisOperator(dictionary(op), op.measure, dest_space(op))
    end
end

function adjoint(op::AnalysisOperator)
    if src_space(op) isa Span && iscompatible(dictionary(src_space(op)), dictionary(op))
        SynthesisOperator(dictionary(op), measure(op))
    else
        error("Adjoint of $(typeof(op)) does not exist if spaces do not match.")
    end
end


text:
export condpot
function condpot{T<:Potential}(p::T,x::IntOrIntArray=[],y::IntOrIntArray=[])
    # compute p(x|y)=p(x,y)/p(y)
    if isempty(x) &&!isempty(y)
        error("first argument of condpot is empty and second argument is non-empty")
    end
    if isempty(x) && isempty(y)
        return p/sum(p)
    else
        if isempty(y)
            return(sum(p,x,SumOver=false)/sum(p))
        else
            pxy=sum(p,[x y],SumOver=false)
            py=sum(pxy,y,SumOver=false)
            return pxy/py
        end
    end
end

import Base.sum
function sum{T<:Potential,I<:Integer}(pots::Dict{I,T})
    # Sum potentials into a single potential
    return sum(convert(Array{T},pots))
end
export sum # is this needed?

import Base.prod
function prod{I<:Integer,P<:Potential}(D::Dict{I,P})
    p=convert(Array{Potential},D)
    return prod(p)
end

function prod(P::Potential)
    return P
end

# converters :
import Base.convert
function convert{T<:Potential,DT<:Potential,I<:Integer}(::Type{Array{T}}, D::Dict{I,DT})
    L=length(collect(keys(D)))
    pot=Array(Potential,L)
    ky=collect(keys(D))
    for k=1:length(keys(D))
        pot[k]=D[ky[k]]
    end
    return pot
end
export convert

export multpots
function multpots{P<:Potential}(pot::Array{P})
    if length(pot)==1
        return pot[1]
        else
            ind=[]
            for i=1:length(pot)
                if isdefined(pot,i) # remove #undef entries (missing)
                    ind=vcat(ind, i)
                end
            end
        if length(ind)==1
            return pot[ind]
        else
            return prod(pot[ind])
        end
    end
end

export standardisevariables
#function standardisevariables{DictOrArray{T<:Potential}}(inpot::T)
    #function standardisevariables{T<:Potential}(inpot::DictOrArray{T})
    # FIXME! NEED TO FIGURE OUT HOW TO DISPATCH ON ONLY DICTORARRY(subtype of Potential)
    function standardisevariables(inpot)
    # This just relabels the first from starting from 1
    pot=DictToArray(deepcopy(inpot))
    variables=potvariables(pot)
    for i=1:length(pot)
        pot[i].variables=memberinds(pot[i].variables,variables)
    end
    return pot, variables
end

export returnvariables
#function returnvariables{DictOrArray{T<:Potential}}(inpot::T,originalvariables)
#function returnvariables{T<:Potential}(inpot::DictOrArray{T},originalvariables)
function returnvariables(inpot,originalvariables)
    pot=DictToArray(deepcopy(inpot))
    variables=potvariables(pot)
    for i=1:length(pot)
        pot[i].variables=originalvariables[memberinds(pot[i].variables,variables)]
    end
    return pot
end

export potvariables
#function [variables nstates con convec]=potvariables(inpot)
function potvariables(inpot)
    # NEED TO EXTEND THIS TO RETURN ADDITIONAL INFO, LIKE THE MATLAB FUNCTION
    #POTVARIABLES Returns information about all variables in a set of potentials
    # [variables nstates con convec]=potvariables(pot)
    #
    # return the variables and their number of states
    # If there is a dimension mismatch in the tables then return con=0
    # convec(i)=0 reports that variable i has conflicting dimension

    pot=deepcopy(inpot) # don't want to change anything in caller
    pot=DictToArray(pot)

    if isempty(pot);return;end

    if !isa(pot,Array)
        tmp=deepcopy(pot)
        pot=Array(Potential,1)
        pot[1]=tmp
    end

    variables=[];
    for p=1:length(pot)
        variables=union(variables,pot[p].variables)
    end

    return variables

end

export dag
function dag(pot)
    #DAG Return the Dictionary list of edge for a Belief Network
    # Assumes that pot[i] contains the distribution p(i|pa(i))

    L=Dict{Integer,Array{Integer}}()
    vars=[];
    for p=1:length(pot)
        vars=vcat(vars, pot[p].variables)
    end
    N=maximum(vars)
    A=zeros(N,N);
    #for p=1:length(pot)
    #    L[p]=setdiff(pot[p].variables,p)
    #end
    for p=1:length(pot) ## might be better to insist that pot is itself a Dict
        pa=setdiff(pot[p].variables,p)
        for i=1:length(pa)
            if haskey(L,pa[i])
                L[pa[i]]=union(L[pa[i]],p)
            else
                L[pa[i]]=union([],p)
            end
        end
    end

    return L

end

export whichpot
function whichpot(pot,variables,n=1)
    #%WHICHPOT Returns potential that contain a set of variables
#% potnum = whichpot(pot,variables,<n>)
    #%
    #% Return potential numbers that contain all the specified variables
    #% The cliques are returned with those containing the smallest number of
    #% variables first.
    #% If optional n is used, returns at most n potential number
    potnum=[]; nvars=[];
    for p=1:length(pot)
        # find the potential that contains variables
        if prod(memberinds(variables,pot[p].variables))>0
            nvars=vcat(nvars,length(pot[p].variables));
            potnum=vcat(potnum,p);
    end
    end
    val=sort(nvars)
    ind=memberinds(val,nvars)
        potnum=potnum[ind[1:n]]
    if n==1
        potnum=potnum[1]
    end

    return potnum
end

export numvariables
function numvariables{T<:Potential}(A::T)
    length(A.variables)
end


text:
"""
    contourf(xs, ys, zs; kwargs...)

Plots a filled contour of the height information in `zs` at horizontal grid positions `xs`
and vertical grid positions `ys`.

The attribute `levels` can be either
- an `Int` that produces n equally wide levels or bands
- an `AbstractVector{<:Real}` that lists n consecutive edges from low to high, which result in n-1 levels or bands

You can also set the `mode` attribute to `:relative`.
In this mode you specify edges by the fraction between minimum and maximum value of `zs`.
This can be used for example to draw bands for the upper 90% while excluding the lower 10% with `levels = 0.1:0.1:1.0, mode = :relative`.

In :normal mode, if you want to show a band from `-Inf` to the low edge, set `extendlow` to `:auto` for the same color as the first level, or specify a different color (default `nothing` means no extended band)
If you want to show a band from the high edge to `Inf`, set `extendhigh` to `:auto` for the same color as the last level, or specify a different color (default `nothing` means no extended band)

## Attributes
$(ATTRIBUTES)
"""
@recipe(Contourf) do scene
    Theme(
        levels = 10,
        mode = :normal,
        colormap = theme(scene, :colormap),
        extendlow = nothing,
        extendhigh = nothing,
        inspectable = theme(scene, :inspectable),
        transparency = false
    )
end

# these attributes are computed dynamically and needed for colorbar e.g.
# _computed_levels
# _computed_colormap
# _computed_extendlow
# _computed_extendhigh

function _get_isoband_levels(levels::Int, mi, ma)
    edges = Float32.(LinRange(mi, ma, levels+1))
end

function _get_isoband_levels(levels::AbstractVector{<:Real}, mi, ma)
    edges = Float32.(levels)
    @assert issorted(edges)
    edges
end

conversion_trait(::Type{<:Contourf}) = ContinuousSurface()

function _get_isoband_levels(::Val{:normal}, levels, values)
    _get_isoband_levels(levels, extrema_nan(values)...)
end

function _get_isoband_levels(::Val{:relative}, levels::AbstractVector, values)
    mi, ma = extrema_nan(values)
    Float32.(levels .* (ma - mi) .+ mi)
end

function Makie.plot!(c::Contourf{<:Tuple{<:AbstractVector{<:Real}, <:AbstractVector{<:Real}, <:AbstractMatrix{<:Real}}})
    xs, ys, zs = c[1:3]

    c.attributes[:_computed_levels] = lift(zs, c.levels, c.mode) do zs, levels, mode
        _get_isoband_levels(Val(mode), levels, vec(zs))
    end

    colorrange = lift(c._computed_levels) do levels
        minimum(levels), maximum(levels)
    end

    c.attributes[:_computed_colormap] = lift(c._computed_levels, c.colormap) do levels, cmap
        levels_scaled = (levels .- minimum(levels)) ./ (maximum(levels) - minimum(levels))
        cgrad(cmap, levels_scaled, categorical = true)
    end

    lowcolor = lift(Union{Nothing, RGBAf}, c.extendlow) do el
        if el === nothing
            nothing
        elseif el === automatic || el == :auto
            RGBAf(get(c._computed_colormap[], 0))
        else
            convert_attribute(el, key"color"())::RGBAf
        end
    end
    c.attributes[:_computed_extendlow] = lowcolor
    is_extended_low = lift(x -> !isnothing(x), lowcolor)

    highcolor = lift(Union{Nothing, RGBAf}, c.extendhigh) do eh
        if eh === nothing
            nothing
        elseif eh === automatic || eh == :auto
            RGBAf(get(c._computed_colormap[], 1))
        else
            convert_attribute(eh, key"color"())::RGBAf
        end
    end
    c.attributes[:_computed_extendhigh] = highcolor
    is_extended_high = lift(x -> !isnothing(x), highcolor)

    PolyType = typeof(Polygon(Point2f[], [Point2f[]]))

    polys = Observable(PolyType[])
    colors = Observable(Float64[])

    function calculate_polys(xs, ys, zs, levels::Vector{Float32}, is_extended_low, is_extended_high)
        empty!(polys[])
        empty!(colors[])

        levels = copy(levels)
        @assert issorted(levels)
        is_extended_low && pushfirst!(levels, -Inf)
        is_extended_high && push!(levels, Inf)
        lows = levels[1:end-1]
        highs = levels[2:end]

        nbands = length(lows)

        # zs needs to be transposed to match rest of makie
        isos = Isoband.isobands(xs, ys, zs', lows, highs)

        allvertices = Point2f[]
        allfaces = NgonFace{3,OffsetInteger{-1,UInt32}}[]
        allids = Int[]
        levelcenters = (highs .+ lows) ./ 2

        for (i, (center, group)) in enumerate(zip(levelcenters, isos))
            points = Point2f.(group.x, group.y)
            polygroups = _group_polys(points, group.id)
            for polygroup in polygroups
                outline = polygroup[1]
                holes = polygroup[2:end]
                push!(polys[], GeometryBasics.Polygon(outline, holes))
                # use contour level center value as color
                center_scaled = (center - colorrange[][1]) / (colorrange[][2] - colorrange[][1])
                push!(colors[], center)
            end
        end
        polys[] = polys[]
        return
    end

    onany(calculate_polys, xs, ys, zs, c._computed_levels, is_extended_low, is_extended_high)
    # onany doesn't get called without a push, so we call
    # it on a first run!
    calculate_polys(xs[], ys[], zs[], c._computed_levels[], is_extended_low[], is_extended_high[])

    poly!(c,
        polys,
        colormap = c._computed_colormap,
        colorrange = colorrange,
        color = colors,
        strokewidth = 0,
        strokecolor = :transparent,
        shading = false,
        inspectable = c.inspectable,
        transparency = c.transparency
    )
end

"""
    _group_polys(points, ids)

Given a vector of polygon vertices, and one vector of group indices, which
are assumed to be returned from the isoband algorithm, return
a vector of groups, where each group has one outer polygon and zero or more
inner polygons which are holes in the outer polygon. It is possible that one
group has multiple outer polygons with multiple holes each.
"""
function _group_polys(points, ids)

    polys = [points[ids .== i] for i in unique(ids)]
    npolys = length(polys)

    polys_lastdouble = [push!(p, first(p)) for p in polys]

    # this matrix stores whether poly i is contained in j
    # because the marching squares algorithm won't give us any
    # intersecting or overlapping polys, it should be enough to
    # check if a single point is contained, saving some computation time
    containment_matrix = [
        p1 != p2 &&
        PolygonOps.inpolygon(first(p1), p2) == 1
        for p1 in polys_lastdouble, p2 in polys_lastdouble]

    unclassified_polyindices = collect(1:size(containment_matrix, 1))
    # @show unclassified_polyindices

    # each group has first an outer polygon, and then its holes
    # TODO: don't specifically type this 2f0?
    groups = Vector{Vector{Point2f}}[]

    # a dict that maps index in `polys` to index in `groups` for outer polys
    outerindex_groupdict = Dict{Int, Int}()

    # all polys have to be classified
    while !isempty(unclassified_polyindices)
        to_keep = ones(Bool, length(unclassified_polyindices))

        # go over unclassifieds and find outer polygons in the remaining containment matrix
        for (ii, i) in enumerate(unclassified_polyindices)
            # an outer polygon is not inside any other polygon of the matrix
            if sum(containment_matrix[ii, :]) == 0
                # an outer polygon
                # println(i, " is an outer polygon")
                push!(groups, [polys_lastdouble[i]])
                outerindex_groupdict[i] = length(groups)
                # delete this poly from further rounds
                to_keep[ii] = false
            end
        end

        # go over unclassifieds and find hole polygons
        for (ii, i) in enumerate(unclassified_polyindices)
            # the hole polygons can only be in one polygon from the current group
            # if they are in more than one, they are "inner outer" or inner hole polys
            # and will be handled in one of the following passes
            if sum(containment_matrix[ii, :]) == 1
                outerpolyindex_of_unclassified = findfirst(containment_matrix[ii, :])
                outerpolyindex = unclassified_polyindices[outerpolyindex_of_unclassified]
                # a hole
                # println(i, " is an inner polygon of ", outerpolyindex)
                groupindex = outerindex_groupdict[outerpolyindex]
                push!(groups[groupindex], polys_lastdouble[i])
                # delete this poly from further rounds
                to_keep[ii] = false
            end
        end

        unclassified_polyindices = unclassified_polyindices[to_keep]
        containment_matrix = containment_matrix[to_keep, to_keep]
    end
    groups
end


text:
using PyPlot, LinearAlgebra, DynamicalSystems, OrdinaryDiffEq

function towel_map(u)
    x, y, z = u
    return [3.8*x*(1.0 - x) - 0.05*(y + 0.35)*(1.0 - 2.0*z),
            0.1*((y + 0.35)*(1.0 - 2.0*z) - 1.0)*(1.0 - 1.9*x),
            3.78*(1.0 - z) + 0.2*y]
end

#num_steps = 1000
#steps = Array{Float64}(undef, 3, num_steps)
#steps[:, 1] = [0.085, -0.121, 0.075]
#for i in 2:num_steps
#    steps[:, i] = towel_map(steps[:, i - 1])
#end
#plot3D(steps[1, :], steps[2, :], steps[3, :])

sys = Systems.towel()

function plot_map()
    steps = trajectory(sys, 10000)
    plot3D(steps[:, 1], steps[:, 2], steps[:, 3], ".", ms=0.5)
    #plot(steps[:, 1], steps[:, 2], ".", ms=0.5)
end

function dysys_lya()
#lya_spec = lyapunovspectrum(sys, 20)
    conv = [lyapunovspectrum(sys, i) for i in 1:40]
    for i in 1:3
        plot(1:size(conv, 1), [conv[j][i] for j in 1:size(conv, 1)])
    end
    xlabel("N")
    ylabel("\$\\lambda_1\$")
end

function computer_lyapunov_spectrum_with_qr_discrete(dds, x0, T, N)
    n = size(x0, 1)
    # this is super inefficient
    function unpack(y)
        x = y[1:n]
        Y = reshape(y[n+1:end], (n, n))
        return x, Y
    end
    function pack(x, Y)
        return SVector{n^2 + n}([x; reshape(Y, n^2)]...)
    end
    function g(y, p, t)
        x, Y = unpack(y)
        next_x = dds.f(x, p, t)
        next_Y = dds.jacobian(x, p, t) * Y
        return pack(next_x, next_Y)
    end
    Y0 = Matrix(1.0I, n, n)
    y0 = pack(x0, Y0)
    Y_dds = DiscreteDynamicalSystem(g, y0, dds.p)
    lyapunov_spectrum = zeros(n)
    for i in 1:N
        y = trajectory(Y_dds, T, y0)[end]
        x, Y = unpack(y)
        Q, R = qr(Y) # the problem is that julia doesn't return the decompostion, where R_jj > 0 forall j

        Y0 = sign.(diag(R)) .* Q
        x0 = x
        y0 = pack(x0, Y0)
        lyapunov_spectrum .+= log.(abs.(diag(R)))
    end
    return lyapunov_spectrum ./ (T*N)
end

computer_lyapunov_spectrum_with_qr_discrete(sys, sys.u0, 3, 40)

max_stars_repo_name: cosmo-jana/numerics-physics-stuff

max_issues_repo_name: cosmo-jana/numerics-physics-stuff

max_forks_repo_name: cosmo-jana/numerics-physics-stuff


text:
using Tweaks.BiMaps

@testset "init" begin
    map = BiMap(Dict('a'=>1, 'b'=>2), Dict(1=>'a', 2=>'b'))
    @test Dict(map) == Dict('a'=>1, 'b'=>2)
    @test Dict(rev(map)) == Dict(1=>'a', 2=>'b')

    map = BiMap(Dict('a'=>1, 'b'=>2))
    @test Dict(map) == Dict('a'=>1, 'b'=>2)
    @test Dict(rev(map)) == Dict(1=>'a', 2=>'b')

    map = BiMap('a'=>1, 'b'=>2)
    @test Dict(map) == Dict('a'=>1, 'b'=>2)
    @test Dict(rev(map)) == Dict(1=>'a', 2=>'b')
end

@testset "rev" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test map == rev(rev(map))
end

@testset "length" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test length(map) == 2
end
@testset "iterate" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test sort(collect(map)) == ['a'=>1, 'b'=>2]
end
@testset "getindex" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test map['a'] == 1
end
@testset "setindex" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test (map['b'] = 99) == 99
    @test map['b'] == 99
    @test rev(map)[99] == 'b'
end
@testset "haskey" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test haskey(map, 'a')
    @test !haskey(map, 'z')
end
@testset "get" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test get(map, 'a', 5) == 1
    @test get(map, 'z', 5) == 5
end
@testset "get_do" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test get(map, 'a') do; 5 end == 1
    @test get(map, 'z') do; 5 end == 5
end
@testset "get!" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test get!(map, 'a', 5) == 1
    @test get!(map, 'z', 5) == 5
    @test map['z'] == 5
    @test rev(map)[5] == 'z'
end
@testset "get!_do" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test get!(map, 'a') do; 5 end == 1
    @test get!(map, 'z') do; 5 end == 5
    @test map['z'] == 5
    @test rev(map)[5] == 'z'
end
@testset "getkey" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test getkey(map, 'a', 'T') == 'a'
    @test getkey(map, 'z', 'T') == 'T'
end
@testset "delete!" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test delete!(map, 'a') === map
    @test !haskey(map, 'a')
    @test !haskey(rev(map), 1)
end
@testset "pop!_default" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test pop!(map, 'b', 99) == 2
    @test !haskey(map, 'b')
    @test !haskey(rev(map), 2)
    @test pop!(map, 'b', 100) == 100
end
@testset "pop!" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test pop!(map, 'b') == 2
    @test !haskey(map, 'b')
    @test !haskey(rev(map), 2)
    @test_throws KeyError pop!(map, 'b')
end
@testset "keys" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test sort!(collect(keys(map))) == ['a', 'b']
end
@testset "values" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test sort!(collect(values(map))) == [1, 2]
end
@testset "pairs" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test sort!(collect(pairs(map))) == ['a'=>1, 'b'=>2]
end
@testset "sizehint!" begin
    map = BiMap('a'=>1, 'b'=>2)
    sizehint!(map, 20)
    # Nothing to test
end
@testset "keytype" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test keytype(map) == Char
end
@testset "valtype" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test valtype(map) == Int
end
@testset "hasval" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test hasval(map, 1)
    @test !hasval(map, 99)
end
@testset "getval" begin
    map = BiMap('a'=>1, 'b'=>2)
    @test getval(map, 2, 99) == 2
    @test getval(map, 9, 99) == 99
end


text:
using DashBootstrapComponents, DashHtmlComponents

spinners = html_div([
    dbc_button(
        dbc_spinner(size="sm"),
        color="primary",
        disabled=true,
        className="mr-1",
    ),
    dbc_button(
        [dbc_spinner(size="sm"), " Loading..."],
        color="primary",
        disabled=true,
    ),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
export RandomizedCD, RandomizedCD2

"""
Randomized coordinate descent

Section 3.7 in Gower, R. M., & Richtárik, P. (2015).
Randomized iterative methods for linear systems.
SIAM Journal on Matrix Analysis and Applications, 36(4), 1660-1690.
"""
function RandomizedCD(stp::AbstractStopping; is_zero_start::Bool = false, kwargs...)
  A, b = get_matrix(stp.pb), get_vector(stp.pb)
  state = stp.current_state
  m, n = size(A)
  if is_zero_start
    state.res .= b
  else
    state.res .= b .- A * state.x
  end
  OK = start!(stp, no_opt_check = true)

  while !OK
    i = Int(floor(rand() * n) + 1)
    Ai = view(A, :, i)

    Aires = @kdot(m, Ai, state.res)
    nAi = @kdot(m, Ai, Ai)
    state.x[i] += Aires / nAi
    state.res .-= Ai .* (Aires / nAi)
    OK = cheap_stop!(stp)
  end

  return stp
end

"""
Randomized coordinate descent for symmetric positive definite matrix

Section 3.4 in Gower, R. M., & Richtárik, P. (2015).
Randomized iterative methods for linear systems.
SIAM Journal on Matrix Analysis and Applications, 36(4), 1660-1690.
"""
function RandomizedCD2(stp::AbstractStopping; is_zero_start::Bool = false, kwargs...)
  A, b = get_matrix(stp.pb), get_vector(stp.pb)
  state = stp.current_state
  n = size(A, 2)
  if is_zero_start
    state.res .= b
  else
    state.res .= b .- A * state.x
  end
  OK = start!(stp, no_opt_check = true)

  while !OK
    i = Int(floor(rand() * n) + 1)
    Ai = view(A, i, :)

    Aix = @kdot(n, Ai, state.x)
    if A[i, i] != 0
      state.x[i] -= (Aix - b[i]) / A[i, i]
    end
    state.res .+= Ai .* ((Aix - b[i]) / A[i, i])

    OK = cheap_stop!(stp)
  end

  return stp
end


text:
# among other things, this file is a workaround to allow simpler
# interactive development without the need to replace Main module
# on each file reload; please, don't merge it into other files

abstract Aggregator  # aggregator of data vectors using custom logic

type VarStats
    n_obs::Int64
    min::Float64
    max::Float64
    var::Float64
    sum::Float64

    VarStats() = new(0, 0, 0, 0)
end

type StatsAggregator <: Aggregator

end

StatsAggregator() = StatsAggregator(Dict{Any, Vector{VarStats}}())

show(io::IO, agg::StatsAggregator) =
    print(io, "StatsAggregator($(keys(agg.groups)))")

type Index
    lookup::Dict{Symbol, Int}
    names::Vector{Symbol}
end

function Index(names::Vector{Symbol})
    lookup = Dict{Symbol, Int}(zip(names, 1:length(names)))
    Index(lookup, names)
end

abstract GroupBy

type StatsGroupBy <: GroupBy
    index::Index
    gvars::Vector{Int}  # indexes of vars to group by
    avars::Vector{Int}  # indexes of vars to aggregate
    agg::StatsAggregator
    groupers::Dict{Symbol, Function}
end

function StatsGroupBy(names::Vector{Symbol}, groupers::Dict{Symbol, Function})
    index = Index(names)
    gnames = filter(x -> in(x, keys(groupers)), names)
    gvars = [index.lookup[name] for name in gnames]
    anames = filter(x -> !in(x, keys(groupers)), names)
    avars = [index.lookup[name] for name in anames]
    agg = StatsAggregator()
    return StatsGroupBy(index, gvars, avars, agg, groupers)
end

function StatsGroupBy(names::Vector{Symbol}, groupby::Vector{Symbol})
    return StatsGroupBy(names, Dict([(name, identity) for name in groupby]))
end

show(io::IO, grb::StatsGroupBy) =
    print(io, "StatsGroupBy($(grb.index.names),$(keys(grb.groupers)))")


text:
struct SegmentedSum{T} <: AggregationFunction
    C::T
end

Flux.@functor SegmentedSum

SegmentedSum(d::Int) = SegmentedSum(zeros(Float32, d))

Base.show(io::IO, sm::SegmentedSum) = print(io, "SegmentedSum($(length(sm.C)))\n")
modelprint(io::IO, sm::SegmentedSum; pad=[]) = paddedprint(io, "SegmentedSum($(length(sm.C)))")

(m::SegmentedSum)(x::MaybeMatrix, bags::AbstractBags, w=nothing) = segmented_sum_forw(x, m.C, bags, w)
function (m::SegmentedSum)(x::AbstractMatrix, bags::AbstractBags, w::AggregationWeights, mask::AbstractVector)
    segmented_sum_forw(x .* mask', m.C, bags, w)
end

segmented_sum_forw(::Missing, C::AbstractVector, bags::AbstractBags, w) = repeat(C, 1, length(bags))
function segmented_sum_forw(x::AbstractMatrix, C::AbstractVector, bags::AbstractBags, w::AggregationWeights)
    y = zeros(eltype(x), size(x, 1), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                y[i, bi] = C[i]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    y[i, bi] += weight(w, i, j) * x[i, j]
                end
            end
        end
    end
    y
end

function segmented_sum_back(Δ, y, x, C, bags, w)
    dx = similar(x)
    dC = zero(C)
    dw = isnothing(w) ? nothing : zero(w)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(C)
                dC[i] += Δ[i, bi]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    dx[i, j] = weight(w, i, j) * Δ[i, bi]
                    ∇dw_segmented_sum!(dw, Δ, x, y, w, i, j, bi)
                end
            end
        end
    end
    dx, dC, nothing, dw
end

function segmented_sum_back(Δ, y, x::Missing, C, bags, w::Nothing)
    dC = zero(C)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(C)
            dC[i] += Δ[i, bi]
        end
    end
    nothing, dC, nothing, nothing
end

∇dw_segmented_sum!(dw::Nothing, Δ, x, y, w::Nothing, i, j, bi) = nothing
function ∇dw_segmented_sum!(dw::AbstractVector, Δ, x, y, w::AbstractVector, i, j, bi)
    dw[j] += Δ[i, bi] * (x[i, j])
end
function ∇dw_segmented_sum!(dw::AbstractMatrix, Δ, x, y, w::AbstractMatrix, i, j, bi)
    dw[i, j] += Δ[i, bi] * (x[i, j])
end

@adjoint function segmented_sum_forw(args...)
    y = segmented_sum_forw(args...)
    grad = Δ -> segmented_sum_back(Δ, y, args...)
    y, grad
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c50a52e2-30fa-11ec-00e9-4f4a81ed30ed
using CSV, DataFrames

# ╔═╡ 9ba9e100-6929-4312-8cd8-75de0cb05a3c
iris = CSV.File("iris_dataset.csv") |> DataFrame

# ╔═╡ b1a1e7fc-58de-45ef-b6c4-355f5e30b7ef
describe(iris)

# ╔═╡ 63668d67-5d58-415d-9fc9-6a891e26ddbd
nrow(iris)

# ╔═╡ 5346265b-2829-4e49-afc7-a016db72db17
unique(iris.Species)

# ╔═╡ 148805ea-c9d6-4bdf-ac79-b9bea8d24f4c
species = groupby(iris, :Species)

# ╔═╡ 7b9ee6d2-95b1-4400-a21d-313b0a97b229
combine(species, nrow)

# ╔═╡ 9f30a14b-eded-4829-bd2a-6b85548ba6b8
combine(species, :SepalLength => minimum => :MinSepalLength)

# ╔═╡ 6e73a20f-860e-4b2e-8dc0-64f1b93f6fa1
colnames = names(iris)

# ╔═╡ 27013be8-56eb-4455-9d2a-53c475d6147a
numeric_columns = colnames[1:4]

# ╔═╡ 4f193ca1-f0a6-472d-a9bf-6deb725c7bb6
transform(iris, numeric_columns .=> (x -> 10x) .=> numeric_columns)

# ╔═╡ 432548a4-67b2-45c9-a0bd-fb3099568f15
begin
	iris[!, :Sample] = 1:nrow(iris)
	long = stack(iris, Not([:Species, :Sample]))
end

# ╔═╡ b94f9bc0-2e63-4cd5-b4b5-b64c41d79270
begin
	long.Dimension = [
		replace(x, r"Sepal|Petal" => "")
		for x in long.variable
		]
	long.Structure = [
		replace(x, r"Length|Width" => "")
		for x in long.variable
		]
	select!(long, Not(:variable))
	unstack(long, :Dimension, :value)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"

[compat]
CSV = "~0.9.10"
DataFrames = "~1.2.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSV]]
deps = ["CodecZlib", "Dates", "FilePathsBase", "InlineStrings", "Mmap", "Parsers", "PooledArrays", "SentinelArrays", "Tables", "Unicode", "WeakRefStrings"]
git-tree-sha1 = "74147e877531d7c172f70b492995bc2b5ca3a843"
uuid = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
version = "0.9.10"

[[CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "d962b5a47b6d191dbcd8ae0db841bc70a05a3f5b"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.13"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[InlineStrings]]
deps = ["Parsers"]
git-tree-sha1 = "19cb49649f8c41de7fea32d089d37de917b553da"
uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
version = "1.0.1"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "a193d6ad9c45ada72c14b731a318bedd3c2f00cf"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.3.0"

[[PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SentinelArrays]]
deps = ["Dates", "Random"]
git-tree-sha1 = "f45b34656397a1f6e729901dc9ef679610bd12b5"
uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
version = "1.3.8"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[WeakRefStrings]]
deps = ["DataAPI", "InlineStrings", "Parsers"]
git-tree-sha1 = "c69f9da3ff2f4f02e811c3323c22e5dfcb584cfa"
uuid = "ea10d353-3f73-51f8-a26c-33c1cb351aa5"
version = "1.4.1"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═c50a52e2-30fa-11ec-00e9-4f4a81ed30ed
# ╠═9ba9e100-6929-4312-8cd8-75de0cb05a3c
# ╠═b1a1e7fc-58de-45ef-b6c4-355f5e30b7ef
# ╠═63668d67-5d58-415d-9fc9-6a891e26ddbd
# ╠═5346265b-2829-4e49-afc7-a016db72db17
# ╠═148805ea-c9d6-4bdf-ac79-b9bea8d24f4c
# ╠═7b9ee6d2-95b1-4400-a21d-313b0a97b229
# ╠═9f30a14b-eded-4829-bd2a-6b85548ba6b8
# ╠═6e73a20f-860e-4b2e-8dc0-64f1b93f6fa1
# ╠═27013be8-56eb-4455-9d2a-53c475d6147a
# ╠═4f193ca1-f0a6-472d-a9bf-6deb725c7bb6
# ╠═432548a4-67b2-45c9-a0bd-fb3099568f15
# ╠═b94f9bc0-2e63-4cd5-b4b5-b64c41d79270
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
const AVAILABLE_THREADS = Base.RefValue{Channel{Int}}()

# Somehow, fetch doesn't do a very good job at preserving
# stacktraces. So, we catch any error in spawn_background
# And return it as a CapturedException, and then use checked_fetch to
# rethrow any exception in that case
function checked_fetch(future)
    value = fetch(future)
    value isa Exception && throw(value)
    return value
end

"""
    spawnbg(f)

Spawn work on any available background thread.
Captures any exception thrown in the thread, to give better stacktraces.

You can use `checked_fetch(spawnbg(f))` to rethrow any exception.

    ** Warning ** this doesn't compose with other ways of scheduling threads
    So, one should use `spawn_background` exclusively in each Julia process.
"""
function spawnbg(f)
    # -1, because we don't spawn on foreground thread 1
    nbackground = Threads.nthreads() - 1
    if nbackground == 0
        # we don't run in threaded mode, so we just run things async
        # to not block forever
        @warn("No threads available, running in foreground thread")
        return @async try
            return f()
        catch e
            # If we don't do this, we get pretty bad stack traces... not sure why!?
            return CapturedException(e, Base.catch_backtrace())
        end
    end
    # Initialize dynamically, could also do this in __init__ but it's nice to keep things in one place
    if !isassigned(AVAILABLE_THREADS)
        # Allocate a Channel with n background threads
        c = Channel{Int}(nbackground)
        AVAILABLE_THREADS[] = c
        # fill queue with available threads
        foreach(i -> put!(c, i + 1), 1:nbackground)
    end
    # take the next free thread... Will block/wait until a thread becomes free
    thread_id = take!(AVAILABLE_THREADS[])

    return ThreadPools.@tspawnat thread_id begin
        try
            return f()
        catch e
            # If we don't do this, we get pretty bad stack traces...
            # not sure why something so basic just doesn't work well \_(ツ)_/¯
            return CapturedException(e, Base.catch_backtrace())
        finally
            # Make thread available again after work is done!
            put!(AVAILABLE_THREADS[], thread_id)
        end
    end
end


text:
export Path
export visitdir

abstract type Path end

"""Represents a filepath"""
struct File <: Path
    path::String
end

"""Represents a dirpath"""
struct Dir <: Path
    path::String
end

function Path(s::AbstractString)
    if isdir(s)
        Dir(s)
    elseif isfile(s)
        File(s)
    else
        error("'$s' is neither a valid directory of file path")
    end
end

"""
    visitdir(fn, dir, args...)

Visits recursively every directory under `dir` and running function
`fn(file, args...)` for every `file` found in each subdirectory.

This means `args...` is an optional list of arguments. You don't have
to provide any. This only visists `.md` files.
"""
function visitdir(fn::Function, dir::AbstractString, args...)
    cd(dir) do
        for entry in readdir()
            if isdir(entry)
                visitdir(fn, entry, args...)
            elseif endswith(entry, ".md")
                fn(entry, args...)
            end
        end
    end
end

max_stars_repo_name: ordovician/AsciiDocTools

max_issues_repo_name: ordovician/AsciiDocTools

max_forks_repo_name: ordovician/AsciiDocTools


text:
abstract type HuffmanTreeNode{T<:SourceSymbol} end

# construct from tuples
HuffmanTreeNode(   t::Tuple{T, <:Real}) where T<:SourceSymbol = HuffmanSymbolNode{T}(t...)
HuffmanTreeNode{T}(t::Tuple{T, <:Real}) where T<:SourceSymbol = HuffmanSymbolNode{T}(t...)
HuffmanTreeNode{T}(t::Tuple{UInt8, <:Real}) where T<:SourceSymbol = HuffmanSymbolNode{T}(t...)
HuffmanTreeNode(t::Tuple{<:HuffmanTreeNode{T}, <:HuffmanTreeNode{T}}) where T<:SourceSymbol =
    HuffmanBranchNode{T}(t...)
HuffmanTreeNode{T}(t::Tuple{<:HuffmanTreeNode{T}, <:HuffmanTreeNode{T}}) where T<:SourceSymbol =
    HuffmanBranchNode{T}(t...)
HuffmanTreeNode{T}(t::Tuple{Tuple, Tuple}) where T<:SourceSymbol =
    HuffmanBranchNode{T}(HuffmanTreeNode{T}(t[1]), HuffmanTreeNode{T}(t[2]))

Base.convert(::Type{T}, t::Tuple) where T <: HuffmanTreeNode = T(t)

Base.show(io::IO, t::HuffmanTreeNode) = print_tree(io, t)
AbstractTrees.nodetype(::HuffmanTreeNode{T}) where T = HuffmanTreeNode{T}

struct HuffmanSymbolNode{T<:SourceSymbol} <: HuffmanTreeNode{T}
    symbol::T
    probability::Float64
    HuffmanSymbolNode{T}(symbol::T, p::Float64=0.0) where T = new(symbol, p)
    HuffmanSymbolNode{T}(symbol::UInt8, p::Float64=0.0) where T = new(symbol, p)
end

const DCSymbolNode = HuffmanSymbolNode{DCSymbol}
const ACSymbolNode = HuffmanSymbolNode{ACSymbol}

HuffmanSymbolNode(symbol::T, p::Float64=0.0) where T = HuffmanSymbolNode{T}(symbol, p)

AbstractTrees.children(s::HuffmanSymbolNode) = ()
AbstractTrees.printnode(io::IO, s::HuffmanSymbolNode) = print(io, s)

function Base.show(io::IO, s::HuffmanSymbolNode{T}) where T <: Union{DCSymbol, ACSymbol}
    compact = get(io, :compact, false)
    if compact
        print(io, T === DCSymbol ? "DCSymbolNode" : "ACSymbolNode")
    else
        print(io, typeof(s))
    end
    print(IOContext(io, :typeinfo=>T), "(", s.symbol, ", ", s.probability, ")")
end

struct HuffmanBranchNode{T<:SourceSymbol} <: HuffmanTreeNode{T}
    n0::HuffmanTreeNode{T} # bit 0
    n1::HuffmanTreeNode{T} # bit 1
    probability::Float64
    HuffmanBranchNode{T}(n0::HuffmanTreeNode{T}, n1::HuffmanTreeNode{T}) where T<:SourceSymbol =
        new(n0, n1, n0.probability + n1.probability)
end
HuffmanBranchNode(n0::HuffmanTreeNode{T}, n1::HuffmanTreeNode{T}) where T<:SourceSymbol =
    HuffmanBranchNode{T}(n0, n1)

AbstractTrees.children(b::HuffmanBranchNode) = (b.n0, b.n1)
AbstractTrees.printnode(io::IO, b::HuffmanBranchNode) = print(io, "・")

struct HuffmanTree{T<:SourceSymbol}
    root::HuffmanTreeNode{T}
end

const DCHuffmanTree = HuffmanTree{DCSymbol}
const ACHuffmanTree = HuffmanTree{ACSymbol}

# The following method creates a huffman tree based on the probabilities.
# Note that the code assignment is not unique and the result is not normalized.
function HuffmanTree(nodes::AbstractVector{S}) where {T, S<:HuffmanTreeNode{T}}
    length(nodes) < 1 && error()
    t = collect(HuffmanTreeNode{T}, nodes)
    while length(t) > 1
        sort!(t, alg=InsertionSort, by=n->n.probability, rev=true)
        # pop the two nodes with the lowest probability
        n1 = pop!(t)
        n0 = pop!(t)
        # push a new branch node with the two nodes as children
        push!(t, HuffmanBranchNode(n0, n1))
    end
    HuffmanTree(t[1])
end

AbstractTrees.children(t::HuffmanTree) = children(t.root)
AbstractTrees.printnode(io::IO, t::HuffmanTree) = print(io, typeof(t))
AbstractTrees.nodetype(::HuffmanTree{T}) where T = HuffmanTreeNode{T}

function Base.show(io::IO, t::HuffmanTree{T}) where T
    function p(io::IO, s::HuffmanSymbolNode{T})
        print(io, "(", s.symbol, ", ", s.probability, ")")
    end
    function p(io::IO, b::HuffmanBranchNode{T})
        print(io, "(")
        p(io, b.n0)
        print(io, ", ")
        p(io, b.n1)
        print(io, ")")
    end
    print(io, typeof(t), "(")
    p(IOContext(io, :compact=>true, :typeinfo=>T), t.root)
    print(io, ")")
end

function Base.show(io::IO, ::MIME"text/plain", t::HuffmanTree{T}) where T
    table = HuffmanTable(t)
    n = 0
    maxdepth = 0
    for li in table.len
        n += li
        n > 20 && break
        maxdepth += 1
    end
    print_tree(IOContext(io, :compact=>true, :typeinfo=>T), t, maxdepth)
end

height(t::HuffmanSymbolNode) = 0
height(t::HuffmanBranchNode) = 1 + max(height(t.n0), height(t.n1))
height(t::HuffmanTree) = height(t.root)

struct HuffmanTable{T<:SourceSymbol}
    id::UInt8 # Th: Huffman table destination identifier
    len::SVector{16, UInt8} # Li
    symbols::AbstractVector{T} # Vij

    function HuffmanTable{T}(len::AbstractVector{<:Integer},
                             symbols::AbstractVector{T};
                             id::Integer=0) where T <: SourceSymbol
        id < 4 || error("`id` must be less than 4.")
        new{T}(id, len, symbols)
    end
end

const DCHuffmanTable = HuffmanTable{DCSymbol}
const ACHuffmanTable = HuffmanTable{ACSymbol}

Base.sizeof(ht::HuffmanTable) = length(ht) + 17
Base.length(ht::HuffmanTable) = length(ht.symbols)

# HuffmanTree <--> HuffmanTable conversions
function HuffmanTree(table::HuffmanTable{T}) where T <: SourceSymbol
    length(table) < 1 && error()
    symbols = collect(zip(table.symbols, codelist(table)))
    # set probability based on the code length
    node(sc) = HuffmanSymbolNode(sc[1], 0.5^length(sc[2]))
    function split(s, depth)
        length(s) == 1 && return node(s[1])
        s0 = filter(sc -> !sc[2][depth], s) # the `depth`-th bit is zero
        s1 = filter(sc ->  sc[2][depth], s) # the `depth`-th bit is one
        n0 = length(s0) == 1 ? node(s0[1]) : split(s0, depth + 1)
        n1 = length(s1) == 1 ? node(s1[1]) : split(s1, depth + 1)
        HuffmanBranchNode(n0, n1)
    end
    HuffmanTree(split(symbols, 1))
end

function HuffmanTable(tree::HuffmanTree{T}, id::Integer=0) where T <: SourceSymbol
    symbols = Tuple{T, Int}[]
    en(t::HuffmanSymbolNode, d::Int) = push!(symbols, (t.symbol, d))
    en(t::HuffmanBranchNode, d::Int) = (en(t.n0, d + 1); en(t.n1, d + 1))
    en(tree.root, 0)
    sort!(symbols, by=(s->s[2]))
    len = UInt8[count(s->s[2]==i, symbols) for i = 1:16]
    HuffmanTable{T}(len, (s->s[1]).(symbols), id=id)
end

function Base.read(io::IO, t::Type{<:HuffmanTable})
    tc, th = u4u4(io)
    len = Vector{UInt8}(undef, 16)
    readbytes!(io, len, 16)
    slen = sum(len)
    symbols = Vector{UInt8}(undef, slen)
    readbytes!(io, symbols, slen)
    t === DCHuffmanTable && tc != 0x0 && error("AC/DC type mismatch")
    t === ACHuffmanTable && tc != 0x1 && error("AC/DC type mismatch")
    if tc == 0x0
        return DCHuffmanTable(len, DCSymbol.(symbols), id=th)
    elseif tc == 0x1
        return ACHuffmanTable(len, ACSymbol.(symbols), id=th)
    else
        error()
    end
end

function codetable(ht::HuffmanTable)
    code = code""
    function u()
        prev = code
        code = inc(code)
        prev
    end
    function v(n)
        code = shl(code)
        CodeBits[u() for j = 1:n]
    end
    SVector{16, Vector{CodeBits}}([v(Int(ht.len[i])) for i = 1:16])
end

codelist(ht::HuffmanTable) = reduce(vcat, codetable(ht), init=[])

function Base.show(io::IO, ht::HuffmanTable)
    print(io, typeof(ht), "([", join(string.(Int.(ht.len)), ", "), "], ")
    print(io, ht.symbols, ", id=", ht.id, ")")
end

function Base.show(io::IO, ::MIME"text/plain", ht::HuffmanTable)
    print(io, length(ht), "-symbol ", typeof(ht), "(id=", ht.id, "):")
    codes = codelist(ht)
    for (s, c) in zip(ht.symbols, codes)
        println(io)
        print(io, " ", s, "  ", c)
    end
end

struct HuffmanDecoder{T<:SourceSymbol}
    symbols::AbstractVector{T}
    codes::AbstractVector{CodeBits}
    index8::SVector{256, Int32} # for codes with length not greater than 8 bits
    indexx::Integer # first index which is out of index8

    function HuffmanDecoder{T}(table::HuffmanTable{T}) where {T<:SourceSymbol}
        symbols = table.symbols
        codes = codelist(table)
        index8 = zeros(Int32, 256)

        length(codes) > 256 && error()
        # e.g.
        # table.len = [0, 1, 5, ..., 0]
        # bitlen=1 -> none
        # bitlen=2 -> code"00"
        # bitlen=3 -> code"010", code"011", code"100", code"101", code"110"
        #
        # 0b00_000000 : 0b00_111111 -> index=1 (w=64, 6 bits)
        #   __ ^^^^^^     __ ^^^^^^
        # 0b010_00000 : 0b010_11111 -> index=2 (w=32, 5 bits)
        # 0b011_00000 : 0b011_11111 -> index=3 (w=32, 5 bits)
        # 0b100_00000 : 0b100_11111 -> index=4 (w=32, 5 bits)
        # 0b101_00000 : 0b101_11111 -> index=5 (w=32, 5 bits)
        # 0b110_00000 : 0b110_11111 -> index=6 (w=32, 5 bits)
        #   ___ ^^^^^     ___ ^^^^^
        index = 1
        p = 0
        for bitlen = 1:8
            w = 0x100 >> bitlen
            for k = 1:table.len[bitlen]
                index8[p+1:p+w] .= index
                index += 1
                p += w
            end
        end
        new{T}(symbols, codes, index8, index)
    end
end

struct BitSequenceError <: Exception end

HuffmanDecoder(table::HuffmanTable{T}) where {T<:SourceSymbol} = HuffmanDecoder{T}(table)

function getcodeindex(d::HuffmanDecoder, bits::UInt32, effective_length::Integer)
    index = d.index8[(bits >> 24) + 1]
    if index > 0 # matched
        c = d.codes[index]
        return effective_length < length(c) ? 0 : index
    end
    for i = d.indexx:length(d.codes)
        c = d.codes[i]
        l = length(c)
        if effective_length < l
            delta = UInt8(l - effective_length)
            (word(c) >> delta) >= (bits >> UInt8(32 - effective_length)) && return 0
            throw(BitSequenceError())
        end
        c == CodeBits(l, bits >> UInt8(32 - l)) && return i
    end
    throw(BitSequenceError())
end


text:
#https://developer.gnome.org/pygtk/stable/class-gtkwidget.html#method-gtkwidget--get-allocation
#https://github.com/JuliaGraphics/Gtk.jl/blob/543a4f13eabf6b62edb0c11c51a715219be37629/src/base.jl

@reexport using Gtk, Cairo
"""
    clear!(task)

Clears all stimuli from window.

# Arguments

- `task`: a task object <: AbstractTask
"""
function clear!(task)
    c = task.canvas
    w = task.width
    @guarded draw(c) do widget
        ctx = getgc(c)
        rectangle(ctx, 0, 0, w, w)
        set_source_rgb(ctx, .8, .8, .8)
        fill(ctx)
    end
    Gtk.showall(c)
    return nothing
end

setup_window(width::Float64, name="") = setup_window(width, width, name)

"""
    setup_window(width::Float64, height::Float64, name="")

Generate a blank window.

# Arguments

- `width`: width of window
- `height`: height of window
- `name=""`: name of window to appear on upper tab
"""
function setup_window(width::Float64, height::Float64, name="")
    canvas = @GtkCanvas()
    window = GtkWindow(canvas, name, width, height)
    Gtk.visible(window, true)
    @guarded draw(canvas) do widget
        ctx = getgc(canvas)
        rectangle(ctx, 0.0, 0.0, width, width)
        set_source_rgb(ctx, .8, .8, .8)
        fill(ctx)
    end
    return canvas,window
end

"""
    draw_attention!(task, actr)

Draws a yellow circle to represent visual attention

# Arguments

- `task`: a task object <: AbstractTask
- `model`: an ACT-R model object
"""
function draw_attention!(task, actr)
    focus = actr.visual.focus
    draw_attention!(task, focus...)
end

"""
    draw_attention!(task, actr)

Draws a yellow circle to represent visual attention

# Arguments

- `task`: a task object <: AbstractTask
- `x`: x coordinate of visual attention
- `y`: y coordinate of visual attention
"""
function draw_attention!(task, x, y)
    c = task.canvas
    @guarded draw(c) do widget
        ctx = getgc(c)
        arc(ctx, x, y, 20, 0, 2π)
        set_source_rgba(ctx, .933, .956, .443, .6)
        fill(ctx)
    end
    Gtk.showall(c)
    return nothing
end

"""
    draw_object!(task, object, x, y)

Draws an object in the window specified at `task.canvas`. The
coordinates of the object are the upper left point of the object

# Arguments

- `task`: a task object <: AbstractTask
- `object`: object to be drawn
- `x`: x coordinate of visual attention
- `y`: y coordinate of visual attention
"""
function draw_object!(task, object, x, y)
    c = task.canvas
    @guarded draw(c) do widget
        ctx = getgc(c)
        select_font_face(ctx, "Arial", Cairo.FONT_SLANT_NORMAL,
             Cairo.FONT_WEIGHT_BOLD);
        set_font_size(ctx, 36)
        set_source_rgb(ctx, 0, 0, 0)
        extents = text_extents(ctx, object)
        x′ = x - (extents[3]/2 + extents[1])
        y′ = y - (extents[4]/2 + extents[2])
        move_to(ctx, x′, y′)
        show_text(ctx, object)
    end
    Gtk.showall(c)
    return nothing
end

function repaint!(task::AbstractTask, actr)
    draw_attention!(task, actr)
end


text:
function warning_callback(message::String)
    global TEST_CALLBACK = true
    error("checking that error are supported")
end

@testset "Errors" begin
    err = ChemfilesError("oops")
    iobuf = IOBuffer()
    show(iobuf, err)
    @test String(iobuf.data[1:(19 + length(err.message))]) == "\"Chemfiles error: oops\""

    Chemfiles.clear_errors()
    @test Chemfiles.last_error() == ""

    remove_chemfiles_warning() do
        @test_throws ChemfilesError Residue(Topology(), 3)
    end
    @test_throws UndefVarError TEST_CALLBACK == false

    @test Chemfiles.last_error() == "Residue index out of bounds in topology: we have 0 residues, but the index is 3"

    Chemfiles.clear_errors()
    @test Chemfiles.last_error() == ""

    Chemfiles.set_warning_callback(warning_callback)
    @test_throws ChemfilesError Residue(Topology(), 3)
    @test TEST_CALLBACK == true

    Chemfiles.set_warning_callback(Chemfiles.__default_warning_callback)
end

@testset "Configuration" begin
    config = joinpath(@__DIR__, "data", "config.toml")
    Chemfiles.add_configuration(config)

    trajectory = joinpath(@__DIR__, "data", "water.xyz")
    frame = read(Trajectory(trajectory))

    @test name(Atom(frame, 9)) == "Oz"
    @test type(Atom(frame, 9)) == "F"
end