text:
################################################################################
# Functions as exported by Clang but removed parameter types.                  #
# This regular expression should help                                          #
# ([A-z_0-9]*)(\{([A-z,_,0-9])*\})*                                          #
################################################################################
function toxav_new(messenger,max_calls)
    ccall((:toxav_new,libtoxcore),Ptr{ToxAv},(Ptr{Tox},Int32),messenger,max_calls)
end

function toxav_kill(av)
    ccall((:toxav_kill,libtoxcore),Void,(Ptr{ToxAv},),av)
end

function toxav_do_interval(av)
    ccall((:toxav_do_interval,libtoxcore),Uint32,(Ptr{ToxAv},),av)
end

function toxav_do(av)
    ccall((:toxav_do,libtoxcore),Void,(Ptr{ToxAv},),av)
end

function toxav_register_callstate_callback(av,cb,id,userdata)
    ccall((:toxav_register_callstate_callback,libtoxcore),Void,(Ptr{ToxAv},ToxAVCallback,ToxAvCallbackID,Ptr{Void}),av,cb,id,userdata)
end

function toxav_register_audio_callback(av,cb,userdata)
    ccall((:toxav_register_audio_callback,libtoxcore),Void,(Ptr{ToxAv},ToxAvAudioCallback,Ptr{Void}),av,cb,userdata)
end

function toxav_register_video_callback(av,cb,userdata)
    ccall((:toxav_register_video_callback,libtoxcore),Void,(Ptr{ToxAv},ToxAvVideoCallback,Ptr{Void}),av,cb,userdata)
end

function toxav_call(av,call_index,friend_id,csettings,ringing_seconds)
    ccall((:toxav_call,libtoxcore),Cint,(Ptr{ToxAv},Ptr{Int32},Cint,Ptr{ToxAvCSettings},Cint),av,call_index,friend_id,csettings,ringing_seconds)
end

function toxav_hangup(av,call_index)
    ccall((:toxav_hangup,libtoxcore),Cint,(Ptr{ToxAv},Int32),av,call_index)
end

function toxav_answer(av,call_index,csettings)
    ccall((:toxav_answer,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{ToxAvCSettings}),av,call_index,csettings)
end

function toxav_reject(av,call_index,reason)
    ccall((:toxav_reject,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{Uint8}),av,call_index,reason)
end

function toxav_cancel(av,call_index,peer_id,reason)
    ccall((:toxav_cancel,libtoxcore),Cint,(Ptr{ToxAv},Int32,Cint,Ptr{Uint8}),av,call_index,peer_id,reason)
end

function toxav_change_settings(av,call_index,csettings)
    ccall((:toxav_change_settings,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{ToxAvCSettings}),av,call_index,csettings)
end

function toxav_stop_call(av,call_index)
    ccall((:toxav_stop_call,libtoxcore),Cint,(Ptr{ToxAv},Int32),av,call_index)
end

function toxav_prepare_transmission(av,call_index,support_video)
    ccall((:toxav_prepare_transmission,libtoxcore),Cint,(Ptr{ToxAv},Int32,Cint),av,call_index,support_video)
end

function toxav_kill_transmission(av,call_index)
    ccall((:toxav_kill_transmission,libtoxcore),Cint,(Ptr{ToxAv},Int32),av,call_index)
end

function toxav_prepare_video_frame(av,call_index,dest,dest_max,input)
    ccall((:toxav_prepare_video_frame,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{Uint8},Cint,Ptr{Cint}),av,call_index,dest,dest_max,input)
end

function toxav_send_video(av,call_index,frame,frame_size)
    ccall((:toxav_send_video,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{Uint8},Uint32),av,call_index,frame,frame_size)
end

function toxav_prepare_audio_frame(av,call_index,dest,dest_max,frame,frame_size)
    ccall((:toxav_prepare_audio_frame,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{Uint8},Cint,Ptr{Int16},Cint),av,call_index,dest,dest_max,frame,frame_size)
end

function toxav_send_audio(av,call_index,frame,size)
    ccall((:toxav_send_audio,libtoxcore),Cint,(Ptr{ToxAv},Int32,Ptr{Uint8},Uint32),av,call_index,frame,size)
end

function toxav_get_peer_csettings(av,call_index,peer,dest)
    ccall((:toxav_get_peer_csettings,libtoxcore),Cint,(Ptr{ToxAv},Int32,Cint,Ptr{ToxAvCSettings}),av,call_index,peer,dest)
end

function toxav_get_peer_id(av,call_index,peer)
    ccall((:toxav_get_peer_id,libtoxcore),Cint,(Ptr{ToxAv},Int32,Cint),av,call_index,peer)
end

function toxav_get_call_state(av,call_index)
    ccall((:toxav_get_call_state,libtoxcore),ToxAvCallState,(Ptr{ToxAv},Int32),av,call_index)
end

function toxav_capability_supported(av,call_index,capability)
    ccall((:toxav_capability_supported,libtoxcore),Cint,(Ptr{ToxAv},Int32,ToxAvCapabilities),av,call_index,capability)
end

function toxav_get_tox(av)
    ccall((:toxav_get_tox,libtoxcore),Ptr{Tox},(Ptr{ToxAv},),av)
end

function toxav_get_active_count(av)
    ccall((:toxav_get_active_count,libtoxcore),Cint,(Ptr{ToxAv},),av)
end

function toxav_add_av_groupchat(tox,audio_callback,userdata)
    ccall((:toxav_add_av_groupchat,libtoxcore),Cint,(Ptr{Tox},Ptr{Void},Ptr{Void}),tox,audio_callback,userdata)
end

function toxav_join_av_groupchat(tox,friendnumber,data,length,audio_callback,userdata)
    ccall((:toxav_join_av_groupchat,libtoxcore),Cint,(Ptr{Tox},Int32,Ptr{Uint8},Uint16,Ptr{Void},Ptr{Void}),tox,friendnumber,data,length,audio_callback,userdata)
end

function toxav_group_send_audio(tox,groupnumber,pcm,samples,channels,sample_rate)
    ccall((:toxav_group_send_audio,libtoxcore),Cint,(Ptr{Tox},Cint,Ptr{Int16},Uint32,Uint8,Uint32),tox,groupnumber,pcm,samples,channels,sample_rate)
end


text:
abstract type Cenum_anonymous <: Cenum end
isanonymous(::Type{<:Cenum_anonymous}) = true
Base.show(io::IO, ::Type{CE}) where {CE<:Cenum_anonymous} = print(io, "<anonymous-enum>")

function Base.show(io::IO, ce::CE) where {CE<:Cenum}
	if !(ce isa get(io, :typeinfo, Nothing))
		show(io, typeof(ce))
		print(io, "(")
	end
	value = convert(eltype(CE), ce)
	found = false
	for (name, val) in Cenumlayout(CE).values
		if value == val
			print(io, "<", name, ">", "(")
			show(io, value)
			print(io, ")")
			found = true
			break
		end
	end
	found || show(io, value)
	if !(ce isa get(io, :typeinfo, Nothing))
		print(io, ")")
	end
end

Base.promote_rule(::Type{T}, ::Type{CE}) where {T<:Integer, CE<:Cenum} = promote_type(T, eltype(CE))

Base.convert(::Type{CE}, x::T) where {CE<:Cenum, T<:Integer} = reinterpret(concrete(CE), convert(eltype(CE), x))
Base.convert(::Type{T}, x::CE) where {CE<:Cenum, T<:Integer} = convert(T, reinterpret(eltype(CE), x))
Base.convert(::Type{CE1}, x::CE2) where {CE1<:Cenum, CE2<:Cenum} = convert(CE1, convert(eltype(CE2), x))

(::Type{CE1})(x::CE2) where {CE1<:Cenum, CE2<:Cenum} = convert(CE1, convert(eltype(CE2), x))
(::Type{CE})(x::T = zero(eltype(CE))) where {CE<:Cenum, T<:Integer} = convert(CE, x)
(::Type{T})(x::CE) where {CE<:Cenum, T<:Integer} = convert(T, x)

Base.eltype(ce::Cenum) = eltype(typeof(ce))
Base.eltype(::Type{CE}) where {CE<:Cenum} = Cenumlayout(CE).type

Base.typemin(::Type{CE}) where {CE<:Cenum} = Cenumlayout(CE).min
Base.typemax(::Type{CE}) where {CE<:Cenum} = Cenumlayout(CE).max

macro cenum(exprs...) return _cenum(__module__, nothing, exprs...) end

function _cenum(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, name::Symbol)
	return _cenum(mod, deps, name, nothing, nothing)
end

function _cenum(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, body::Expr)
	return _cenum(mod, deps, nothing, body, nothing)
end

function _cenum(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, body::Expr, strategy::Symbol)
	return _cenum(mod, deps, nothing, body, strategy)
end

function _cenum(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, name::Symbol, body::Expr)
	return _cenum(mod, deps, name, body, nothing)
end

function _cenum(mod::Module, deps::Union{Vector{Pair{Symbol, Expr}}, Nothing}, name::Union{Symbol, Nothing}, body::Union{Expr, Nothing}, strategy::Union{Symbol, Nothing})
	isnothing(body) || Base.is_expr(body, :braces) || Base.is_expr(body, :bracescat) || error("Expected @cenum to have a `{ ... }` expression for the body of the type, but found `$(body)`")

	strategy = isnothing(strategy) ? :(CBinding.ALIGN_NATIVE) : :(Calignment{$(QuoteNode(Symbol(String(strategy)[3:end-2])))})
	isanon = isnothing(name)
	super = isanon ? :(Cenum_anonymous) : :(Cenum)
	name = isanon ? gensym("anonymous-cenum") : name
	escName = esc(name)
	concreteName = esc(Symbol(name, "\u200B"))  # this is so incredibly evil: appending a 0-width space to the type name

	isOuter = isnothing(deps)
	deps = isOuter ? Pair{Symbol, Expr}[] : deps
	if isnothing(body)
		push!(deps, name => quote
			abstract type $(escName) <: $(super) end
		end)
	else
		fields = []
		values = []
		for arg in body.args
			arg = _expand(mod, deps, arg)

			lastArg = isempty(values) ? :(0) : :($(last(last(values).args)) + 1)
			if arg isa Symbol
				push!(fields, :(const $(arg) = ($(lastArg))::Integer))
				push!(values, :(Pair{$(QuoteNode(arg)), $(arg)}))
			elseif Base.is_expr(arg, :escape, 1) && arg.args[1] isa Symbol
				push!(fields, :(const $(arg) = ($(lastArg))::Integer))
				push!(values, :(Pair{$(QuoteNode(arg.args[1])), $(arg)}))
			elseif Base.is_expr(arg, :(=), 2) && arg.args[1] isa Symbol
				push!(fields, :(const $(arg.args[1]) = ($(arg.args[2]))::Integer))
				push!(values, :(Pair{$(QuoteNode(arg.args[1])), $(arg.args[1])}))
			elseif Base.is_expr(arg, :(=), 2) && Base.is_expr(arg.args[1], :escape, 1) && arg.args[1].args[1] isa Symbol
				push!(fields, :(const $(arg.args[1]) = ($(arg.args[2]))::Integer))
				push!(values, :(Pair{$(QuoteNode(arg.args[1].args[1])), $(arg.args[1])}))
			else
				error("Expected @cenum to have a member symbol (`VALUE`) or member assignment expression (`VALUE = value_expr`) in its body")
			end
		end

		isempty(values) && error("Expected @cenum to have at least 1 value")

		push!(deps, name => quote
			$(fields...)
			abstract type $(escName) <: $(super) end
			primitive type $(concreteName) <: $(escName) sizeof(Cenumlayout(Ctypespec($(super), $(strategy), Tuple{$(values...)})).type)*8 end
			#=
				EnumSpec = Tuple{
					Pair{:name, value}  # each enum constant's name and value
				}
			=#
			CBinding.concrete(::Type{$(escName)}) = $(concreteName)
			CBinding.concrete(::Type{$(concreteName)}) = $(concreteName)
			CBinding.strategy(::Type{$(concreteName)}) = $(strategy)
			CBinding.specification(::Type{$(concreteName)})  = Tuple{$(values...)}
			Base.sizeof(::Type{$(escName)}) = sizeof(CBinding.concrete($(concreteName)))
		end)
	end

	return isOuter ? quote $(map(last, deps)...) ; $(escName) end : escName
end

ext: jl


text:
using Base.Dates

today = today()

Aug16 = Dates.days(Date(2016, 7, 27)-today)/365
Sep16 = Dates.days(Date(2016, 8, 26)-today)/365
Oct16 = Dates.days(Date(2016, 9, 26)-today)/365
Nov16 = Dates.days(Date(2016, 10, 27)-today)/365
Dec16 = Dates.days(Date(2016, 11, 25)-today)/365
Jan17 = Dates.days(Date(2016, 12, 23)-today)/365
Feb17 = Dates.days(Date(2017, 1, 27)-today)/365
Mar17 = Dates.days(Date(2017, 2, 24)-today)/365
Apr17 = Dates.days(Date(2017, 3, 27)-today)/365
May17 = Dates.days(Date(2017, 4, 26)-today)/365
Jun17 = Dates.days(Date(2017, 5, 26)-today)/365
Jul17 = Dates.days(Date(2017, 6, 26)-today)/365
Aug17 = Dates.days(Date(2017, 7, 27)-today)/365
Sep17 = Dates.days(Date(2017, 8, 25)-today)/365
Oct17 = Dates.days(Date(2017, 9, 26)-today)/365
Nov17 = Dates.days(Date(2017, 10, 27)-today)/365
Dec17 = Dates.days(Date(2017, 11, 24)-today)/365
Jan18 = Dates.days(Date(2017, 12, 27)-today)/365
Feb18 = Dates.days(Date(2018, 1, 26)-today)/365
Mar18 = Dates.days(Date(2018, 2, 23)-today)/365
Apr18 = Dates.days(Date(2018, 3, 27)-today)/365
May18 = Dates.days(Date(2018, 4, 26)-today)/365
Jun18 = Dates.days(Date(2018, 5, 25)-today)/365
Jul18 = Dates.days(Date(2018, 6, 26)-today)/365
Aug18 = Dates.days(Date(2018, 7, 27)-today)/365
Sep18 = Dates.days(Date(2018, 8, 27)-today)/365
Oct18 = Dates.days(Date(2018, 9, 26)-today)/365
Nov18 = Dates.days(Date(2018, 10, 26)-today)/365
Dec18 = Dates.days(Date(2018, 11, 26)-today)/365

Q416 = [Oct16, Nov16, Dec16]
Q117 = [Jan17, Feb17, Mar17]
Q217 = [Apr17, May17, Jun17]
Q317 = [Jul17, Aug17, Sep17]
Q417 = [Oct17, Nov17, Dec17]
Q118 = [Jan18, Feb18, Mar18]
Q218 = [Apr18, May18, Jun18]
Q318 = [Jul18, Aug18, Sep18]
Q418 = [Oct18, Nov18, Dec18]

W16 = [Q416, Q117]
S17 = [Q217, Q317]
W17 = [Q417, Q118]
S18 = [Q218, Q318]

# Date(2017,12) - Dates.Day(7)
# Dates.dayofweek(Date(2017,10) - Dates.Day(5))
# Date(2017,7)
#
# Dates.days(Date(2016,7) - Dates.Day(5)
#
# if (Dates.dayofweek(Date(2016,7) - Dates.Day(5)) == 7)


text:
using Test
using StatsModels

@testset "Grouping" begin
    g = Grouping()
    @test isnothing(g.levels)
end

@testset "Grouping pseudo-contrasts" begin
    d = (y = rand(2_000_000), grp=string.([1:1_000_000; 1:1_000_000]))
    ## OOM seems to result in the process being killed on Mac so this messes up CI
    # @test_throws OutOfMemoryError schema(d)
    sch = schema(d, Dict(:grp => Grouping()))
    t = sch[term(:grp)]
    @test t isa CategoricalTerm{Grouping}
    @test size(t.contrasts.matrix) == (0,0)
    @test length(t.contrasts.levels) == 1_000_000
    @test_throws ErrorException StatsModels.modelcols(t, (a = 1.,))

    levs = sort(string.(1:1_000_000))

    @test all(t.contrasts.invindex[lev] == i for (i,lev) in enumerate(levs))
    @test all(t.contrasts.levels[i] == lev for (i,lev) in enumerate(levs))

    # @test_throws OutOfMemoryError fit(MixedModel, @formula(y ~ 1 + (1 | grp)), d)
    @test fit(MixedModel, @formula(y ~ 1 + (1 | grp)), d, contrasts=Dict(:grp => Grouping()), progress=false) isa LinearMixedModel
end


text:
using BinaryProvider
using Compat
using Compat.Test
using SHA

# The platform we're running on
const platform = platform_key()

# Useful command to launch `sh` on any platform
const sh = gen_sh_cmd

# Output of a few scripts we are going to run
const simple_out = "1\n2\n3\n4\n"
const long_out = join(["$(idx)\n" for idx in 1:100], "")
const newlines_out = join(["marco$d\npolo$d\n" for d in ("","\r","\r\n")], "")

# Explicitly probe platform engines in verbose mode to get coverage and make
# CI debugging easier
BinaryProvider.probe_platform_engines!(;verbose=true)

@testset "OutputCollector" begin
    cd("output_tests") do
        # Collect the output of `simple.sh``
        oc = OutputCollector(sh(`./simple.sh`))

        # Ensure we can wait on it and it exited properly
        @test wait(oc)

        # Ensure further waits are fast and still return 0
        let
            tstart = time()
            @test wait(oc)
            @test time() - tstart < 0.1
        end

        # Test that we can merge properly
        @test merge(oc) == simple_out

        # Test that merging twice works
        @test merge(oc) == simple_out

        # Test that `tail()` gives the same output as well
        @test tail(oc) == simple_out

        # Test that colorization works
        let
            red = Base.text_colors[:red]
            def = Base.text_colors[:default]
            gt = "1\n$(red)2\n$(def)3\n4\n"
            @test merge(oc; colored=true) == gt
            @test tail(oc; colored=true) == gt
        end

        # Test that we can grab stdout and stderr separately
        @test stdout(oc) == "1\n3\n4\n"
        @test stderr(oc) == "2\n"
    end

    # Next test a much longer output program
    cd("output_tests") do
        oc = OutputCollector(sh(`./long.sh`))

        # Test that it worked, we can read it, and tail() works
        @test wait(oc)
        @test merge(oc) == long_out
        @test tail(oc; len=10) == join(["$(idx)\n" for idx in 91:100], "")
    end

    # Next, test a command that fails
    cd("output_tests") do
        oc = OutputCollector(sh(`./fail.sh`))

        @test !wait(oc)
        @test merge(oc) == "1\n2\n"
    end

    # Next, test a command that kills itself (NOTE: This doesn't work on windows.  sigh.)
    @static if !iswindows()
        cd("output_tests") do
            oc = OutputCollector(sh(`./kill.sh`))

            @test !wait(oc)
            @test stdout(oc) == "1\n2\n"
        end
    end

    # Next, test reading the output of a pipeline()
    grepline = pipeline(sh(`-c 'printf "Hello\nWorld\nJulia"'`), `grep ul`)
    oc = OutputCollector(grepline)

    @test wait(oc)
    @test merge(oc) == "Julia\n"

    # Next, test that \r and \r\n are treated like \n
    cd("output_tests") do
        oc = OutputCollector(sh(`./newlines.sh`))

        @test wait(oc)
        @test stdout(oc) == newlines_out
    end

    # Next, test that tee'ing to a stream works
    cd("output_tests") do
        ios = IOBuffer()
        oc = OutputCollector(sh(`./simple.sh`); tee_stream=ios, verbose=true)
        @test wait(oc)
        @test merge(oc) == simple_out
        seekstart(ios)
        println(readstring(ios))
    end
end

@testset "PlatformNames" begin
    # Ensure the platform type constructors are well behaved
    @test_throws ArgumentError Linux(:not_a_platform)
    @test_throws ArgumentError MacOS(:i686)
    @test_throws ArgumentError Windows(:armv7l)
    @test_throws ArgumentError Linux(:x86_64, :crazy_libc)

    # Test that our platform_dlext stuff works
    @test platform_dlext(Linux(:x86_64)) == platform_dlext(Linux(:i686))
    @test platform_dlext(Windows(:x86_64)) == platform_dlext(Windows(:i686))
    @test platform_dlext(MacOS()) != platform_dlext(Linux(:armv7l))

    # Test some valid dynamic library paths
    @test valid_dl_path("libfoo.so.1.2.3", Linux(:x86_64))
    @test valid_dl_path("libfoo-1.dll", Windows(:x86_64))
    @test valid_dl_path("libfoo.1.2.3.dylib", MacOS())
    @test !valid_dl_path("libfoo.dylib", Linux(:x86_64))
    @test !valid_dl_path("libfoo.so", Windows(:x86_64))

    # Make sure the platform_key() with explicit triplet works or doesn't
    @test platform_key("x86_64-linux-gnu") == Linux(:x86_64)
    @test platform_key("i686-unknown-linux-gnu") == Linux(:i686)
    @test platform_key("x86_64-apple-darwin14") == MacOS()
    @test platform_key("armv7l-pc-linux-gnueabihf") == Linux(:armv7l)
    @test platform_key("aarch64-unknown-linux-gnu") == Linux(:aarch64)
    @test platform_key("powerpc64le-linux-gnu") == Linux(:ppc64le)
    @test platform_key("x86_64-w64-mingw32") == Windows(:x86_64)
    @test platform_key("i686-w64-mingw32") == Windows(:i686)
    @test_throws ArgumentError platform_key("invalid-triplet-yo")
    @test_throws ArgumentError platform_key("aarch64-unknown-gnueabihf")
    @test_throws ArgumentError platform_key("x86_64-w32-mingw64")

    # Test that we can indeed ask if something is linux or windows, etc...
    @test Compat.Sys.islinux(Linux(:aarch64))
    @test !Compat.Sys.islinux(Windows(:x86_64))
    @test iswindows(Windows(:i686))
    @test !iswindows(Linux(:x86_64))
    @test Compat.Sys.isapple(MacOS())
    @test !Compat.Sys.isapple(Linux(:ppc64le))

    # Test that every supported platform is _something_
    if isdefined(Base.Sys, :isapple)
        isbasesomething(p) = Sys.islinux(p) || Sys.iswindows(p) || Sys.isapple(p)
        @test all(isbasesomething, supported_platforms())
    end
    issomething(p) = Compat.Sys.islinux(p) || iswindows(p) ||
                     Compat.Sys.isapple(p)
    @test all(issomething, supported_platforms())

    @test wordsize(Linux(:i686)) == wordsize(Linux(:armv7l)) == 32
    @test wordsize(MacOS()) == wordsize(Linux(:aarch64)) == 64

    @test triplet(Windows(:i686)) == "i686-w64-mingw32"
    @test triplet(Linux(:x86_64, :musl)) == "x86_64-linux-musl"
    @test triplet(Linux(:armv7l, :musl)) == "arm-linux-musleabihf"
    @test triplet(Linux(:x86_64)) == "x86_64-linux-gnu"
    @test triplet(Linux(:armv7l)) == "arm-linux-gnueabihf"
    @test triplet(MacOS()) == "x86_64-apple-darwin14"
end

@testset "Prefix" begin
    mktempdir() do temp_dir
        prefix = Prefix(temp_dir)

        # Test that it's taking the absolute path
        @test prefix.path == abspath(temp_dir)

        # Test that `bindir()`, `libdir()` and `includedir()` all work
        for dir in unique([bindir(prefix), libdir(prefix), includedir(prefix)])
            @test !isdir(dir)
            mkpath(dir)
        end

        # Create a little script within the bindir to ensure we can run it
        ppt_path = joinpath(bindir(prefix), "prefix_path_test.sh")
        open(ppt_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        chmod(ppt_path, 0o775)

        # Test that activation adds certain paths to our environment variables
        activate(prefix)

        # PATH[1] should be "<prefix>/bin" now
        @test BinaryProvider.split_PATH()[1] == bindir(prefix)
        @test Libdl.DL_LOAD_PATH[1] == libdir(prefix)

        # Test we can run the script we dropped within this prefix.  Once again,
        # something about Windows | busybox | Julia won't pick this up even though
        # the path clearly points to the file.  :(
        @static if !iswindows()
            @test success(sh(`$(ppt_path)`))
            @test success(sh(`prefix_path_test.sh`))
        end

        # Now deactivate and make sure that all traces are gone
        deactivate(prefix)
        @test BinaryProvider.split_PATH()[1] != bindir(prefix)
        @test Libdl.DL_LOAD_PATH[1] != libdir(prefix)
    end
end

@testset "Products" begin
    temp_prefix() do prefix
        # Test that basic satisfication is not guaranteed
        e_path = joinpath(bindir(prefix), "fooifier")
        l_path = joinpath(libdir(prefix), "libfoo.$(Libdl.dlext)")
        e = ExecutableProduct(prefix, "fooifier")
        ef = FileProduct(e_path)
        l = LibraryProduct(prefix, "libfoo")
        lf = FileProduct(l_path)

        @test !satisfied(e; verbose=true)
        @test !satisfied(ef; verbose=true)
        @test !satisfied(l, verbose=true)
        @test !satisfied(lf, verbose=true)

        # Test that simply creating a file that is not executable doesn't
        # satisfy an Executable Product (and say it's on Linux so it doesn't
        # complain about the lack of an .exe extension)
        mkpath(bindir(prefix))
        touch(e_path)
        @test satisfied(ef, verbose=true)
        @static if !iswindows()
            # Windows doesn't care about executable bit, grumble grumble
            @test !satisfied(e, verbose=true, platform=Linux(:x86_64))
        end

        # Make it executable and ensure this does satisfy the Executable
        chmod(e_path, 0o777)
        @test satisfied(e, verbose=true, platform=Linux(:x86_64))

        # Remove it and add a `$(path).exe` version to check again, this
        # time saying it's a Windows executable
        rm(e_path; force=true)
        touch("$(e_path).exe")
        chmod("$(e_path).exe", 0o777)
        @test locate(e, platform=Windows(:x86_64)) == "$(e_path).exe"

        # Test that simply creating a library file doesn't satisfy it if we are
        # testing something that matches the current platform's dynamic library
        # naming scheme, because it must be `dlopen()`able.
        mkpath(libdir(prefix))
        touch(l_path)
        @test satisfied(lf, verbose=true)
        @test !satisfied(l, verbose=true)

        # But if it is from a different platform, simple existence will be
        # enough to satisfy a LibraryProduct
        @static if iswindows()
            l_path = joinpath(libdir(prefix), "libfoo.so")
            touch(l_path)
            @test satisfied(l, verbose=true, platform=Linux(:x86_64))
        else
            l_path = joinpath(libdir(prefix), "libfoo.dll")
            touch(l_path)
            @test satisfied(l, verbose=true, platform=Windows(:x86_64))
        end
    end

    # Ensure that the test suite thinks that these libraries are foreign
    # so that it doesn't try to `dlopen()` them:
    foreign_platform = @static if platform_key() == Linux(:aarch64)
        # Arbitrary architecture that is not dlopen()'able
        Linux(:ppc64le)
    else
        # If we're not Linux(:aarch64), then say the libraries are
        Linux(:aarch64)
    end

    # Test for valid library name permutations
    for ext in ["1.so", "so", "so.1", "so.1.2", "so.1.2.3"]
        temp_prefix() do prefix
            l_path = joinpath(libdir(prefix), "libfoo.$ext")
            l = LibraryProduct(prefix, "libfoo")
            mkdir(dirname(l_path))
            touch(l_path)
            @test satisfied(l; verbose=true, platform=foreign_platform)
        end
    end

    # Test for invalid library name permutations
    for ext in ["so.1.2.3a", "so.1.a"]
        temp_prefix() do prefix
            l_path = joinpath(libdir(prefix), "libfoo.$ext")
            l = LibraryProduct(prefix, "libfoo")
            mkdir(dirname(l_path))
            touch(l_path)
            @test !satisfied(l; verbose=true, platform=foreign_platform)
        end
    end
end

@testset "Packaging" begin
    # Clear out previous build products
    for f in readdir(".")
        if !endswith(f, ".tar.gz") || !endswith(f, ".sha256")
            continue
        end
        rm(f; force=true)
    end

    # Gotta set this guy up beforehand
    tarball_path = nothing
    tarball_hash = nothing

    temp_prefix() do prefix
        # Create random files
        mkpath(bindir(prefix))
        mkpath(libdir(prefix))
        bar_path = joinpath(bindir(prefix), "bar.sh")
        open(bar_path, "w") do f
            write(f, "#!/bin/sh\n")
            write(f, "echo yolo\n")
        end
        baz_path = joinpath(libdir(prefix), "baz.so")
        open(baz_path, "w") do f
            write(f, "this is not an actual .so\n")
        end

        # Next, package it up as a .tar.gz file
        tarball_path, tarball_hash = package(prefix, "./libfoo"; verbose=true)
        @test isfile(tarball_path)

        # Check that we are calculating the hash properly
        tarball_hash_check = open(tarball_path, "r") do f
            bytes2hex(sha256(f))
        end
        @test tarball_hash_check == tarball_hash

        # Test that packaging into a file that already exists fails
        @test_throws ErrorException package(prefix, "./libfoo")
    end

    # Test that we can inspect the contents of the tarball
    contents = list_tarball_files(tarball_path)
    const libdir_name = iswindows() ? "bin" : "lib"
    @test joinpath("bin", "bar.sh") in contents
    @test joinpath(libdir_name, "baz.so") in contents

    # Install it within a new Prefix
    temp_prefix() do prefix
        # Install the thing
        @test install(tarball_path, tarball_hash; prefix=prefix, verbose=true)

        # Ensure we can use it
        bar_path = joinpath(bindir(prefix), "bar.sh")
        baz_path = joinpath(libdir(prefix), "baz.so")

        # Ask for the manifest that contains these files to ensure it works
        manifest_path = manifest_for_file(bar_path; prefix=prefix)
        @test isfile(manifest_path)
        manifest_path = manifest_for_file(baz_path; prefix=prefix)
        @test isfile(manifest_path)

        # Ensure that manifest_for_file doesn't work on nonexistant files
        @test_throws ErrorException manifest_for_file("nonexistant"; prefix=prefix)

        # Ensure that manifest_for_file doesn't work on orphan files
        orphan_path = joinpath(bindir(prefix), "orphan_file")
        touch(orphan_path)
        @test isfile(orphan_path)
        @test_throws ErrorException manifest_for_file(orphan_path; prefix=prefix)

        # Ensure that trying to install again over our existing files is an error
        @test_throws ErrorException install(tarball_path, tarball_path; prefix=prefix)

        # Ensure we can uninstall this tarball
        @test uninstall(manifest_path; verbose=true)
        @test !isfile(bar_path)
        @test !isfile(baz_path)
        @test !isfile(manifest_path)

        # Ensure that we don't want to install tarballs from other platforms
        cp(tarball_path, "./libfoo_juliaos64.tar.gz")
        @test_throws ArgumentError install("./libfoo_juliaos64.tar.gz", tarball_hash; prefix=prefix)
        rm("./libfoo_juliaos64.tar.gz"; force=true)

        # Ensure that hash mismatches throw errors
        fake_hash = reverse(tarball_hash)
        @test_throws ErrorException install(tarball_path, fake_hash; prefix=prefix)
    end

    rm(tarball_path; force=true)
    rm("$(tarball_path).sha256"; force=true)
end

@testset "Verification" begin
    temp_prefix() do prefix
        foo_path = joinpath(prefix, "foo")
        open(foo_path, "w") do file
            write(file, "test")
        end
        foo_hash = bytes2hex(sha256("test"))

        # Check that verifying with the right hash works
        info("This should say; no hash cache found")
        @test verify(foo_path, foo_hash; verbose=true)

        # Check that it created a .sha256 file
        @test isfile("$(foo_path).sha256")

        # Check that it verifies the second time around properly
        info("This should say; hash cache is consistent")
        @test verify(foo_path, foo_hash; verbose=true)

        # Sleep for imprecise filesystems
        sleep(2)

        # Get coverage of messing with different parts of the verification chain
        touch(foo_path)
        info("This should say; file has been modified")
        @test verify(foo_path, foo_hash; verbose=true)
        @test_throws ErrorException verify(foo_path, "0"^32; verbose=true)
        touch(foo_path)
        @test verify(foo_path, foo_hash; verbose=true)
        open("$(foo_path).sha256", "w") do file
            write(file, "this is not the right hash")
        end
        info("This should say; hash has changed")
        @test verify(foo_path, foo_hash; verbose=true)
    end
end

# Use `build_libfoo_tarball.jl` in the BinDeps2.jl repository to generate more of these
const bin_prefix = "https://github.com/staticfloat/small_bin/raw/74b7fd81e3fbc8963b14b0ebbe5421e270d8bdcf"
const libfoo_downloads = Dict(
    Linux(:i686) =>     ("$bin_prefix/libfoo.i686-linux-gnu.tar.gz", "1398353bcbbd88338189ece9c1d6e7c508df120bc4f93afbaed362a9f91358ff"),
    Linux(:x86_64) =>   ("$bin_prefix/libfoo.x86_64-linux-gnu.tar.gz", "b9d57a6e032a56b1f8641771fa707523caa72f1a2e322ab99eeeb011f13ad9f3"),
    Linux(:aarch64) =>  ("$bin_prefix/libfoo.aarch64-linux-gnu.tar.gz", "19d9da0e6e7fb506bf4889eb91e936fda43493a39cd4fd7bd5d65506cede6f95"),
    Linux(:armv7l) =>   ("$bin_prefix/libfoo.arm-linux-gnueabihf.tar.gz", "8e33c1a0e091e6e5b8fcb902e5d45329791bb57763ee9cbcde49c1ec9bd8532a"),
    Linux(:ppc64le) =>  ("$bin_prefix/libfoo.powerpc64le-linux-gnu.tar.gz", "b48a64d48be994ec99b1a9fb60e0af7f4415a57596518cb90a340987b79fad81"),
    MacOS() =>          ("$bin_prefix/libfoo.x86_64-apple-darwin14.tar.gz", "661b71edb433ab334b0fef70db3b5c45d35f2b3bee0d244f54875f1ec899c10f"),
    Windows(:i686) =>   ("$bin_prefix/libfoo.i686-w64-mingw32.tar.gz", "3d4a8d4bf0169007a42d809a1d560083635b1540a1bc4a42108841dcb6d2aaea"),
    Windows(:x86_64) => ("$bin_prefix/libfoo.x86_64-w64-mingw32.tar.gz", "2d08fbc9a534cd021f36b6bbe86ddabb2dafbedeb589581240aa4a8c5b896055"),
)

# Test manually downloading and using libfoo
@testset "Downloading" begin
    temp_prefix() do prefix
        if !haskey(libfoo_downloads, platform)
            warn("Platform $platform does not have a libfoo download, skipping download tests")
        else
            # Test a good download works
            url, hash = libfoo_downloads[platform]
            @test install(url, hash; prefix=prefix, verbose=true)

            fooifier = ExecutableProduct(prefix, "fooifier")
            libfoo = LibraryProduct(prefix, "libfoo")

            @test satisfied(fooifier; verbose=true)
            @test satisfied(libfoo; verbose=true)

            fooifier_path = locate(fooifier)
            libfoo_path = locate(libfoo)

            # We know that foo(a, b) returns 2*a^2 - b
            result = 2*2.2^2 - 1.1

            # Test that we can invoke fooifier
            @test !success(`$fooifier_path`)
            @test success(`$fooifier_path 1.5 2.0`)
            @test parse(Float64,readchomp(`$fooifier_path 2.2 1.1`)) ≈ result

            # Test that we can dlopen() libfoo and invoke it directly
            hdl = Libdl.dlopen_e(libfoo_path)
            @test hdl != C_NULL
            foo = Libdl.dlsym_e(hdl, :foo)
            @test foo != C_NULL
            @test ccall(foo, Cdouble, (Cdouble, Cdouble), 2.2, 1.1) ≈ result
            Libdl.dlclose(hdl)

            # Test uninstallation
            @test uninstall(manifest_from_url(url; prefix=prefix); verbose=true)

            # Test that download_verify_unpack() works
            download_verify_unpack(url, hash, prefix.path)
            @test satisfied(fooifier; verbose=true)
            @test satisfied(libfoo; verbose=true)

            # Test that download_verify twice in a row works, and that mucking
            # with the file causes a redownload if `force` is true:
            tmpfile = joinpath(prefix, "libfoo.tar.gz")
            @test download_verify(url, hash, tmpfile; verbose=true)
            @test download_verify(url, hash, tmpfile; verbose=true)

            # We sleep for at least a second here so that filesystems with low
            # precision in their mtime implementations don't get confused
            sleep(2)

            open(tmpfile, "w") do f
                write(f, "hehehehe")
            end

            @test_throws ErrorException download_verify(url, hash, tmpfile; verbose=true)
            @test download_verify(url, hash, tmpfile; verbose=true, force=true)

        end

        # Test a bad download fails properly
        bad_url = "http://localhost:1/this_is_not_a_file.x86_64-linux-gnu.tar.gz"
        bad_hash = "0"^64
        @test_throws ErrorException install(bad_url, bad_hash; prefix=prefix, verbose=true)
    end
end

# Test the same as the above, but using BinaryPackage abstraction
@testset "BinaryPackage" begin
    temp_prefix() do prefix
        if !haskey(libfoo_downloads, platform)
            warn("Platform $platform does not have a libfoo download, skipping download tests")
        else
            url, hash = libfoo_downloads[platform]
            fooifier = ExecutableProduct(prefix, "fooifier")
            libfoo = LibraryProduct(prefix, "libfoo")
            binpkg = BinaryPackage(url, hash, platform, [fooifier, libfoo])

            # Test installation and uninstallation
            @test install(binpkg; prefix=prefix, verbose=true)
            @test uninstall(binpkg; prefix=prefix, verbose=true)

            # Now test that we can uninstall even if we don't have the right `url`:
            @test install(binpkg; prefix=prefix, verbose=true)
            binpkg2 = BinaryPackage("fakeurl", hash, platform, [fooifier, libfoo])
            @test uninstall(binpkg2; prefix=prefix, verbose=true)

            # Test that we can't uninstall from the wrong prefix
            temp_prefix() do wrong_prefix
                @test_throws ErrorException uninstall(binpkg; prefix=wrong_prefix, verbose=true)
            end

            # Test that we can't guess a manifest path from a package with the wrong
            # url and no products:
            binpkg3 = BinaryPackage("fakeurl", hash, platform)
            @test_throws ErrorException uninstall(binpkg3; prefix=prefix, verbose=true)
        end
    end
end

# Test installation and failure modes of the bundled LibFoo.jl
@testset "LibFoo.jl" begin
    const color="--color=$(Base.have_color ? "yes" : "no")"
    cd("LibFoo.jl") do
        rm("./deps/deps.jl"; force=true)
        rm("./deps/usr"; force=true, recursive=true)

        # Install `libfoo` and build the `deps.jl` file for `LibFoo.jl`
        run(`$(Base.julia_cmd()) $(color) deps/build.jl`)

        # Ensure `deps.jl` was actually created
        @test isfile("deps/deps.jl")
    end

    cd("LibFoo.jl/test") do
        # Now, run `LibFoo.jl`'s tests, adding `LibFoo.jl` to the LOAD_PATH
        # so that the tests can pick up the `LibFoo` module
        withenv("JULIA_LOAD_PATH"=>joinpath(pwd(),"..","src")) do
            run(`$(Base.julia_cmd()) $(color) runtests.jl`)
        end
    end
end


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result = 0
    for i ∈ 1:trunc(Int, log10(n))
        for j ∈ 10 ^ trunc(Int, ceil(i / 2 - 1)):10 ^ trunc(Int, ceil(i / 2)) - 1
            palindrome = parse(Int, string(string(j)[1:end - i % 2], reverse(string(j))))
            binary = string(palindrome, base = 2)
            if binary == reverse(binary)
                result += palindrome
            end
        end
    end
    return result
end

compute(1_000_000)

@benchmark compute(1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
module HermitianDiagonalization

using LinearAlgebra, SparseArrays, LinearMaps
using Requires

export diagonalizer, diagonalizer!, reset!, upper, lower, AbstractEigenMethod,
       Direct, Arpack_IRAM, ArnoldiMethod_IRAM, KrylovKit_IRAM, IterativeSolvers_LOBPCG

abstract type AbstractEigenMethod{Tv} end

function __init__()
    @require Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2" include("pardiso.jl")
    @require KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" include("krylovkit.jl")
    @require Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" include("arpack.jl")
    @require ArnoldiMethod = "ec485272-7323-5ecc-a04f-4719b315124d" include("arnoldimethod.jl")
    @require IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153" include("iterativesolvers.jl")
end

############################################################
# LinearMap
############################################################

struct SpectrumEdge
    upper::Bool
end

const upper = SpectrumEdge(true)
const lower = SpectrumEdge(false)

linearmap(h::AbstractArray, ::SpectrumEdge) = LinearMap(h, ishermitian = true), missing

function linearmap(h::AbstractArray{Tv}, shift::Number) where {Tv}
    fac = lu(h - Tv(shift) * I)
    lmap = let fac = fac
        LinearMap{Tv}((x, y) -> ldiv!(x, fac, y), size(h)...,
                      ismutating = true, ishermitian = true)
    end
    return lmap, fac
end

############################################################
# Diagonalizer
############################################################

mutable struct Diagonalizer{M<:AbstractEigenMethod,Tv,A<:AbstractArray{Tv},L<:LinearMap{Tv},C,E}
    matrix::A
    method::M
    lmap::L
    point::Float64
    engine::E   # Optional support for lmap (e.g. Pardiso solver or factorization)
    codiag::C
end

Base.show(io::IO, d::Diagonalizer{M,Tv}) where {M,Tv} = print(io,
    "Diagonaliser{$M} for $(size(d.matrix)) Hermitian matrix around point $(d.point)")

"""
    diagonalizer(h::AbstractMatrix, method = Direct; point = 0.0, codiag = missing)

Create and initialize a `Diagonalizer` of Hermitian matrix `h` (no need wrap it in
`Hermitian`) with random initial preconditioners. Diagonalization is performed around the
specified `point` in the spectrum using the specified `method`, to choose amongst `Direct`,
`Arpack_IRAM`, `ArnoldiMethod_IRAM`, `KrylovKit_IRAM` and `IterativeSolvers_LOBPCG`. The
corresponding package needs to be loaded (e.g. `using Arpack`) for the method to become
available, except for `Direct` which uses the `eigen` method in `LinearAlgebra`.

See also `ParidoShift` for options to use advanced shift-and-invert methods using tne
MKL Pardiso library.

To compute `nev::Integer` eigenvectors and eigenvalues of `h` using a diagonalizer
`d = diagonalizer(h, ...)` run `d(nev; kw...)`. The result is given in the form of an `Eigen`
object (see LinearAlgebra stdlib). For allowed keywords `kw` see the documentation of the
different `method`s.

# Examples
```jldoctest
julia> using Arpack, SparseArrays

julia> h = sprand(10^3, 10^3, 10^-2); h = h + h'; d = diagonalizer(h, Arpack_IRAM, point = 0.1)
Diagonaliser{Arpack_IRAM{Float64}} for (1000, 1000) Hermitian matrix around point 0.1

julia> d(4).values
4-element Array{Float64,1}:
 0.10165764052779397
 0.09520587982702707
 0.10605045570396891
 0.08574843542711974
```
"""
diagonalizer(h::AbstractArray, ::Type{S} = Direct; kw...) where {S<:AbstractEigenMethod} =
    diagonalizer(h, S(h); kw...)
function diagonalizer(h::AbstractArray{Tv}, method::S,
                      point = 0.0, codiag = missing) where {Tv,S<:AbstractEigenMethod}
    ishermitian(h) || error("Matrix is non-Hermitian")
    lmap, engine = linearmap(h, point)
    return Diagonalizer(h, method, lmap, getpoint(point), engine, codiag)
end

getpoint(point::Number) = Float64(point)
getpoint(p::SpectrumEdge) = p.upper ? Inf : -Inf

"""
    reset!(d::Diagonalizer)

Resets the preconditioners in `d` to a random value.
"""
reset!(d::Diagonalizer{M}) where {M} = (d.method = M(d.matrix); d)

"""
    diagonalizer!(d::Diagonalizer, h; kw...)

Modifies `d` with a new matrix `h` and keyword arguments `kw`, keeping its predonditioners
if they exist. See `diagonalizer` for `kw` options.
"""
function diagonalizer!(d::Diagonalizer, h; kw...)
    e = diagonalizer(h, d.method; kw...)
    d.matrix, d.lmap, d.point, d.engine, d.codiag = e.matrix, e.lmap, e.point, e.engine, e.codiag
    return d
end

############################################################
# Direct diagonalizer
############################################################

struct Direct{Tv} <: AbstractEigenMethod{Tv}
    dense::Matrix{Tv}
end
(::Type{<:Direct})(h::SparseMatrixCSC{Tv}) where {Tv} = Direct{Tv}(Matrix(h))

(d::Diagonalizer{<:Direct})(; kw...) =
    eigen(Hermitian(d.method.dense); sortby = sortfunc(d.point), kw...)

function (d::Diagonalizer{<:Direct})(nev::Integer; kw...)
    e = d(; kw...)
    return Eigen(e.values[1:nev], e.vectors[:, 1:nev])
end

sortfunc(point) = isfinite(point) ? (λ -> abs(λ - point)) : (point > 0 ? reverse : identity)

reset!(d::Diagonalizer{<:Direct}) = d

end # module

max_stars_repo_name: pablosanjose/HermitianDiagonalization

max_issues_repo_name: pablosanjose/HermitianDiagonalization

max_forks_repo_name: pablosanjose/HermitianDiagonalization


text:
#
# types for Simulate.jl
#

"""
    Timing

Enumeration type for scheduling events and timed conditions:

- `at`: schedule an event at a given time,
- `after`: schedule an event a given time after current time,
- `every`: schedule an event every given time from now on,
- `before`: a timed condition is true before a given time,
- `until`: delay until t.
"""
@enum Timing at after every before until

"""
```
SimFunction(func::Function, arg...; kw...)
SF(func::Function, arg...; kw...)
```
Prepare a function for being called as an event in a simulation.

# Arguments
- `func::Function`: function to be executed at a later simulation time
- `arg...`: arguments to the function
- `kw...`: keyword arguments

!!! note
    Be aware that, if the variables stored in a SimFunction are composite types,
    they can change until they are evaluated later by `func`. But that's the nature
    of simulation.

# Example
```jldoctest
julia> using Simulate

julia> f(a,b,c; d=4, e=5) = a+b+c+d+e  # define a function
f (generic function with 1 method)

julia> sf = SimFunction(f, 10, 20, 30, d=14, e=15)  # store it as SimFunction
SimFunction(f, (10, 20, 30), Base.Iterators.Pairs(:d => 14,:e => 15))

julia> sf.func(sf.arg...; sf.kw...)  # and it can be executed later
89

julia> d = Dict(:a => 1, :b => 2) # now we set up a dictionary
Dict{Symbol,Int64} with 2 entries:
  :a => 1
  :b => 2

julia> f(t) = t[:a] + t[:b] # and a function adding :a and :b
f (generic function with 2 methods)

julia> f(d)  # our add function gives 3
3

julia> ff = SimFunction(f, d)   # we set up a SimFunction
SimFunction(f, (Dict(:a => 1,:b => 2),), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}())

julia> d[:a] = 10  # later somehow we need to change d
10

julia> ff  # our SimFunction ff has changed too
SimFunction(f, (Dict(:a => 10,:b => 2),), Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}())

julia> ff.func(ff.arg...; ff.kw...)  # and calling it gives a different result
12
```
"""
struct SimFunction
    func::Function
    arg::Tuple
    kw::Base.Iterators.Pairs

    SimFunction(func, arg...; kw...) = new(func, arg, kw)
end
const SF = SimFunction

"""
    SimExpr = Union{Expr, SimFunction}

A type which is either a `SimFunction` or any Julia expression `Expr`.
"""
SimExpr = Union{Expr, SimFunction}

"""
    sconvert(ex::Union{SimExpr,Array,Tuple})::Array{SimExpr,1}

convert a SimExpr or an array or a tuple of it to an Array{SimExpr,1}
"""
function sconvert(ex::Union{SimExpr,Array,Tuple})::Array{SimExpr,1}
    if isa(ex, SimExpr)
        return convert(Array{SimExpr,1}, [ex])
    elseif isa(ex, Array)
        return convert(Array{SimExpr,1}, ex)
    else
        return convert(Array{SimExpr,1}, [i for i in ex])
    end
end

"""
```
SimEvent(ex::Array{SimExpr, 1}, scope::Module, t::Float64, Δt::Float64)
```
Create a simulation event: a SimExpr or an array of SimExpr to be
executed at event time.

# Arguments
- `ex::Array{SimExpr, 1}`: an array of SimExpr to be evaluated at event time,
- `scope::Module`: evaluation scope,
- `t::Float64`: event time,
- `Δt::Float64`: repeat rate with which the event gets repeated.
"""
struct SimEvent
    "expression or SimFunction to be evaluated at event time"
    ex::Array{SimExpr, 1}
    "evaluation scope"
    scope::Module
    "event time"
    t::Float64
    "repeat time"
    Δt::Float64

    SimEvent(ex::Array{SimExpr, 1}, scope::Module, t::Number, Δt::Number) =
        new(ex, scope, t, Δt)
end

"""
    SimCond(cond::Array{SimExpr, 1}, ex::Array{SimExpr, 1}, scope::Module)

create a condition to be evaluated repeatedly with expressions or functions
to be executed if conditions are met.

# Arguments
- `cond::Array{SimExpr, 1}`: Expr or SFs to be evaluated as conditions
- `ex::Array{SimExpr, 1}`: Expr or SFs to be evaluated if conditions are all true
- `scope::Module`: evaluation scope
"""
struct SimCond
    "Expr or SFs to be evaluated as conditions"
    cond::Array{SimExpr, 1}
    "Expr or SFs to be evaluated if all conditions are True"
    ex::Array{SimExpr, 1}
    "evaluation scope"
    scope::Module
end

"""
    Sample(ex::SimExpr, scope::Module)

Create a sampling expression.

# Arguments
- `ex::SimExpr`: expression or SimFunction to be called at sample time
- `scope::Module`: evaluation scope
"""
struct Sample
    "expression or function to be called at sample time"
    ex::SimExpr
    "evaluation scope"
    scope::Module
end

"""
    SimException(ev::SEvent, value=nothing)

Define a SimException, which can be thrown to processes.

# Parameters
- `ev::SEvent`: delivers an event to the interrupted task
- `value=nothing`: deliver some other value
"""
struct SimException <: Exception
  ev::Any
  value::Any
  SimException(ev::SEvent, value=nothing) = new(ev, value)
end

"""
```
SimProcess( id, func::Function, arg...; kw...)
SP(id, func::Function, arg...; kw...)
```
Prepare a function to run as a process in a simulation.

# Arguments
- `id`: some unique identification, it should get registered with
- `func::Function`: a function `f(arg...; kw...)`
- `arg...`: further arguments to `f`
- `kw...`: keyword arguments to `f`

!!! note
    A function as a SimProcess most often runs in a loop. It has to
    give back control by e.g. doing a `take!(input)` or by calling
    `delay!` etc., which will `yield` it. Otherwise it will starve everything else!

# Examples
```jldoctest
julia> using Simulate

```
"""
mutable struct SimProcess
    id
    task
    sim
    state::SState
    func::Function
    arg::Tuple
    kw::Base.Iterators.Pairs

    SimProcess( id, func::Function, arg...; kw...) =
        new(id, nothing, nothing, Undefined(), func, arg, kw)
end
const SP = SimProcess

"""
```
Clock(Δt::Number=0; t0::Number=0, unit::FreeUnits=NoUnits)
```
Create a new simulation clock.

# Arguments
- `Δt::Number=0`: time increment. If no Δt is given, the simulation doesn't tick,
    but jumps from event to event. Δt can be set later with `sample_time!`.
- `t0::Number=0`: start time for simulation
- `unit::FreeUnits=NoUnits`: clock time unit. Units can be set explicitely by
    setting e.g. `unit=minute` or implicitly by giving Δt as a time or else setting
    t0 to a time, e.g. `t0=60s`.

# Examples
```jldoctest
julia> using Simulate, Unitful

julia> import Unitful: s, minute, hr

julia> c = Clock()
Clock: state=Simulate.Undefined(), time=0.0, unit=, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.0
julia> init!(c)
Simulate.Idle()
julia> c = Clock(1s, unit=minute)
Clock: state=Simulate.Undefined(), time=0.0, unit=minute, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.016666666666666666
julia> c = Clock(1s)
Clock: state=Simulate.Undefined(), time=0.0, unit=s, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=1.0
julia> c = Clock(t0=60s)
Clock: state=Simulate.Undefined(), time=60.0, unit=s, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=0.0
julia> c = Clock(1s, t0=1hr)
Clock: state=Simulate.Undefined(), time=3600.0, unit=s, events: 0, cevents: 0, processes: 0, sampling: 0, sample rate Δt=1.0
```
"""
mutable struct Clock <: SEngine
    "clock state"
    state::SState
    "clock time"
    time::Float64
    "time unit"
    unit::FreeUnits
    "scheduled events"
    events::PriorityQueue{SimEvent,Float64}
    "conditional events"
    cevents::Array{SimCond,1}
    "registered processes"
    processes::Dict{Any, SimProcess}
    "process lock"
    lock::ReentrantLock
    "end time for simulation"
    end_time::Float64
    "event evcount"
    evcount::Int64
    "sample count"
    scount::Int64
    "next event time"
    tev::Float64

    "sampling time, timestep between ticks"
    Δt::Float64
    "Array of sampling expressions to evaluate at each tick"
    sexpr::Array{Sample,1}
    "next sample time"
    tsa::Float64

    function Clock(Δt::Number=0;
                   t0::Number=0, unit::FreeUnits=NoUnits)
        if 1unit isa Time
            Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
            t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
        elseif Δt isa Time
            unit = Unitful.unit(Δt)
            t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
            Δt = Δt.val
        elseif t0 isa Time
            unit = Unitful.unit(t0)
            t0 = t0.val
        else
            nothing
        end
        new(Undefined(), t0, unit, PriorityQueue{SimEvent,Float64}(), SimCond[],
            Dict{Any, SimProcess}(), ReentrantLock(), t0, 0, 0, t0, Δt, Sample[],
            t0 + Δt)
    end
end

function show(io::IO, sim::Clock)
    s1::String = "Clock: "
    s2::String = "state=$(sim.state), "
    s3::String = "time=$(sim.time), "
    s4::String = "unit=$(sim.unit), "
    s5::String = "events: $(length(sim.events)), "
    s6::String = "cevents: $(length(sim.cevents)), "
    s7::String = "processes: $(length(sim.processes)), "
    s8::String = "sampling: $(length(sim.sexpr)), "
    s9::String = "sample rate Δt=$(sim.Δt)"
    print(io, s1 * s2 * s3 * s4 * s5 * s6 * s7 * s8 * s9)
end


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

"""
    WeightedSampler(size, [weights]; replace=false)

A method for weighted sampling from spatial objects that produces
samples of given `size` based on `weights` with or without replacement
depending on the option `replace`. By default weights are uniform.
"""
struct WeightedSampler{W} <: AbstractSampler
  size::Int
  weights::W
  replace::Bool
end

WeightedSampler(size, weights=nothing; replace=false) =
  WeightedSampler(size, weights, replace)

function sample(object, sampler::WeightedSampler)
  n = nelms(object)
  s = sampler.size
  w = sampler.weights
  r = sampler.replace
  if s > n && r == false
    @error "invalid sample size for sampling without replacement"
  end
  ws = isnothing(w) ? fill(1/n, n) : collect(w)
  @assert length(ws) == n "invalid number of weights for object"

  view(object, sample(1:n, Weights(ws), s, replace=r))
end


text:
using BenchmarkTools
using CUDA
using Yota

const SUITE = BenchmarkGroup()

#################################### TRACE ####################################

SUITE["trace"] = BenchmarkGroup()

# x = rand(8, 8)

# const TRACE_FUNCTIONS = [
#     "simple" => x -> sum(x),
#     "loop" => x -> (for _=1:100 x = x'x end; x),
# ]

# for (name, f) in TRACE_FUNCTIONS
#     SUITE["trace"][name] = @benchmarkable Yota.trace($f, $x) samples=10 seconds=60
# end

################################### HELPERS ###################################

# SUITE["helpers"] = BenchmarkGroup()

# x = rand(4, 5)
# dx = zero(x)
# dy = ones(4, 3)

# SUITE["helpers"]["ungetindex!"] = @benchmarkable Yota.ungetindex!($dx, $x, $dy, :, $[1, 3, 1])
# SUITE["helpers"]["ungetindex"] = @benchmarkable Yota.ungetindex($x, $dy, :, $[1, 3, 1])

# if CUDA.functional()
#     x = rand(4, 5) |> cu
#     dx = zero(x) |> cu
#     dy = ones(4, 3) |> cu

#     SUITE["helpers"]["cu:ungetindex!"] = @benchmarkable Yota.ungetindex!($dx, $x, $dy, :, $(cu([1, 3, 1])))
#     SUITE["helpers"]["cu:ungetindex"] = @benchmarkable Yota.ungetindex($x, $dy, :, $(cu([1, 3, 1])))
# end


text:
# Simple static kd-trees.

abstract type KDNode end

struct KDLeafNode <: KDNode
end

struct KDInternalNode{T <: AbstractFloat} <: KDNode
    j::Int             # dimension on which the data is split
    med::T             # median value where the split occours
    leftnode::KDNode
    rightnode::KDNode
end

struct KDTree{T <: AbstractFloat}
    xs::AbstractMatrix{T} # A matrix of n, m-dimensional observations
    perm::Vector{Int}     # permutation of data to avoid modifying xs
    root::KDNode          # root node
    verts::Set{Vector{T}}
    bounds::Matrix{T}     # Top-level bounding box
end

"""
    KDTree(xs, leaf_size_factor, leaf_diameter factor)

Construct a kd-tree

Args:
  - `xs`: an `n` x `m` matrix containing `n`, `m`-dimensional observations.
  - `leaf_size_factor`: Stop spliting if a node contains
      fewer than `leaf_size_factor * n` elements.
  - `leaf_diameter_factor`: Stop spliting if a node's bounding
      hypercube has diameter less that `leaf_diameter_factor`
      times the diameter of the root node's bounding hypercube.

Returns:
  A `KDTree` object

"""
function KDTree(xs::AbstractMatrix{T},
                leaf_size_factor=0.05,
                leaf_diameter_factor=0.0) where T <: AbstractFloat

    n, m = size(xs)
    perm = collect(1:n)

    bounds = Array{T}(undef, 2, m)
    for j in 1:m
        col = xs[:,j]
        bounds[1, j] = minimum(col)
        bounds[2, j] = maximum(col)
    end

    diam = diameter(bounds)

    leaf_size_cutoff = ceil(Int, leaf_size_factor * n)
    leaf_diameter_cutoff = leaf_diameter_factor * diam
    verts = Set{Vector{T}}()

    # Add a vertex for each corner of the hypercube
    for vert in Iterators.product([bounds[:,j] for j in 1:m]...)
        push!(verts, T[vert...])
    end

    root = build_kdtree(xs, perm, bounds, leaf_size_cutoff, leaf_diameter_cutoff, verts)

    KDTree(xs, collect(1:n), root, verts, bounds)
end

"""
    diameter(bounds)

Compute the diamater of a hypercube (i.e. the maximum distance between any 2 points in the hypercube) defined by `bounds`.

Args:
  - `bounds`: A 2 by `n` matrix where `bounds[1,i]` gives the lower bound in
      the `i`th dimension and `bounds[2,j]` the upper bound.

Returns:
  Computed diameter

"""
function diameter(bounds::Matrix)
    euclidean(vec(bounds[1,:]), vec(bounds[2,:]))
end

"""
    build_kdtree(xs, perm, bounds, leaf_size_cutoff, leaf_diameter_cutoff, verts)

Recursively build a kd-tree

Args:
  - `xs`: Data being orginized.
  - `perm`: Permutation of the data, used to avoid
      directly sorting or modifying `xs`.
  - `bounds`: Bounding hypercube of the node.
  - `leaf_size_cutoff`: stop spliting on nodes with more than
      this many values.
  - `leaf_diameter_cutoff`: stop splitting on nodes with less
       than this diameter.
  - `verts`: current set of vertexes

Modifies:
  `perm`, `verts`

Returns:
  Either a `KDLeafNode` or a `KDInternalNode`
"""
function build_kdtree(xs::AbstractMatrix{T},
                      perm::AbstractArray,
                      bounds::Matrix{T},
                      leaf_size_cutoff::Int,
                      leaf_diameter_cutoff::T,
                      verts::Set{Vector{T}}) where T
    n, m = size(xs)

    if length(perm) <= leaf_size_cutoff || diameter(bounds) <= leaf_diameter_cutoff
        return KDLeafNode()
    end

    # split on the dimension with the largest spread
    # maxspread, j = findmax(maximum(xs[perm, k]) - minimum(xs[perm, k]) for k in 1:m)
    j = 1
    maxspread = 0
    for k in 1:m
        xmin = Inf
        xmax = -Inf
        for i in perm
            xmin = min(xmin, xs[i, k])
            xmax = max(xmax, xs[i, k])
        end
        if xmax - xmin > maxspread
            maxspread = xmax - xmin
            j = k
        end
    end

    # find the median and partition
    if isodd(length(perm))
        mid = length(perm) ÷ 2
        partialsort!(perm, mid, by=i -> xs[i, j])
        med = xs[perm[mid], j]
        mid1 = mid
        mid2 = mid + 1
    else
        mid1 = length(perm) ÷ 2
        mid2 = mid1 + 1
        partialsort!(perm, mid1:mid2, by=i -> xs[i, j])
        med = (xs[perm[mid1], j] + xs[perm[mid2], j]) / 2
    end

    leftbounds = copy(bounds)
    leftbounds[2, j] = med
    leftnode = build_kdtree(xs, view(perm,1:mid1), leftbounds,
                            leaf_size_cutoff, leaf_diameter_cutoff, verts)

    rightbounds = copy(bounds)
    rightbounds[1, j] = med
    rightnode = build_kdtree(xs, view(perm,mid2:length(perm)), rightbounds,
                             leaf_size_cutoff, leaf_diameter_cutoff, verts)

    coords = Array{Array}(undef, m)
    for i in 1:m
        if i == j
            coords[i] = [med]
        else
            coords[i] = bounds[:, i]
        end
    end

    for vert in Iterators.product(coords...)
        push!(verts, T[vert...])
    end

    KDInternalNode{T}(j, med, leftnode, rightnode)
end

"""
    bounds_verts(bounds)

Given a bounding hypecube `bounds`, return its verticies
"""
function bounds_verts(bounds::Matrix)
    collect(Iterators.product([bounds[:, i] for i in 1:size(bounds, 2)]...))
end

"""
    traverse(kdtree, x)

Traverse the tree `kdtree` to the bottom and return the verticies of
the bounding hypercube of the leaf node containing the point `x`.
"""
function traverse(kdtree::KDTree{T}, x::AbstractVector{T}) where T
    m = size(kdtree.bounds, 2)

    if length(x) != m
        throw(DimensionMismatch("$(m)-dimensional kd-tree searched with a length $(length(x)) vector."))
    end

    for j in 1:m
        if x[j] < kdtree.bounds[1, j] || x[j] > kdtree.bounds[2, j]
            error(
                  """
                  Loess cannot perform extrapolation. Predict can only be applied
                  to points within the bounding hypercube of the data used to train
                  the model.
                  """)
        end
    end

    bounds = copy(kdtree.bounds)
    node = kdtree.root
    while !isa(node, KDLeafNode)
        if x[node.j] <= node.med
            bounds[2, node.j] = node.med
            node = node.leftnode
        else
            bounds[1, node.j] = node.med
            node = node.rightnode
        end
    end

    bounds_verts(bounds)
end


text:
using DifferentialEquations
using Distributions
using Plots
using DelimitedFiles

#-----------------------------------
#specify the signaling topology
#-----------------------------------

topology = [0 1 1
            0 0 1
            0 0 0]

nCell = size(topology, 1)

#-----------------------------------
#specify the number of iterations
#-----------------------------------

nits = 20

#-----------------------------------
# Range of A0 values
#-----------------------------------

A0_values = 0.68:0.01:0.74

#-----------------------------------
#specify the values of lambda and kappa
#-----------------------------------

lambda = 1.0
kappa = 1.0

#-------------------------------------------------
#specify the mean wait time and pulse time
#-------------------------------------------------

mean = 50.0

pulse = 5.0

#-----------------------------------
# internal ODE system
#-----------------------------------

function GataPu1(du,u,p,t)
    a = p[1:5]
    b = p[6:8]
    g = p[9:17]
    L = p[18:20] #L = [A,B,C]
    du[1] = -b[1]*u[1] + (a[1]*L[1] + a[2]*u[1])/(1 + g[1]*L[1] + g[2]*u[1] + g[3]*u[1]*u[2])
    du[2] = -b[2]*u[2] + (a[3]*L[2] + a[4]*u[2])/(1 + g[4]*L[2] + g[5]*u[2] + g[6]*u[1]*u[2] + g[7]*u[1]*u[3])
    du[3] = -b[3]*u[3] + (a[5]*u[1])/(1 + g[8]*u[1] + g[9]*L[3])
end

#-----------------------------------
# parameter vector p
#-----------------------------------

p = zeros(0)
a = [1.0,0.25,1.0,0.25,0.01]
b = [0.01,0.01,0.01]
g = [1.0,0.25,1.0,1.0,0.25,1.0,0.13,0.01,10]
L = [1.0,0.5,0] #L = [A,B,C] = [A,0.5,0]
append!(p,a)
append!(p,b)
append!(p,g)
append!(p,L)

#-----------------------------------
# wait time distribution
#-----------------------------------

function rand_exponential(mean)
    if mean <= 0.0
        error("mean must be positive")
    end
    -mean*log(rand())
end

function wait_time(mean)
    round(rand_exponential(mean))
end

#-----------------------------------
#Simulation Code
#-----------------------------------

GHighCounts = zeros(length(A0_values), nCell+1)

for l = 1:length(A0_values)
    A0 = A0_values[l]
    for j = 1:nits

        CellValues = [[20.0, 1.0, 20.0, A0] for i=1:nCell, n=1:1]
        CellIts = zeros(nCell)
        CellWaitTimes = [wait_time(mean) for i=1:nCell, n=1:1]
        CellState = ["wait" for i=1:nCell, n=1:1]
        PulseCoeffs = zeros(Float64, nCell, nCell)

        while min(CellIts...) < 1000 + 40*nCell
            t = min(CellWaitTimes...)

            for i = 1:nCell
                u0 = [CellValues[i][1], CellValues[i][2], CellValues[i][3]]
                p[18] = CellValues[i][4]
                tspan = (0.0, t)
                prob = ODEProblem(GataPu1,u0,tspan,p)
                sol = solve(prob)
                CellValues[i][1] = sol[length(sol)][1]
                CellValues[i][2] = sol[length(sol)][2]
                CellValues[i][3] = sol[length(sol)][3]
            end

            CellWaitTimes = CellWaitTimes .- t

            for i = 1:nCell
                if CellWaitTimes[i] == 0
                    if CellState[i] == "wait"
                        for m = 1:nCell
                            if topology[i,m] == 1
                                PulseCoeffs[i,m] = (CellValues[i][1]/(lambda + CellValues[i][2]))
                                CellValues[m][4] = CellValues[m][4] * PulseCoeffs[i,m]
                            elseif topology[i,m] == -1
                                PulseCoeffs[i,m] = (CellValues[i][2]/(kappa + CellValues[i][1]))
                                CellValues[m][4] = CellValues[m][4] * PulseCoeffs[i,m]
                            end
                        end
                        CellWaitTimes[i] = pulse
                        CellState[i] = "pulse"
                        CellIts[i] = CellIts[1] + 1
                    else
                        for m = 1:nCell
                            if topology[i,m] == 1
                                CellValues[m][4] = CellValues[m][4] / PulseCoeffs[i,m]
                            elseif topology[i,m] == -1
                                CellValues[m][4] = CellValues[m][4] / PulseCoeffs[i,m]
                            end
                        end

                        CellWaitTimes[i] = wait_time(mean)
                        CellState[i] = "wait"
                    end
                end
            end
        end

        for i = 1:nCell
            if CellValues[i][1] > 30
                GHighCounts[l, i+1] = GHighCounts[l, i+1] + 1
                GHighCounts[l, 1] = A0
            end
        end
    end
    for i = 1:nCell
        GHighCounts[l, i+1] = GHighCounts[l, i+1]/nits
    end
end

writedlm("../results/Topology_$nCell"*"_cells.txt", GHighCounts, ", ")

max_stars_repo_name: alavendelm/Cell-Cell-Communication

max_issues_repo_name: alavendelm/Cell-Cell-Communication

max_forks_repo_name: alavendelm/Cell-Cell-Communication


text:
"""
    spzeros(A, E, B, C, D; fast = false, atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the (finite and infinite) Smith zeros of the structured matrix pencil `M-λN`

              | A-λE | B |
     M - λN = |------|---|
              |  C   | D |

in `val`, information on the multiplicities of infinite zeros in `iz` and
information on the Kronecker-structure in the `KRInfo` object.

The information on the multiplicities of infinite zeros is provided in the vector `iz`,
where each `i`-th element `iz[i]` is equal to `k-1`, where `k` is the order of an infinite elementary divisor with `k > 0`.
The number of infinite zeros contained in `val` is the sum of the components of `iz`.

The information on the Kronecker-structure consists of the right Kronecker indices `rki`, left Kronecker indices `lki`,
infinite elementary divisors `id`, the number of finite eigenvalues `nf`, the normal rank `nrank`, and can be obtained
from `KRInfo` as `KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`pkstruct`](@ref).

The computation of the zeros is performed by reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF)
which exhibits the splitting of the infinite and finite eigenvalue,
the multiplicities of the infinite eigenvalues, the left and right Kronecker indices and the normal rank.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguements `atol1`, `atol2`  and `rtol` specify the absolute tolerance for the nonzero
elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
"""
function spzeros(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,UniformScaling{Bool},Missing},
    B::Union{AbstractVecOrMat,Missing}, C::Union{AbstractMatrix,Missing}, D::Union{AbstractVecOrMat,Missing};
    fast::Bool = true, atol1::Real = ismissing(A) ? zero(real(eltype(D))) : zero(real(eltype(A))),
    atol2::Real = ismissing(A) ? zero(real(eltype(D))) : zero(real(eltype(A))),
    rtol::Real = (ismissing(A) ? 1 : min(size(A)...))*eps(real(float(ismissing(A) ? one(real(eltype(D))) : one(real(eltype(A))))))*iszero(min(atol1,atol2)))

   # Step 0: Reduce to the standard form
   M, N, Q, Z, n, m, p = sreduceBF(A, E, B, C, D, atol = atol2, rtol = rtol, fast = fast, withQ = false, withZ = false)
   mM, nM = size(M)
   T = eltype(M)

   maxmn = max(mM,nM)
   rki = Vector{Int}(undef,maxmn)
   id = Vector{Int}(undef,maxmn)
   lki = Vector{Int}(undef,maxmn)

   tol1 = max(atol1, rtol*opnorm(M,1))
   mrinf = 0
   nrinf = 0
   niz = 0
   i = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      ired = mrinf+1:mM
      jred = nrinf+1:nM
      τ, ρ = _preduce1!(n, m, p, view(M,ired,jred), view(N,ired,jred), Q, Z, tol1;
                        fast = fast, withQ = false, withZ = false)
      i += 1
      mui = ρ+τ
      rki[i] = m - mui
      id[i]  = τ
      mrinf += mui
      nrinf += m
      i > 1 && (niz += τ*(i-1))
      n -= ρ
      m = ρ
      p -= τ
   end

   r = mrinf + n

   rtrail = 0
   ctrail = 0
   j = 0
   while p > 0
      # Step 3: Particular case of the dual PREDUCE algorithm
      ired = mrinf+1:mM-rtrail
      jred = nrinf+1:nM-ctrail
      ρ = _preduce4!(n, 0, p, view(M,ired,jred),view(N,ired,jred), Q, Z, tol1, fast = fast, withQ = false, withZ = false)
      j += 1
      rtrail += p
      ctrail += ρ
      lki[j] = p - ρ
      n -= ρ
      p = ρ
   end
   if1 = mrinf+1:mrinf+n
   jf1 = nrinf+1:nrinf+n
   return [eigvalsnosort!(view(M,if1,jf1),view(N,if1,jf1)); Inf*ones(real(T),niz) ], minf(id[2:i]),
           KRInfo(kroni(rki[1:i]), kroni(lki[1:j]), minf(id[1:i]), n, r)
end
"""
    speigvals(A, E, B, C, D; fast = false, atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the (finite and infinite) eigenvalues of the structured matrix pencil `M-λN`

              | A-λE | B |
     M - λN = |------|---|
              |  C   | D |

in `val` and information on the  Kronecker-structure in the `KRInfo` object.

The information on the  Kronecker-structure consists of the right Kronecker indices `rki`, left Kronecker indices `lki`,
infinite elementary divisors `id`, the number of finite eigenvalues `nf`, the normal rank `nrank`, and can be obtained
from `KRInfo` as `KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`pkstruct`](@ref).

The computation of the eigenvalues
is performed by reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF)
which exhibits the splitting of the infinite and finite eigenvalue,
the multiplicities of the infinite eigenvalues, the left and right Kronecker indices and the normal rank.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguements `atol1`, `atol2`  and `rtol` specify the absolute tolerance for the nonzero
elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
"""
function speigvals(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,UniformScaling{Bool},Missing},
   B::Union{AbstractVecOrMat,Missing}, C::Union{AbstractMatrix,Missing}, D::Union{AbstractVecOrMat,Missing};
   fast::Bool = true, atol1::Real = ismissing(A) ? zero(real(eltype(D))) : zero(real(eltype(A))),
   atol2::Real = ismissing(A) ? zero(real(eltype(D))) : zero(real(eltype(A))),
   rtol::Real = (ismissing(A) ? 1 : min(size(A)...))*eps(real(float(ismissing(A) ? one(real(eltype(D))) : one(real(eltype(A))))))*iszero(min(atol1,atol2)))

   # Step 0: Reduce to the standard form
   M, N, Q, Z, n, m, p = sreduceBF(A, E, B, C, D, atol = atol2, rtol = rtol, fast = fast, withQ = false, withZ = false)
   mM, nM = size(M)
   T = eltype(M)

   maxmn = max(mM,nM)
   rki = Vector{Int}(undef,maxmn)
   id = Vector{Int}(undef,maxmn)
   lki = Vector{Int}(undef,maxmn)

   tol1 = max(atol1, rtol*opnorm(M,1))
   mrinf = 0
   nrinf = 0
   ni = 0
   i = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      ired = mrinf+1:mM
      jred = nrinf+1:nM
      τ, ρ = _preduce1!(n, m, p, view(M,ired,jred), view(N,ired,jred), Q, Z, tol1;
                        fast = fast, withQ = false, withZ = false)
      i += 1
      mui = ρ+τ
      rki[i] = m - mui
      id[i]  = τ
      mrinf += mui
      nrinf += m
      ni += τ*i
      n -= ρ
      m = ρ
      p -= τ
   end

   r = mrinf + n

   rtrail = 0
   ctrail = 0
   j = 0
   while p > 0
      # Step 3: Particular case of the dual PREDUCE algorithm
      ired = mrinf+1:mM-rtrail
      jred = nrinf+1:nM-ctrail
      ρ = _preduce4!(n, 0, p, view(M,ired,jred),view(N,ired,jred), Q, Z, tol1, fast = fast, withQ = false, withZ = false)
      rtrail += p
      ctrail += ρ
      j += 1
      lki[j] = p - ρ
      n -= ρ
      p = ρ
   end
   if1 = mrinf+1:mrinf+n
   jf1 = nrinf+1:nrinf+n
   return [eigvalsnosort!(view(M,if1,jf1),view(N,if1,jf1)); Inf*ones(real(T),ni) ],
          KRInfo(kroni(rki[1:i]), kroni(lki[1:j]), minf(id[1:i]), n, r)
end
"""
    spkstruct(A, E, B, C, D; fast = false, atol1::Real = 0, atol2::Real = 0, rtol::Real=min(atol1,atol2)>0 ? 0 : n*ϵ) -> KRInfo

Determine the Kronecker-structure information of the structured matrix pencil `M-λN`

              | A-λE | B |
     M - λN = |------|---|
              |  C   | D |

and return an `KRInfo` object.

The information on the  Kronecker-structure consists of the right Kronecker indices `rki`, left Kronecker indices `lki`,
infinite elementary divisors `id`, the number of finite eigenvalues `nf`, the normal rank `nrank`, and can be obtained
from `KRInfo` as `KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`pkstruct`](@ref).

The right Kronecker indices are provided in the integer vector `rki`, where each `i`-th element `rki[i]` is
the row dimension `k` of an elementary Kronecker block of size `k x (k+1)`.
The number of elements of `rki` is the dimension of
the right nullspace of the pencil `M-λN` and their sum is the least degree of a right polynomial nullspace basis.

The left Kronecker indices are provided in the integer vector `lki`, where each `i`-th element `lki[i]` is
the column dimension `k` of an elementary Kronecker block of size `(k+1) x k`.
The number of elements of `lki` is the dimension of
the left nullspace of the pencil `M-λN` and their sum is the least degree of a left polynomial nullspace basis.

The multiplicities of infinite eigenvalues are provided in the integer vector `id`, where each `i`-th element `id[i]` is
the order of an infinite elementary divisor (i.e., the multiplicity of an infinite eigenvalue).

The determination of the Kronecker-structure information is performed by reducing the pencil `M-λN` to an
appropriate Kronecker-like form (KLF) which exhibits the number of finite eigenvalues,
the multiplicities of the infinite eigenvalues, the left and right Kronecker indices and the normal rank.
The reduction is performed using orthonal similarity transformations and involves rank decisions based
on rank revealing QR-decompositions with column pivoting, if `fast = true`, or, the more reliable,
SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguements `atol1`, `atol2`  and `rtol` specify the absolute tolerance for the nonzero
elements of `M`, the absolute tolerance for the nonzero elements of `N`, and the relative tolerance for the nonzero elements of `M` and `N`, respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
"""
function spkstruct(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,UniformScaling{Bool},Missing},
   B::Union{AbstractVecOrMat,Missing}, C::Union{AbstractMatrix,Missing}, D::Union{AbstractVecOrMat,Missing};
   fast::Bool = true, atol1::Real = ismissing(A) ? zero(real(eltype(D))) : zero(real(eltype(A))),
   atol2::Real = ismissing(A) ? zero(real(eltype(D))) : zero(real(eltype(A))),
   rtol::Real = (ismissing(A) ? 1 : min(size(A)...))*eps(real(float(ismissing(A) ? one(real(eltype(D))) : one(real(eltype(A))))))*iszero(min(atol1,atol2)))

   # Step 0: Reduce to the standard form
   M, N, Q, Z, n, m, p = sreduceBF(A, E, B, C, D, atol = atol2, rtol = rtol, fast = fast, withQ = false, withZ = false)
   mM, nM = size(M)
   T = eltype(M)

   maxmn = max(mM,nM)
   rki = Vector{Int}(undef,maxmn)
   id = Vector{Int}(undef,maxmn)
   lki = Vector{Int}(undef,maxmn)

   tol1 = max(atol1, rtol*opnorm(M,1))
   mrinf = 0
   nrinf = 0
   ni = 0
   i = 0
   while m > 0
      # Steps 1 & 2: Standard algorithm PREDUCE
      ired = mrinf+1:mM
      jred = nrinf+1:nM
      τ, ρ = _preduce1!(n, m, p, view(M,ired,jred), view(N,ired,jred), Q, Z, tol1;
                        fast = fast, withQ = false, withZ = false)
      i += 1
      mui = ρ+τ
      rki[i] = m - mui
      id[i]  = τ
      mrinf += mui
      nrinf += m
      n -= ρ
      m = ρ
      p -= τ
   end

   r = mrinf + n

   rtrail = 0
   ctrail = 0
   j = 0
   while p > 0
      # Step 3: Particular case of the dual PREDUCE algorithm
      ired = mrinf+1:mM-rtrail
      jred = nrinf+1:nM-ctrail
      ρ = _preduce4!(n, 0, p, view(M,ired,jred),view(N,ired,jred), Q, Z, tol1, fast = fast, withQ = false, withZ = false)
      j += 1
      rtrail += p
      ctrail += ρ
      lki[j] = p - ρ
      n -= ρ
      p = ρ
   end
   return KRInfo(kroni(rki[1:i]), kroni(lki[1:j]), minf(id[1:i]), n, r)
end
"""
    sprank(A, E, B, C, D; fastrank = true, atol1 = 0, atol2 = 0, rtol = min(atol1,atol2)>0 ? 0 : n*ϵ)

Compute the normal rank of the structured matrix pencil `M - λN`

              | A-λE | B |
     M - λN = |------|---|.
              |  C   | D |

If `fastrank = true`, the rank is evaluated by counting how many singular values of `M - γ N` have magnitude greater than `max(max(atol1,atol2), rtol*σ₁)`,
where `σ₁` is the largest singular value of `M - γ N` and `γ` is a randomly generated value. If `fastrank = false`,
the rank is evaluated as `nr + ni + nf + nl`, where `nr` and `nl` are the sums of right and left Kronecker indices,
respectively, while `ni` and `nf` are the number of infinite and finite eigenvalues, respectively. The sums `nr+ni` and
`nf+nl`, are determined from an appropriate Kronecker-like form (KLF) exhibiting the spliting of the right and left structures
of the pencil `M - λN`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of `M`, the absolute tolerance for the nonzero elements of `N`,  and the relative tolerance
for the nonzero elements of `M` and `N`. The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `M`, and `ϵ` is the
machine epsilon of the element type of `M`.
For efficiency purpose, the reduction to the relevant KLF is only partially performed
using rank decisions based on rank revealing SVD-decompositions.

!!! compat "Julia 1.1"
    The use of `atol` and `rtol` keyword arguments in rank determinations requires at least Julia 1.1.
    To enforce compatibility with Julia 1.0, the newer function rank in Julia 1.1 has been explicitly included.
"""
function sprank(A::Union{AbstractMatrix,Missing}, E::Union{AbstractMatrix,UniformScaling{Bool},Missing},
   B::Union{AbstractVecOrMat,Missing}, C::Union{AbstractMatrix,Missing}, D::Union{AbstractVecOrMat,Missing};
   fastrank::Bool = true, atol1::Real = ismissing(A) ? (ismissing(D) ? zero(1.) : zero(real(eltype(D)))) : zero(real(eltype(A))),
   atol2::Real = ismissing(A) ? (ismissing(D) ? zero(1.) : zero(real(eltype(D)))) : zero(real(eltype(A))),
   rtol::Real = (ismissing(A) ? 1 : min(size(A)...))*eps(real(float(ismissing(A) ? (ismissing(D) ? one(1.) : one(real(eltype(D)))) : one(real(eltype(A))))))*iszero(min(atol1,atol2)))

   if fastrank
      xor(ismissing(A),ismissing(E)) && error("A and E must be both either present or missing")
      ismissing(A) && !ismissing(B) && error("B can not be present if A is missing")
      ismissing(A) && !ismissing(C) && error("C can not be present if A is missing")
      !ismissing(D) && !ismissing(B) && ismissing(C)  && error("D can not be present if C is missing")
      !ismissing(D) && !ismissing(C) && ismissing(B)  && error("D can not be present if B is missing")
      eident = (typeof(E) == UniformScaling{Bool})
      if ismissing(A) && ismissing(D)
         return 0
      elseif ismissing(A)
         return rank(D,atol = atol1, rtol = rtol)
      elseif ismissing(B) && ismissing(C)
         isa(A,Adjoint) && (A = copy(A))
         isa(E,Adjoint) && (E = copy(E))
         ndx, nx = size(A)
         eident || (ndx,nx) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
         eident ? T = eltype(A) : T = promote_type(eltype(A), eltype(E))
         T <: BlasFloat || (T = promote_type(Float64,T))
         nu = 0
         B = zeros(T,ndx,nu)
         ny = 0
         C = zeros(T,ny,nx)
         D = zeros(T,ny,nu)
      elseif ismissing(B)
         isa(A,Adjoint) && (A = copy(A))
         isa(E,Adjoint) && (E = copy(E))
         ndx, nx = size(A)
         eident || (ndx,nx) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
         nx == size(C,2) || throw(DimensionMismatch("A and C must have the same number of columns"))
         T = promote_type(eltype(A), eltype(E), eltype(C) )
         T <: BlasFloat || (T = promote_type(Float64,T))
         nu = 0
         B = zeros(T,ndx,nu)
         ny = size(C,1)
         D = zeros(T,ny,nu)
      elseif ismissing(C)
         isa(A,Adjoint) && (A = copy(A))
         isa(E,Adjoint) && (E = copy(E))
         ndx, nx = size(A)
         eident || (ndx,nx) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
         ndx1, nu = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
         ndx == ndx1 || throw(DimensionMismatch("A and B must have the same number of rows"))
         T = promote_type(eltype(A), eltype(E), eltype(B), eltype(E), eltype(A), eltype(E), )
         T <: BlasFloat || (T = promote_type(Float64,T))
         ny = 0
         C = zeros(T,ny,nx)
         D = zeros(T,ny,nu)
      else
         isa(A,Adjoint) && (A = copy(A))
         isa(E,Adjoint) && (E = copy(E))
         ndx, nx = size(A)
         T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
         eident || (ndx,nx) == size(E) || throw(DimensionMismatch("A and E must have the same dimensions"))
         ny, nu = typeof(D) <: AbstractVector ? (length(D),1) : size(D)
         ndx1, nu1 = typeof(B) <: AbstractVector ? (length(B),1) : size(B)
         (ndx,nu) == (ndx1, nu1) || throw(DimensionMismatch("A, B and D must have compatible dimensions"))
         (ny,nx) == size(C) || throw(DimensionMismatch("A, C and D must have compatible dimensions"))
         T = promote_type(eltype(A), eltype(E), eltype(B), eltype(C), eltype(D))
         T <: BlasFloat || (T = promote_type(Float64,T))
      end

      eident || (T = promote_type(T,eltype(E)))
      T <: BlasFloat || (T = promote_type(Float64,T))
      if eident
         ndx == nx || throw(DimensionMismatch("A must be a square matrix"))
      end
      # fast finishing if possible
      max(ndx,nx,nu,ny) == 0 && (return 0)

      (!ismissing(A) && eltype(A) != T) && (A = convert(Matrix{T},A))
      (!eident && !ismissing(E) && eltype(E) != T) && (E = convert(Matrix{T},E))
      (!ismissing(B) && eltype(B) != T) && (B = convert(Matrix{T},B))
      (!ismissing(C) && eltype(C) != T) && (C = convert(Matrix{T},C))
      (!ismissing(D) && eltype(D) != T) && (D = convert(Matrix{T},D))
      if eident
         scale = opnorm(A,1)*rand()
         return rank([A+scale*E B; C D], atol = atol1, rtol = rtol)
      else
         nrmN = opnorm(E,1)
         nrmN == zero(nrmN) && (return rank([A B; C D], atol = atol1, rtol = rtol))
         nrmM = max(opnorm(A,1), typeof(B) <: AbstractVector ? norm(B,1) : opnorm(B,1), opnorm(C,Inf), typeof(D) <: AbstractVector ? norm(D,1) : opnorm(D,1))
         scale = nrmM/nrmN*rand()
         return rank([A+scale*E B; C D], atol = max(atol1,atol2), rtol = rtol)
      end
   else
      M, N, Q, Z, n, m, p = sreduceBF(A, E, B, C, D, atol = atol2, rtol = rtol,
                                      fast = false, withQ = false, withZ = false)
      mM, nM = size(M)

      n == min(mM,nM) && (return n)
      tol1 = max(atol1, rtol*opnorm(M,1))
      mrinf = 0
      nrinf = 0
      while m > 0
         # Steps 1 & 2: Standard algorithm PREDUCE
         ired = mrinf+1:mM
         jred = nrinf+1:nM
         τ, ρ = _preduce1!(n, m, p, view(M,ired,jred), view(N,ired,jred), Q, Z, tol1;
                           fast = false, withQ = false, withZ = false)
         mrinf += ρ+τ
         nrinf += m
         n -= ρ
         m = ρ
         p -= τ
      end
      return mrinf+n
   end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Wigxjpf"
version = v"1.11.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("http://fy.chalmers.se/subatom/wigxjpf/wigxjpf-latest.tar.gz", "5d078bbbf87c917d0df3c5e2204cbd8a51041517630ac84bdca728611ea2f12f")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/wigxjpf-1.11

mkdir -p $WORKSPACE/destdir/lib
mkdir -p $WORKSPACE/destdir/bin
mkdir -p $WORKSPACE/destdir/shared/licenses

if [ "$(uname)" == "Darwin" ]; then
    make lib/libwigxjpf_shared.dylib
    cp lib/libwigxjpf_shared.dylib $WORKSPACE/destdir/lib/
else
    make lib/libwigxjpf_shared.so
    cp lib/libwigxjpf_shared.so $WORKSPACE/destdir/lib/
    cp lib/libwigxjpf_shared.so $WORKSPACE/destdir/bin/libwigxjpf_shared.dll
fi

cp README $WORKSPACE/destdir/shared/licenses/LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(;experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libwigxjpf_shared", :wigxjpf)
]

# Dependencies that must be installed before this package can be built
dependencies = [Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Helper data for color parsing
include("names_data.jl")

# Color Parsing
# -------------

const col_pat_hex  = r"^\s*(#|0x)([[:xdigit:]]{3,8})\s*$"
const col_pat_rgb  = r"^\s*rgb\(\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*\)\s*$"
const col_pat_hsl  = r"^\s*hsl\(\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*\)\s*$"
const col_pat_rgba = r"^\s*rgba?\(\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*[,/]\s*((?:\d+|(?=\.\d))(?:\.\d*)?%?)\s*\)\s*$"
const col_pat_hsla = r"^\s*hsla?\(\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*[,\s]\s*(\d+%?)\s*[,/]\s*((?:\d+|(?=\.\d))(?:\.\d*)?%?)\s*\)\s*$"

chop1(x) = SubString(x, 1, lastindex(x) - 1) # `chop` is slightly slow

# Parse a number used in the "rgb()" or "hsl()" color.
function parse_rgb(num::AbstractString)
    if @inbounds num[end] == '%'
        return N0f8(clamp(parse(Int, chop1(num), base=10) / 100, 0, 1))
    else
        v = clamp(parse(Int, num, base=10), 0, 255)
        return reinterpret(N0f8, unsafe_trunc(UInt8, v))
    end
end

function parse_hsl_hue(num::AbstractString)
    if @inbounds num[end] == '%'
        error("hue cannot end in %")
    else
        return parse(Int, num, base=10)
    end
end

function parse_hsl_sl(num::AbstractString)
    if @inbounds num[end] != '%'
        error("saturation and lightness must end in %")
    else
        return parse(Int, chop1(num), base=10) / 100
    end
end

# Parse a number used in the alpha field of "rgba()" and "hsla()".
function parse_alpha_num(num::AbstractString)
    if @inbounds num[end] == '%'
        return parse(Int, chop1(num), base=10) / 100f0
    else
        # `parse(Float32, num)` is somewhat slow on Windows(x86_64-w64-mingw32).
        # However, the following has the opposite effect on Linux.
        # m = match(r"0?\.(\d{1,9})", num)
        # if m != nothing
        #     d = m.captures[1]
        #     return parse(Int, d, base=10) / Float32(exp10(length(d)))
        # end
        return parse(Float32, num)
    end
end

function _parse_colorant(desc::String)
    n0f8(x) = reinterpret(N0f8, unsafe_trunc(UInt8, x))
    mat = match(col_pat_hex, desc)
    if mat !== nothing
        prefix, len = mat.captures[1], length(mat.captures[2])
        digits = parse(UInt32, mat.captures[2], base=16)
        if len == 6
            return convert(RGB{N0f8}, reinterpret(RGB24, digits))
        elseif len == 3
            return RGB(n0f8((digits>>8) & 0xF * 0x11),
                       n0f8((digits>>4) & 0xF * 0x11),
                       n0f8((digits>>0) & 0xF * 0x11))
        elseif len == 8
            if prefix[1] == '0'
                return ARGB{N0f8}(n0f8(digits>>16),
                                  n0f8(digits>> 8),
                                  n0f8(digits>> 0),
                                  n0f8(digits>>24))
            else
                return RGBA{N0f8}(n0f8(digits>>24),
                                  n0f8(digits>>16),
                                  n0f8(digits>> 8),
                                  n0f8(digits>> 0))
            end
        elseif len == 4
            if prefix[1] == '0'
                return ARGB{N0f8}(n0f8((digits>> 8) & 0xF * 0x11),
                                  n0f8((digits>> 4) & 0xF * 0x11),
                                  n0f8((digits>> 0) & 0xF * 0x11),
                                  n0f8((digits>>12) & 0xF * 0x11))
            else
                return RGBA{N0f8}(n0f8((digits>>12) & 0xF * 0x11),
                                  n0f8((digits>> 8) & 0xF * 0x11),
                                  n0f8((digits>> 4) & 0xF * 0x11),
                                  n0f8((digits>> 0) & 0xF * 0x11))
            end
        end
    end
    mat = match(col_pat_rgb, desc)
    if mat !== nothing
        return RGB{N0f8}(parse_rgb(mat.captures[1]),
                         parse_rgb(mat.captures[2]),
                         parse_rgb(mat.captures[3]))
    end

    mat = match(col_pat_hsl, desc)
    if mat !== nothing
        T = ColorTypes.eltype_default(HSL)
        return HSL{T}(parse_hsl_hue(mat.captures[1]),
                      parse_hsl_sl(mat.captures[2]),
                      parse_hsl_sl(mat.captures[3]))
    end

    mat = match(col_pat_rgba, desc)
    if mat !== nothing
        return RGBA{N0f8}(parse_rgb(mat.captures[1]),
                          parse_rgb(mat.captures[2]),
                          parse_rgb(mat.captures[3]),
                          parse_alpha_num(mat.captures[4]))
    end

    mat = match(col_pat_hsla, desc)
    if mat !== nothing
        T = ColorTypes.eltype_default(HSLA)
        return HSLA{T}(parse_hsl_hue(mat.captures[1]),
                       parse_hsl_sl(mat.captures[2]),
                       parse_hsl_sl(mat.captures[3]),
                       parse_alpha_num(mat.captures[4]))
    end

    sdesc = strip(desc)
    c = get(color_names, sdesc, nothing)
    c !== nothing && return RGB{N0f8}(n0f8(c[1]), n0f8(c[2]), n0f8(c[3]))

    # since `lowercase` is slightly slow, it is applied only when needed
    ldesc = lowercase(sdesc)
    c = get(color_names, ldesc, nothing)
    c !== nothing && return RGB{N0f8}(n0f8(c[1]), n0f8(c[2]), n0f8(c[3]))

    ldesc == "transparent" && return RGBA{N0f8}(0,0,0,0)

    wo_spaces = replace(ldesc, r"(?<=[^ ]{3}) (?=[^ ]{3})" => "")
    c = get(color_names, wo_spaces, nothing)
    if c !== nothing
        camel = replace(titlecase(ldesc), " " => "")
        Base.depwarn(
            """
            The X11 color names with spaces are not recommended because they are not allowed in the SVG/CSS.
            Use "$camel" or "$wo_spaces" instead.
            """, :parse)
        return RGB{N0f8}(n0f8(c[1]), n0f8(c[2]), n0f8(c[3]))
    end

    error("Unknown color: ", desc)
end

"""
    parse(Colorant, desc)

Parse a color description.

This parses a subset of HTML/CSS color specifications. In particular, everything
is supported but: `currentColor`.

It does support named colors (though it uses X11 named colors, which are
slightly different than W3C named colors in some cases), `rgb()`, `hsl()`,
`#RGB`, and `#RRGGBB` syntax.

# Arguments

- `Colorant`: literal Colorant
- `desc`: color name or description

A literal Colorant will parse according to the `desc` string (usually returning
an `RGB`); any more specific choice will return a color of the specified type.

# Returns

- an `RGB{N0f8}` color, or

- an `HSL` color if `hsl(h, s, l)` was used

- an `RGBA` color if `rgba(r, g, b, a)` was used

- an `HSLA` color if `hsla(h, s, l, a)` was used

- an `ARGB{N0f8}` color if `0xAARRGGBB`/`0xARGB` was used

- a specific `Colorant` type as specified in the first argument

!!! note "Note for X11 named colors"
    The X11 color names with spaces (e.g. "sea green") are not recommended
    because they are not allowed in the SVG/CSS.

!!! note "Note for hex notations"
    You can parse not only the CSS-style hex notations `#RRGGBB`/`#RGB`, but
    also `0xRRGGBB`/`0xRGB`.

    You can also parse the 8-digit or 4-digit hex notation into an RGB color
    with alpha. However, the result depends on the prefix (i.e. `#` or `0x`).
    ```@example
    julia> parse(Colorant, "#FF8800AA") # transparent orange
    RGBA{N0f8}(1.0,0.533,0.0,0.667)

    julia> parse(Colorant, "0xFF8800AA") # opaque purple
    ARGB{N0f8}(0.533,0.0,0.667,1.0)
    ```
"""
function Base.parse(::Type{C}, desc::AbstractString) where {C<:Colorant}
    c = _parse_colorant(String(desc))
    C === Colorant && return c
    c = convert(C, c)
    return isconcretetype(C) ? c::C : c
end

Base.parse(::Type{C}, desc::Symbol) where {C<:Colorant} = parse(C, string(desc))
Base.parse(::Type{C}, c::Colorant) where {C<:Colorant} = c

"""
    @colorant_str(ex)

Parse a literal color name as a Colorant.
"""
macro colorant_str(ex)
    isa(ex, AbstractString) || error("colorant requires literal strings")
    col = parse(Colorant, ex)
    :($col)
end

@noinline function ColorTypes.color(str::AbstractString)
    Base.depwarn("color(\"$str\") is deprecated, use colorant\"$str\" or parse(Colorant, \"$str\")", :color)
    parse(Colorant, str)
end


text:
export Agent

import Functors: functor
using Setfield: @set

"""
    Agent(;kwargs...)

A wrapper of an `AbstractPolicy`. Generally speaking, it does nothing but to
update the trajectory and policy appropriately in different stages.

# Keywords & Fields

- `policy`::[`AbstractPolicy`](@ref): the policy to use
- `trajectory`::[`AbstractTrajectory`](@ref): used to store transitions between an agent and an environment
"""
Base.@kwdef struct Agent{P<:AbstractPolicy,T<:AbstractTrajectory} <: AbstractPolicy
    policy::P
    trajectory::T
end

functor(x::Agent) = (policy = x.policy,), y -> @set x.policy = y.policy

(agent::Agent)(env) = agent.policy(env)

function check(agent::Agent, env::AbstractEnv)
    if ActionStyle(env) === FULL_ACTION_SET &&
       !haskey(agent.trajectory, :legal_actions_mask)
        #     @warn "The env[$(nameof(env))] is of FULL_ACTION_SET, but I can not find a trace named :legal_actions_mask in the trajectory"
    end
    check(agent.policy, env)
end

Base.nameof(agent::Agent) = nameof(agent.policy)

#####
# Default behaviors
#####

"""
Here we extend the definition of `(p::AbstractPolicy)(::AbstractEnv)` in
`RLBase` to accept an `AbstractStage` as the first argument. Algorithm designers
may customize these behaviors respectively by implementing:

- `(p::YourPolicy)(::AbstractStage, ::AbstractEnv)`
- `(p::YourPolicy)(::PreActStage, ::AbstractEnv, action)`

The default behaviors for `Agent` are:

1. Update the inner `trajectory` given the context of `policy`, `env`, and
   `stage`.
  1. By default we do nothing.
  2. In `PreActStage`, we `push!` the current **state** and the **action** into
     the `trajectory`.
  3. In `PostActStage`, we query the `reward` and `is_terminated` info from
     `env` and push them into `trajectory`.
  4. In the `PosEpisodeStage`, we push the `state` at the end of an episode and
     a dummy action into the `trajectory`.
  5. In the `PreEpisodeStage`, we pop out the latest `state` and `action` pair
     (which are dummy ones) from `trajectory`.

2. Update the inner `policy` given the context of `trajectory`, `env`, and
   `stage`.
  1. By default, we only `update!` the `policy` in the `PreActStage`. And it's
     dispatched to `update!(policy, trajectory, env, stage)`.
"""
function (agent::Agent)(stage::AbstractStage, env::AbstractEnv)
    update!(agent.trajectory, agent.policy, env, stage)
    update!(agent.policy, agent.trajectory, env, stage)
end

function (agent::Agent)(stage::PreActStage, env::AbstractEnv, action)
    update!(agent.trajectory, agent.policy, env, stage, action)
    update!(agent.policy, agent.trajectory, env, stage)
end

function RLBase.update!(
    ::AbstractPolicy,
    ::AbstractTrajectory,
    ::AbstractEnv,
    ::AbstractStage,
) end

#####
# Default behaviors for known trajectories
#####

function RLBase.update!(
    ::AbstractTrajectory,
    ::AbstractPolicy,
    ::AbstractEnv,
    ::AbstractStage,
) end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    ::AbstractPolicy,
    ::AbstractEnv,
    ::PreEpisodeStage,
)
    if length(trajectory) > 0
        pop!(trajectory[:state])
        pop!(trajectory[:action])
        if haskey(trajectory, :legal_actions_mask)
            pop!(trajectory[:legal_actions_mask])
        end
    end
end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PreActStage,
    action,
)
    s = policy isa NamedPolicy ? state(env, nameof(policy)) : state(env)
    push!(trajectory[:state], s)
    push!(trajectory[:action], action)
    if haskey(trajectory, :legal_actions_mask)
        lasm =
            policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) :
            legal_action_space_mask(env)
        push!(trajectory[:legal_actions_mask], lasm)
    end
end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PostEpisodeStage,
)
    # Note that for trajectories like `CircularArraySARTTrajectory`, data are
    # stored in a SARSA format, which means we still need to generate a dummy
    # action at the end of an episode. Here we simply select a random one using
    # the global rng. In theory it shouldn't affect the performance of specific
    # algorithm.
    # TODO: how to inject a local rng here to avoid polluting the global rng

    s = policy isa NamedPolicy ? state(env, nameof(policy)) : state(env)
    a = policy isa NamedPolicy ? rand(action_space(env, nameof(policy))) : rand(action_space(env))
    push!(trajectory[:state], s)
    push!(trajectory[:action], a)
    if haskey(trajectory, :legal_actions_mask)
        lasm =
            policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) :
            legal_action_space_mask(env)
        push!(trajectory[:legal_actions_mask], lasm)
    end
end

function RLBase.update!(
    trajectory::AbstractTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PostActStage,
)
    r = policy isa NamedPolicy ? reward(env, nameof(policy)) : reward(env)
    push!(trajectory[:reward], r)
    push!(trajectory[:terminal], is_terminated(env))
end

#####
# Pre-training
#####
function (agent::Agent)(stage::PreExperimentStage, env::AbstractEnv)
    update!(agent.policy, agent.trajectory, env, stage)
end