text:
using contact_sim

function main()
    #walker inits
    q1 = 0.2; u1 = -0.25;
    q2 = -0.4; u2 = 0.2;
    z0 = [q1 u1 q2 u2]
    steps = 4
    fps = 20
    walker = Walker(1.0, 0.5, 0.02, 1.0, 0.5, 1.0, 0.01, true, 1.0, 100.0, 0.28564, z0[3], z0[4], 0.0)
    z_init = [ 0.142819999999995  -0.326813112785275  -0.285640000000000   0.068243824367047]
    zpert = z_init + [0. 0.05 -0.1 0.2]

    #ball inits
    l = 0.4
    μ = 0.3
    ϵ = 0.5
    dt = 0.01
    m = 0.04
    g = [0., -9.81, 0.]
    R = (1.0/12.) * (2 * l * l)
    M = [[m 0. 0.];
        [0. m 0.];
        [0. 0. m * R]]
    Mi = [[1.0/m  0.  0.];
        [0. 1.0/m 0.];
        [0. 0. 1.0/(m*R)]]
    Δ = 0.001
    T = 150
    q₀ = [0.0, 0.2, π/180.0*30.]

    Torque = 27.2 #hinge motor torque (Nm)
    Fe = Torque*walker.l #kick force (N)
    th = 0.24 #kick angle (rad)
    v₀ = dt*Mi*[Fe*cos(th); Fe*sin(th); 0.0]
    println("Torque: ",Torque," vx: ", v₀[1], " vy: ", v₀[2])
    # v₀ = [4.2, 3.2, 0.0]

    # body = Object(l, μ, ϵ, m, R, M)
    # physics = Physics(dt, g, Δ, T)

    # zs, ts = simulate_walking(zpert,walker,steps)
    # qs, v = simulate(q₀, v₀, body, physics, get_ball_contact_jacobian)
    # render_kick_trajectory(ts, zs, qs, walker, steps, fps, body, physics, true)
end

main()


text:
struct TypedMessages end

struct Pinger
    pings::UInt
end

mutable struct Ponger
    pongs::UInt
end

struct Ping!{R <: Id}
    re::R
end

struct Pong!
    pongs::UInt
end

hear(s::Scene{Pinger}, msg::Ping!) = let pings = my(s).pings
    local pong

    for _ in 1:pings
        pong = ask(s, msg.re, Ping!(me(s)), Pong!)
    end

    say(s, stage(s), pong)
end

hear(s::Scene{Ponger}, msg::Ping!) = let pongs = (my(s).pongs += 1)
    say(s, msg.re, Pong!(pongs))
end

function hear(s::Scene{TypedMessages}, ::Genesis!)
    big = 100_000 # 1_000_000
    little = 220

    enter_pingpong = (pings) -> (enter!(s, Pinger(pings)), enter!(s, Ponger(0)))

    (pinger, ponger) = enter_pingpong(little)
    pong = ask(s, pinger, Ping!(ponger), Pong!)
    @test pong.pongs == little

    (pinger, ponger) = enter_pingpong(big)
    any_time = @elapsed ask(s, pinger, Ping!(ponger), Pong!)

    enter_pingpong =
        (pings) -> (enter!(s, Pinger(pings), Union{Leave!, Ping!, Pong!}),
                    enter!(s, Ponger(0), Union{Leave!, Ping!}))

    (pinger, ponger) = enter_pingpong(little)
    pong = ask(s, pinger, Ping!(ponger), Pong!)
    @test pong.pongs == little

    (pinger, ponger) = enter_pingpong(big)
    union_time = @elapsed ask(s, pinger, Ping!(ponger), Pong!)

    @debug "Typed message times" any_time union_time

    leave!(s)
end

@testset LuvvyTestSet "Typed Messages" begin
    play!(TypedMessages())
end

ext: jl


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

"""
    RegularGridDiscretizer(dims)

Discretize spatial region into regular grid of given
dimensions `dims`.
"""
struct RegularGridDiscretizer{N} <: AbstractDiscretizer
  dims::Dims{N}
end

discretize(r::Rectangle, d::RegularGridDiscretizer) =
  RegularGrid(extrema(r)..., dims=d.dims)


text:
module Remark

import Literate
import Documenter, DocumenterMarkdown
import DefaultApplication
using Random

export slideshow

const _pkg_assets = joinpath(dirname(@__DIR__), "assets")

const deps = [
    "http://gnab.github.io/remark/downloads/remark-latest.min.js",
    "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.9.0/katex.min.js",
    "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.9.0/contrib/auto-render.min.js",
    "https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.9.0/katex.min.css"
]

const depnames =  ["remark.min.js", "katex.min.js", "auto-render.min.js", "katex.min.css"]
const depfiles = joinpath.(_pkg_assets, depnames)

const styles_css = joinpath(_pkg_assets, "styles.css")

function slideshow(inputfile, outputdir = dirname(inputfile); documenter = true, css = styles_css)
    inputfile = realpath(abspath(inputfile))
    outputdir = realpath(abspath(outputdir))
    css = realpath(abspath(css))
    mkpath.(joinpath.(outputdir, ("src", "build")))
    mk_file = _create_index_md(inputfile, outputdir; documenter = documenter)
    _create_index_html(outputdir, mk_file)
    cp(css, joinpath(outputdir, "build", "styles.css"), force=true)
    rm(mk_file)
    return outputdir
end

function _create_index_md(inputfile, outputdir; documenter = true)
    if occursin(r".jl$", inputfile)
        Literate.markdown(inputfile, joinpath(outputdir, "src"), name = "index")
    else
        cp(inputfile, joinpath(outputdir, "src", "index.md"), force=true)
    end

    s1 = randstring('a':'z', 50)
    s2 = randstring('a':'z', 50)
    outputfile = joinpath(outputdir, "build", "index.md")
    if documenter
        _replace_line(
            joinpath(outputdir, "src", "index.md"),
            r"^(\s)*(--)(\s)*$" => s1,
            r"^(\s)*(\$\$)(\s)*$" => s2,
        )
        Documenter.makedocs(format = DocumenterMarkdown.Markdown(), root = outputdir)
        _replace_line(
            outputfile,
            Regex("^($s1)\$") => "--",
            s2 => raw"$$",
            r"^<a id=.*$" => "",
        )
    else
        cp(joinpath(outputdir, "src", "index.md"), outputfile, force=true)
    end
    outputfile
end

function _create_index_html(outputdir, md_file)

    Base.open(joinpath(outputdir, "build", "index.html"), "w") do f
        template = Base.open(joinpath(_pkg_assets, "indextemplate.html"))
        for line in eachline(template, keep=true)
            occursin(r"^(\s)*sfTiCgvZnilxkAh6ccwvfYSrKb4PmBKK", line) ? copytobuffer!(f, md_file) : write(f, line)
        end
        close(template)
    end
    for (name, file) in zip(depnames, depfiles)
        dest = joinpath(outputdir, "build", name)
        isfile(dest) || cp(file, dest)
    end
    dest = joinpath(outputdir, "build", "fonts")
    isdir(dest) || cp(joinpath(_pkg_assets, "fonts"), dest)
    joinpath(outputdir, "build", "index.html")
end

open(outputdir) =
    DefaultApplication.open(joinpath(outputdir, "build", "index.html"))

function copytobuffer!(f, filename)
    for line in eachline(filename, keep=true)
        write(f, line)
    end
end

function _replace_line(filename, pairs...)
    f = Base.open(filename)
    (tmp, tmpstream) = mktemp()
    for line in eachline(f, keep=false)
        newline = foldl(replace, pairs, init=line)
        write(tmpstream, newline)
        write(tmpstream, '\n')
    end
    close(f)
    close(tmpstream)
    mv(tmp, filename, force=true)
end

end # module


text:
# Uncomment the code to enable SearchLight support

#=
using SearchLight, SearchLight.QueryBuilder

Core.eval(SearchLight, :(config.db_config_settings = SearchLight.Configuration.load_db_connection()))

SearchLight.Loggers.empty_log_queue()

if SearchLight.config.db_config_settings["adapter"] != nothing
  SearchLight.Database.setup_adapter()
  SearchLight.Database.connect()
  SearchLight.load_resources()
end
=#


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ a0b3813e-adab-11eb-2983-616cf2bb6f5e
using LinearAlgebra, DifferentialEquations, PlutoUI,Plots

# ╔═╡ d84dbb08-b1a4-11eb-2d3b-0711fddd1347
html"""

Supplementary Materials

 Forecasting virus outbreakes with Neural ODEs.

"""

# ╔═╡ ef8d6690-720d-4772-a41f-b260d306b5b2
TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)

# ╔═╡ 30969341-9079-4732-bf55-d6bba2c2c16c
md"""
`sign(x)` returns 0 if x is 0, or ±1 for positive/negative x.
"""

# ╔═╡ 6139554e-c6c9-4252-9d64-042074f68391
begin
	f(y,a) =   sign(y)  + a+1 -  y
	f(y,a,t) = f(y,a) # just for the difeq solver
end

# ╔═╡ 61960e99-2932-4a8f-9e87-f64a7a043489
md"""
a = $(@bind a Slider(-5:.1:5, show_value=true, default=0) )

"""

# ╔═╡ 6f61180d-6900-48fa-998d-36110e79d2dc
gr()

# ╔═╡ 5027e1f8-8c50-4538-949e-6c95c550016e
md"""
### Flux argument

"""

# ╔═╡ 465f637c-0555-498b-a881-a2f6e5714cbb
md"""
y₀ = $(@bind y₀ Slider(-6:.1:6, show_value=true, default=2.0) )
"""

# ╔═╡ a09a0064-c6ab-4912-bc9f-96ab72b8bbca
sol = solve(  ODEProblem( f, y₀, (0,10.0), a ));

# ╔═╡ fb6e34b0-b4eb-11eb-3687-257bc18c3e52
md"""
### More

"""

# ╔═╡ bc96abb0-b4ea-11eb-21a1-fdd01945f9c6
md"""
### Auxiliary Code

"""

# ╔═╡ f8ee2373-6af0-4d81-98fb-23bde10198ef
function plotit(y₀, a)

	sol = solve(  ODEProblem( f, y₀, (0,10.0), a ));
	p = plot( sol , legend=false, background_color_inside=:black , ylims=(-7,7), lw=3, c=:red)

	# plot direction field:
		xs = Float64[]
		ys = Float64[]

	    lrx = LinRange( xlims(p)..., 30)
		for x in  lrx
			for y in LinRange( ylims(p)..., 30)
				v = [1, f(y, a ,x) ]
				v ./=  (100* (lrx[2]-lrx[1]))

				push!(xs, x - v[1], x + v[1], NaN)
				push!(ys, y - v[2], y + v[2], NaN)

			end
		end

	if a<1 hline!( [-1+a],c=:white,ls=:dash); end
	if a>-1  hline!( [1+a],c=:white,ls=:dash);end
	    hline!( [0 ],c=:pink,ls=:dash)
		plot!(xs, ys, alpha=0.7, c=:yellow)
	    ylabel!("y")
	    annotate!(-.5,y₀,text("y₀",color=:red))
	if a>-1
	   annotate!(5,1+a,text(round(1+a,digits=3),color=:white))
	end
	if a<1
	   annotate!(5,-1+a-.4,text(round(-1+a,digits=3),color=:white))
	end
    title!("Solution to y'=f(y,a)")
return(p)
end

# ╔═╡ a4686bca-90d6-4e02-961c-59f08fc37553
plotit(y₀, a)

# ╔═╡ Cell order:
# ╟─d84dbb08-b1a4-11eb-2d3b-0711fddd1347
# ╟─a0b3813e-adab-11eb-2983-616cf2bb6f5e
# ╟─ef8d6690-720d-4772-a41f-b260d306b5b2
# ╟─30969341-9079-4732-bf55-d6bba2c2c16c
# ╟─6139554e-c6c9-4252-9d64-042074f68391
# ╟─61960e99-2932-4a8f-9e87-f64a7a043489
# ╟─6f61180d-6900-48fa-998d-36110e79d2dc
# ╟─a09a0064-c6ab-4912-bc9f-96ab72b8bbca
# ╟─5027e1f8-8c50-4538-949e-6c95c550016e
# ╟─465f637c-0555-498b-a881-a2f6e5714cbb
# ╟─a4686bca-90d6-4e02-961c-59f08fc37553
# ╟─fb6e34b0-b4eb-11eb-3687-257bc18c3e52
# ╟─bc96abb0-b4ea-11eb-21a1-fdd01945f9c6
# ╟─f8ee2373-6af0-4d81-98fb-23bde10198ef


text:
# code reflection entry-points

export
    @code_lowered, @code_typed, @code_warntype,
    code_llvm, code_ptx, code_sass, @code_llvm, @code_ptx, @code_sass

#
# code_* replacements
#

# NOTE: default capability is a sane one for testing purposes

function code_llvm(io::IO, func::ANY, types::ANY=Tuple;
                   optimize::Bool=true, dump_module::Bool=false, cap::VersionNumber=v"2.0")
    mod, entry = irgen(func, types)
    if optimize
        optimize!(mod, cap)
    end
    if dump_module
        show(io, mod)
    else
        show(io, entry)
    end
end
code_llvm(func::ANY, types::ANY=Tuple; kwargs...) = code_llvm(STDOUT, func, types; kwargs...)

"""
    code_ptx([io], f, types, cap=v"2.0")

Prints the PTX assembly generated for the method matching the given generic function and
type signature to `io` which defaults to `STDOUT`.
"""
function code_ptx(io::IO, func::ANY, types::ANY=Tuple;
                  cap::VersionNumber=v"2.0", kernel::Bool=false)
    @assert isa(func, Core.Function)
    tt = Base.to_tuple_type(types)
    kernel && check_kernel(func, tt)

    ptx,_ = compile_function(func, tt, cap; kernel=kernel)
    # TODO: this code contains all the functions in the call chain,
    #       is it possible to implement `dump_module`?
    print(io, ptx)
end
code_ptx(func::ANY, types::ANY=Tuple; kwargs...) = code_ptx(STDOUT, func, types; kwargs...)

"""
    code_sass([io], f, types, cap=v"2.0")

Prints the SASS code generated for the method matching the given generic function and type
signature to `io` which defaults to `STDOUT`.

Note that the method needs to be a valid entry-point kernel, ie. it should not return any
values.
"""
function code_sass(io::IO, func::ANY, types::ANY=Tuple; cap::VersionNumber=v"2.0")
    @assert isa(func, Core.Function)
    tt = Base.to_tuple_type(types)
    check_kernel(func, tt)

    ptx,_ = compile_function(func, tt, cap)

    fn = tempname()
    gpu = "sm_$(cap.major)$(cap.minor)"
    # NOTE: this might not match what is being executed, due to the PTX->SASS conversion
    #       by the driver possibly not matching what `ptxas` (part of the toolkit) does.
    # TODO: see how `nvvp` extracts SASS code when doing PC sampling, and copy that.
    Base.run(`ptxas --gpu-name $gpu --output-file $fn --input-as-string $ptx`)
    try
        print(io, readstring(`cuobjdump --dump-sass $fn`))
    finally
        rm(fn)
    end
end
code_sass(func::ANY, types::ANY=Tuple; kwargs...) = code_sass(STDOUT, func, types; kwargs...)

#
# @code_* replacements
#

for fname in [:code_lowered, :code_typed, :code_warntype, :code_llvm, :code_ptx, :code_sass]
    # types often need to be converted (eg. CuArray -> CuDeviceArray),
    # so generate a type-converting wrapper, and a macro to call it
    fname_wrapper = Symbol(fname, :_cputyped)

    @eval begin
        function $fname_wrapper(func, types)
            _, codegen_types, _ =
                convert_arguments(fill(Symbol(), length(types.parameters)),
                                  types.parameters)
            $fname(func, codegen_types)
        end

        macro ($fname)(ex0)
            if ex0.head == :macrocall
                # @cuda (...) f()
                ex0 = ex0.args[3]
            end
            Base.gen_call_with_extracted_types($(Expr(:quote,fname_wrapper)), ex0)
        end
    end
end


text:
using Useless # let's say our project depends on this useless package: https://github.com/crstnbr/Useless.jl

include("somedefinitions.jl")

println(life_universe_everything("What is the meaning of life?"))

# Done.

max_stars_repo_path: 2 Tips and tricks/3_share-and-instantiate-projects/my-project/main.jl

max_stars_repo_name: crstnbr/julia-workshop


text:
@system C3Base begin
    Ci: intercellular_co2 ~ hold
    I2: effective_irradiance ~ hold
    T: leaf_temperature ~ hold

    Tk(T): absolute_leaf_temperature ~ track(u"K")

    # Arrhenius equation
    Tb: base_temperature => 25 ~ preserve(u"°C", parameter)
    Tbk(Tb): absolute_base_temperature ~ preserve(u"K")
    T_dep(T, Tk, Tb, Tbk; Ea(u"kJ/mol")): temperature_dependence_rate => begin
        exp(Ea * (T - Tb) / (Tbk * u"R" * Tk))
    end ~ call
end

@system C3c(C3Base) begin
    # Michaelis constant of rubisco for CO2 of C3 plants, ubar, from Bernacchi et al. (2001)
    Kc25: rubisco_constant_for_co2_at_25 => 404.9 ~ preserve(u"μbar", parameter)
    # Activation energy for Kc, Bernacchi (2001)
    Eac: activation_energy_for_co2 => 79.43 ~ preserve(u"kJ/mol", parameter)
    Kc(T_dep, Kc25, Eac): rubisco_constant_for_co2 => begin
        Kc25 * T_dep(Eac)
    end ~ track(u"μbar")

    # Michaelis constant of rubisco for O2, mbar, from Bernacchi et al., (2001)
    Ko25: rubisco_constant_for_o2_at_25 => 278.4 ~ preserve(u"mbar", parameter)
    # Activation energy for Ko, Bernacchi (2001)
    Eao: activation_energy_for_o2 => 36.38 ~ preserve(u"kJ/mol", parameter)
    Ko(T_dep, Ko25, Eao): rubisco_constant_for_o2 => begin
        Ko25 * T_dep(Eao)
    end ~ track(u"mbar")

    # mesophyll O2 partial pressure
    Om: mesophyll_o2_partial_pressure => 210 ~ preserve(u"mbar", parameter)
    # effective M-M constant for Kc in the presence of O2
    Km(Kc, Om, Ko): rubisco_constant_for_co2_with_o2 => begin
        Kc * (1 + Om / Ko)
    end ~ track(u"μbar")

    Vcm25: maximum_carboxylation_rate_at_25 => 108.4 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    EaVc: activation_energy_for_carboxylation => 52.1573 ~ preserve(u"kJ/mol", parameter)
    Vcmax(T_dep, Vcm25, EaVc): maximum_carboxylation_rate => begin
        Vcm25 * T_dep(EaVc)
    end ~ track(u"μmol/m^2/s" #= CO2 =#)
end

@system C3j(C3Base) begin
    Jm25: maximum_electron_transport_rate_at_25 => 169.0 ~ preserve(u"μmol/m^2/s" #= Electron =#, parameter)
    Eaj: activation_energy_for_electron_transport => 23.9976 ~ preserve(u"kJ/mol", parameter)
    Sj: electron_transport_temperature_response => 616.4 ~ preserve(u"J/mol/K", parameter)
    Hj: electron_transport_curvature => 200 ~ preserve(u"kJ/mol", parameter)
    Jmax(Tk, Tbk, T_dep, Jm25, Eaj, Sj, Hj): maximum_electron_transport_rate => begin
        R = u"R"
        Jm25 * begin
            T_dep(Eaj) *
            (1 + exp((Sj*Tbk - Hj) / (R*Tbk))) /
            (1 + exp((Sj*Tk  - Hj) / (R*Tk)))
        end
    end ~ track(u"μmol/m^2/s" #= Electron =#)

    # θ: sharpness of transition from light limitation to light saturation
    θ: light_transition_sharpness => 0.7 ~ preserve(parameter)
    J(I2, Jmax, θ): electron_transport_rate => begin
        a = θ
        b = -(I2+Jmax)
        c = I2*Jmax
        a*J^2 + b*J + c ⩵ 0
    end ~ solve(lower=0, upper=Jmax, u"μmol/m^2/s")
end

@system C3p(C3Base) begin
    Tp25: triose_phosphate_limitation_at_25 => 16.03 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    EaTp: activation_energy_for_Tp => 47.10 ~ preserve(u"kJ/mol", parameter)
    Tp(T_dep, Tp25, EaTp): triose_phosphate_limitation => begin
        Tp25 * T_dep(EaTp)
    end ~ track(u"μmol/m^2/s" #= CO2 =#)
end

@system C3r(C3Base) begin
    Rd25: dark_respiration_at_25 => 1.08 ~ preserve(u"μmol/m^2/s" #= CO2 =#, parameter)
    Ear: activation_energy_for_respiration => 49.39 ~ preserve(u"kJ/mol", parameter)
    Rd(T_dep, Rd25, Ear): dark_respiration => begin
        Rd25 * T_dep(Ear)
    end ~ track(u"μmol/m^2/s")
    #Rm(Rd) => 0.5Rd ~ track(u"μmol/m^2/s")

    # CO2 compensation point in the absence of day respiration, value from Bernacchi (2001)
    Γ25: co2_compensation_point_at_25 => 42.75 ~ preserve(u"μbar", parameter)
    Eag: activation_energy_for_co2_compensation_point => 37.83 ~ preserve(u"kJ/mol", parameter)
    Γ(T_dep, Γ25, Eag): co2_compensation_point => begin
        Γ25 * T_dep(Eag)
    end ~ track(u"μbar")
end

@system C3Rate(C3c, C3j, C3p, C3r) begin
    Ac(Vcmax, Ci, Γ, Km, Rd): enzyme_limited_photosynthesis_rate => begin
        Vcmax * (Ci - Γ) / (Ci + Km) - Rd
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    # light and electron transport limited A mediated by J
    Aj(J, Ci, Γ, Rd): transport_limited_photosynthesis_rate => begin
        J * (Ci - Γ) / 4(Ci + 2Γ) - Rd
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    Ap(Tp, Rd): triose_phosphate_limited_photosynthesis_rate => begin
        3Tp - Rd
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    A_net(Ac, Aj, Ap): net_photosynthesis => begin
        min(Ac, Aj, Ap)
    end ~ track(u"μmol/m^2/s" #= CO2 =#)

    A_gross(A_net, Rd): gross_photosynthesis => begin
        # gets negative when PFD = 0, Rd needs to be examined, 10/25/04, SK
        A_gross = A_net + Rd
        #max(A_gross, zero(A_gross))
    end ~ track(u"μmol/m^2/s" #= CO2 =#)
end

@system C3(C3Rate)


text:
struct OrnsteinUhlenbeck{T1,T2,T3}
  Θ::T1
  μ::T2
  σ::T3
end
# http://www.math.ku.dk/~susanne/StatDiff/Overheads1b.pdf
function (p::OrnsteinUhlenbeck)(W,dt,rng) #dist
  if typeof(W.dW) <: AbstractArray
    rand_val = wiener_randn(rng,W.dW)
  else
    rand_val = wiener_randn(rng,typeof(W.dW))
  end
  drift = p.μ .+ (W[end] .- p.μ) .* exp.(-p.Θ*dt)
  diffusion = p.σ .* sqrt.((1 .- exp.(-2p.Θ*dt))./(2p.Θ))
  drift .+ rand_val .* diffusion .- W[end]
end

#=
http://www.tandfonline.com/doi/pdf/10.1080/14697688.2014.941913?needAccess=true
http://www.tandfonline.com/doi/full/10.1080/14697688.2014.941913?src=recsys
=#

#=
(qX + r) = Θ(μ-x) = Θμ - Θx
q = -Θ
r = Θμ
https://arxiv.org/pdf/1011.0067.pdf page 18
=#
function ou_bridge(ou,W,W0,Wh,q,h) end
function ou_bridge!(rand_vec,ou,W,W0,Wh,q,h) end

function OrnsteinUhlenbeckProcess(Θ,μ,σ,t0,W0,Z0=nothing;kwargs...)
  ou = OrnsteinUhlenbeck(Θ,μ,σ)
  NoiseProcess(t0,W0,Z0,ou,nothing;kwargs...)
end

struct OrnsteinUhlenbeck!{T1,T2,T3}
  Θ::T1
  μ::T2
  σ::T3
end

function (p::OrnsteinUhlenbeck!)(rand_vec,W,dt,rng) #dist!
  wiener_randn!(rng,rand_vec)
  @.. rand_vec = p.μ+(W[end]-p.μ)*exp(-p.Θ*dt) + rand_vec*p.σ*sqrt((1-exp.(-2*p.Θ.*dt))/(2*p.Θ)) - W[end]
end
function OrnsteinUhlenbeckProcess!(Θ,μ,σ,t0,W0,Z0=nothing;kwargs...)
  ou = OrnsteinUhlenbeck!(Θ,μ,σ)
  NoiseProcess(t0,W0,Z0,ou,nothing;kwargs...)
end


text:
import Clustering
import SparseArrays

"""
    Mapper result type
"""
struct Mapper
    adj::AbstractMatrix{<:Integer}
    filter::Vector{<:Real}
    patches::Vector{Vector{<:Integer}}
    centers::Matrix{<:Real}
end

Base.show(io::IO, mpr::Mapper) = show(io, "Mapper[$(length(mpr.patches))]")

"""
    mapper(X::AbstractMatrix{<:Real}[; filter = eccentricity])

Calculate a simplicial complex of `X` dataset using Mapper algorithm.
"""
function mapper(X::AbstractMatrix{<:Real}; kwargs...)

    # setup parameters
    filterfn = eccentricity
    coverfn = balanced
    clusterfn = Clustering.kmeans
    clusterselectionfn = silhouette
    for (p,v) in kwargs
        if p == :filter
            @assert isa(v, Function) "`$p` parameter must be function"
            filterfn = v
        elseif p == :cover
            @assert isa(v, Function) "`$p` parameter must be function"
            coverfn = v
        elseif p == :clustering
            @assert isa(v, Function) "`$p` parameter must be function"
            clusterfn = v
        elseif p == :clustselection
            @assert isa(v, Function) "`$p` parameter must be function"
            clusterselectionfn = v
        elseif p == :seed
            Random.seed!(v)
        end
    end

    # calculate filter range
    flt = filterfn(X)

    # construct cover of the filter range
    covering = coverfn(flt; kwargs...)

    # using clustering algorithm create patches from cover elements
    patches = Vector{Int}[]
    for c in covering
        println(c)
        if length(c) == 1
            push!(patches, c)
        end
        lbls = clusterselectionfn(clusterfn, view(X, :, c); kwargs...)
        plen = length(unique(lbls))
        println(plen)
        for i in 1:plen
            push!(patches, c[findall(isequal(i), lbls)])
        end
    end

    # combine all patches & determine which has overlaps
    P = length(patches)
    adj = SparseArrays.spzeros(UInt8, P, P)
    for i in 1:P
        println("$i => ", patches[i])
        for j in i+1:P
            overlap = intersect(patches[i], patches[j])
            if length(overlap) > 0
                adj[i,j] = 0x01
            end
        end
    end

    # calculate centers of cover patches
    cntrs = hcat((mean(view(X, : ,p), dims=2) for p in patches)...)

    return Mapper(adj, flt, patches, cntrs)
end

@recipe function f(mpr::Mapper; complex_layout = circular_layout,
                   minvsize = 15, maxvsize = 35)

    xpos, ypos = complex_layout(mpr)

    # set image limits
    xlims --> extrema(xpos) .* 1.2
    ylims --> extrema(ypos) .* 1.2

    # show 1-skeleton
    for i in 1:size(mpr.adj, 1)
        idxs = findall(e->e>0, view(mpr.adj, i, :))
        push!(idxs, i)
        @series begin
            seriestype := :path
            linewidth --> 2
            linecolor --> :black
            label --> ""
            xpos[idxs], ypos[idxs]
        end
    end

    # calculate vertex attribues
    n = length(mpr.patches)
    xcrd = zeros(n)
    ycrd = zeros(n)
    zcol = zeros(n)
    msize = fill(1,n)
    for (i, p) in enumerate(mpr.patches)
        zcol[i] = mean(mpr.filter[p])
        msize[i] = length(p)
        xcrd[i] = xpos[i]
        ycrd[i] = ypos[i]
    end
    manno = map(string, msize)
    smin, smax = extrema(msize)
    srng = smax-smin
    msize = (maxvsize-minvsize).*(msize .- smin)./srng .+ minvsize

    # show nodes
    @series begin
        seriestype := :scatter
        markersize := msize
        label --> ""
        zcolor := zcol
        series_annotations := manno
        xcrd, ycrd
    end
end


text:
module CreateTableHeres

import SearchLight.Migrations: create_table, column, primary_key, add_index, drop_table
using SearchLight

function up()
  SearchLight.query("""
    CREATE TABLE heres (
      id serial,
      name varchar(255),
      lnglat geometry(point, 4326),
      primary key(id)
    );
  """)
  # create_table(:heres) do
  #   [
  #     primary_key()
  #     column(:name, :string)
  #     column(:lnglat, :geometry)
  #   ]
  # end

  # add_index(:heres, :name)
end

function down()
  drop_table(:heres)
end

end

max_stars_repo_name: huaxk/rlservice

max_issues_repo_name: huaxk/rlservice

max_forks_repo_name: huaxk/rlservice


text:
# Helper data for CIE observer functions
include("cie_data.jl")

# readonly static 3x3 matrix
struct Mat3x3{T} <: AbstractMatrix{T}
    e::NTuple{9, T}
end
Mat3x3(mat::Matrix{T}) where {T} = Mat3x3{T}(Tuple(mat))
Base.IndexStyle(::Type{<:Mat3x3}) = IndexLinear()
Base.size(::Mat3x3) = (3, 3)
Base.getindex(M::Mat3x3{T}, i::Int) where {T} = M.e[i]

macro mul3x3(C, M, c1, c2, c3)
    esc(quote
        F = typeof(0.5f0 * $c1) === Float32 ? Float32 : eltype($M)
        if $c1 isa N0f8 && $c2 isa N0f8 && $c3 isa N0f8
            s = 0xffffff
            c1x = Int32(reinterpret($c1)) * 0x10101
            c2x = Int32(reinterpret($c2)) * 0x10101
            c3x = Int32(reinterpret($c3)) * 0x10101
        else
            s = true
            c1x, c2x, c3x = $c1, $c2, $c3
        end
        @inbounds ret = $C(
            muladd(F($M[1,1] / s), c1x, muladd(F($M[1,2] / s), c2x, F($M[1,3] / s) * c3x)),
            muladd(F($M[2,1] / s), c1x, muladd(F($M[2,2] / s), c2x, F($M[2,3] / s) * c3x)),
            muladd(F($M[3,1] / s), c1x, muladd(F($M[3,2] / s), c2x, F($M[3,3] / s) * c3x)))
        ret
    end)
end

# for optimization
div60(x) = x / 60
_div60(x::T) where T = muladd(x, T(1/960), x * T(0x1p-6))
if reduce(max, _div60.((90.0f0,))) == 1.5f0
    div60(x::T) where T <: Union{Float32, Float64} = _div60(x)
else
    # force two-step multiplication
    div60(x::T) where T <: Union{Float32, Float64} = x * T(0x1p-6) + x * T(1/960)
end

# mod6 supports the input `x` in [-2^28, 2^29]
mod6(::Type{T}, x::Int32) where T = unsafe_trunc(T, x - 6 * ((widemul(x, 0x2aaaaaaa) + Int64(0x20000000)) >> 0x20))

# Approximation of the reciprocal of the cube root, x^(-1/3).
# assuming that x > 0.003, the conditional branches are omitted.
@inline function rcbrt(x::Float64)
    ix = reinterpret(UInt64, x)
    e0 = (ix >> 0x34) % UInt32
    ed = e0 ÷ 0x3
    er = e0 - ed * 0x3
    a = 0x000b_f2d7 - 0x0005_5718 * er
    e = (UInt32(1363) - ed) << 0x14 | a
    t1 = reinterpret(Float64, UInt64(e) << 0x20)
    h1 = muladd(t1^2, -x * t1, 1.0)
    t2 = muladd(@evalpoly(h1, 1/3, 2/9, 14/81), h1 * t1, t1)
    h2 = muladd(t2^2, -x * t2, 1.0)
    t3 = muladd(muladd(2/9, h2, 1/3), h2 * t2, t2)
    reinterpret(Float64, reinterpret(UInt64, t3) & 0xffff_ffff_8000_0000)
end
@inline function rcbrt(x::Float32)
    ix = reinterpret(UInt32, x)
    e0 = ix >> 0x17 + 0x2
    ed = e0 ÷ 0x3
    er = e0 - ed * 0x3
    a = 0x005f_9cbe - 0x002a_bd7d * er
    t1 = reinterpret(Float32, (UInt32(169) - ed) << 0x17 | a)
    h1 = muladd(t1^2, -x * t1, 1.0f0)
    t2 = muladd(muladd(2/9f0, h1, 1/3f0), h1 * t1, t1)
    h2 = muladd(t2^2, -x * t2, 1.0f0)
    t3 = muladd(1/3f0, h2 * t2, t2)
    reinterpret(Float32, reinterpret(UInt32, t3) & 0xffff_f000)
end

cbrt01(x) = cbrt(x)
@inline function cbrt01(x::Float64)
    r = rcbrt(x) # x^(-1/3)
    h = muladd(r^2, -x * r, 1.0)
    e = muladd(2/9, h, 1/3) * h * r
    muladd(r, x * r, x * e * (r + r + e)) # x * x^(-2/3)
end
@inline function cbrt01(x::Float32)
    r = Float64(rcbrt(x)) # x^(-1/3)
    h = muladd(r^2, -Float64(x) * r, 1.0)
    e = muladd(muladd(14/81, h, 2/9), h, 1/3) * h
    Float32(1 / muladd(r, e, r))
end

pow3_4(x) = (y = @fastmath(sqrt(x)); y*@fastmath(sqrt(y))) # x^(3/4)

# `pow5_12` is called from `srgb_compand`.
pow5_12(x) = pow3_4(x) / cbrt(x) # 5/12 == 1/2 + 1/4 - 1/3 == 3/4 - 1/3
@inline function pow5_12(x::Float64)
    p3_4 = pow3_4(x)
    # x^(-1/6)
    if x < 0.02
        t0 = @evalpoly(x, 3.1366722556806232,
            -221.51395962221136, 19788.889459114234, -905934.6541469148, 1.5928561711645417e7)
    elseif x < 0.12
        t0 = @evalpoly(x, 2.3135905865468644,
            -26.43664640894651, 385.0146581045545, -2890.0920682466267, 8366.343115590817)
    elseif x < 1.2
        t0 = @evalpoly(x, 1.7047813285940905, -3.1261253501167308,
            7.498744828350077, -10.100319516746419, 6.820601476522508, -1.7978894213531524)
    else
        return p3_4 / cbrt01(x)
    end
    # x^(-1/3)
    t1 = t0 * t0
    h1 = muladd(t1^2, -x * t1, 1.0)
    t2 = muladd(h1, 1/3 * t1, t1)
    h2 = muladd(t2^2, -x * t2, 1.0)
    t2h = @evalpoly(h2, 1/3, 2/9, 14/81) * h2 * t2 # Taylor series of (1-h)^(-1/3)
    # x^(3/4) * x^(-1/3)
    muladd(p3_4, t2, p3_4 * t2h)
end
@inline function pow5_12(x::Float32)
    # x^(-1/3)
    rc = rcbrt(x)
    rcx = -rc * x
    rch = muladd(muladd(rc, x, rcx), -rc^2, muladd(rc^2, rcx, 1.0f0)) # 1 - x * rc^3
    rce = muladd(2/9f0, rch, 1/3f0) * rch * rc
    # x^(3/4)
    p3_4_f64 = pow3_4(Float64(x))
    p3_4r = reinterpret(Float64, reinterpret(UInt64, p3_4_f64) & 0xffffffff_e0000000)
    p3_4 = Float32(p3_4r)
    p3_4e = Float32(p3_4_f64 - p3_4r)
    # x^(3/4) * x^(-1/3)
    muladd(p3_4, rc, muladd(p3_4, rce, p3_4e * rc))
end

# `pow12_5` is called from `invert_srgb_compand`.
pow12_5(x) = pow12_5(Float64(x))
pow12_5(x::BigFloat) = x^big"2.4"
@inline function pow12_5(x::Float64)
    # x^0.4
    t1 = @evalpoly(@fastmath(min(x, 1.75)), 0.24295462640373672,
        1.7489099720303518, -1.9919942887850166, 1.3197188815160004, -0.3257258790067756)
    t2 = muladd(2/5, muladd(x / t1^2, @fastmath(sqrt(t1)), -t1), t1) # Newton's method
    t3 = muladd(2/5, muladd(x / t2^2, @fastmath(sqrt(t2)), -t2), t2)
    t4 = muladd(2/5, muladd(x / t3^2, @fastmath(sqrt(t3)), -t3), t3)
    # x^0.4 * x^2
    rx = reinterpret(Float64, reinterpret(UInt64, x) & 0xffffffff_f8000000) # hi
    e = x - rx # lo
    muladd(t4, rx^2, t4 * (rx + rx + e) * e)
end

pow7(x) = (y = x*x*x; y*y*x)

# TODO: migrate to ColorTypes.jl
@inline function atan360_kernel(t::Float64)
    @evalpoly(t^2,
        0.8952465548919112,  -0.2984155182972285,   0.1790493109637673,   -0.12789236387151418,
        0.09947179554077099, -0.08138502549008089,  0.06884985860325084,  -0.059532134551367105,
        0.05164982834653973, -0.042514874197367984, 0.028600080170923112, -0.011001809246782802)
end
@inline function atan360_kernel(t::Float32)
    @evalpoly(t^2,
        0.89524657f0, -0.2984143f0, 0.17899014f0, -0.12683357f0, 0.090689056f0, -0.045563098f0)
end

# a variant of `atand` returning the angle in the range of [0, 360]
atan360(y, x) = (a = atand(y, x); signbit(a) ? oftype(a, a + 360) : a)
@inline function atan360(y::T, x::T) where T <: Union{Float32, Float64}
    (isnan(x) | isnan(y)) && return T(NaN)
    ax, ay = abs(x), abs(y)
    n, m = @fastmath minmax(ax, ay)
    if m == T(Inf)
        d0 = n == T(Inf) ? T(45) : T(0)
    else
        m = m == T(0) ? T(0.5) : m
        ta = (n + n) > m ? T(0.5) : T(0) # 1-step CORDIC
        # ro=(n + n) > m ? T(atand(0.5) / 64) : T(0)
        ro = @fastmath max(T(0), ta - T(0.5 - 0.4150789246418436))
        n1 = n - ta * m
        m1 = m + ta * n
        t = n1 / m1 # in [0, 0.5]
        p = atan360_kernel(t)
        d0 = muladd(t, p, ro) * T(64)
    end
    b1 = T( 90) + flipsign(T( -90), x)
    b2 = T(180) + flipsign(T(-180), y)
    d1 = ay > ax ? T(90) - d0 : d0
    d2 = b1 + flipsign(d1, x) # signbit(x) ? T(180) - d1 : d1
    dd = b2 + flipsign(d2, y) # signbit(y) ? T(360) - d2 : d2
    return dd
end

# override only the `Lab` and `Luv` versions just for now
@inline ColorTypes.hue(c::Lab) = atan360(c.b, c.a)
@inline ColorTypes.hue(c::Luv) = atan360(c.v, c.u)

@inline function sin_kernel(x::Float64)
    x * @evalpoly(x^2,
        1.117010721276371, -0.23228479064016105, 0.014491237085286733, -0.00043049771889962576,
        7.460244157055791e-6, -8.462038405688494e-8, 6.767827147797153e-10, -3.987482394639226e-12)
end
@inline function sin_kernel(x::Float32)
    y = @evalpoly(x^2, 0.11701072f0, -0.23228478f0, 0.014491233f0, -0.0004304645f0, 7.368049f-6)
    muladd(x, y, x)
end
@inline function cos_kernel(x::Float64)
    @evalpoly(x^2, 1.0,
        -0.6238564757231793, 0.06486615038362423, -0.002697811198135598, 6.010882091788964e-5,
        -8.333171603045294e-7, 7.876495580576226e-9, -5.351642293798961e-11)
end
@inline function cos_kernel(x::Float32)
    @evalpoly(x^2, 1.0f0, -0.6238565f0,0.06486606f0,-0.002697325f0, 5.904168f-5)
end

sincos360(x) = sincos(deg2rad(x))
@inline function sincos360(x::T) where T <: Union{Float32, Float64}
    isfinite(x) || return T(NaN), T(NaN)
    t = x - round(x * T(1/360)) * 360 # [-180, 180]
    a0 = @fastmath abs(t)              # [0, 180]
    a1 = a0 <= 90 ? a0 : 180 - a0      # [0, 90]
    a2 = a1 <= 45 ? a1 : 90 - a1       # [0, 45]
    ax = a2 * T(1/64)
    sn, cs  = sin_kernel(ax), cos_kernel(ax)
    sn1, cs1 = a1 === a2 ? (sn, cs) : (cs, sn)
    return @fastmath flipsign(sn1, t), flipsign(cs1, 90 - a0)
end

"""
    x, y = polar_to_cartesian(r, theta)

Convert a polar coordinate of radius `r` and angle `theta` (in degrees) to a
Cartesian coordinate.
"""
@inline function polar_to_cartesian(r, theta)
    y, x = r .* sincos360(theta)
    return x, y
end

# Linear interpolation in [a, b] where x is in [0,1],
# or coerced to be if not.
function lerp(x, a, b)
    a + (b - a) * max(min(x, one(x)), zero(x))
end

clamp01(v::T) where {T<:Fractional} = ifelse(v < zero(T), zero(T), ifelse(v > oneunit(T), oneunit(T), v))
clamp01(v::T) where {T<:Union{Bool,N0f8,N0f16,N0f32,N0f64}} = v

"""
    HexNotation{C, A, N}

This is a private type for specifying the style of hex notations. It is not
recommended to use this type and its derived types in user scripts or other
packages, since they may change in the future without notice.

# Arguments
- `C`: a base colorant type.
- `A`: a symbol (`:upper` or `:lower`) to specify the letter casing.
- `N`: a total number of digits.
"""
abstract type HexNotation{C, A, N} end
abstract type HexAuto <: HexNotation{Colorant,:upper,0} end
abstract type HexShort{A} <: HexNotation{Colorant,A,0} end

"""
    hex(c::Colorant)
    hex(c::Colorant, style::Symbol)

Convert a color to a hexadecimal string, optionally specifying its style.

# Arguments
- `c`: a target color.
- `style`: a symbol to specify the hexadecimal notation. Spesifying the
  uppercase symbols means the return values are in uppercase. The following
  symbols are available:
  - `:AUTO`: notation automatically selected according to the type of `c`
  - `:RRGGBB`/`:rrggbb`: 6-digit opaque notation
  - `:AARRGGBB`/`:aarrggbb`: 8-digit notation with alpha at the head
  - `:RRGGBBAA`/`:rrggbbaa`: 8-digit notation with alpha at the tail
  - `:RGB`/`:rgb`/`:ARGB`/`:argb`/`:RGBA`/`:rgba`: 3-digit or 4-digit noatation
  - `:S`/`:s`: short notation if available

# Examples
```jldoctest; setup = :(using Colors)
julia> hex(RGB(1,0.5,0))
"FF8000"

julia> hex(ARGB(1,0.5,0,0.25))
"40FF8000"

julia> hex(HSV(30,1.0,1.0), :AARRGGBB)
"FFFF8000"

julia> hex(ARGB(1,0.533,0,0.267), :rrggbbaa)
"ff880044"

julia> hex(ARGB(1,0.533,0,0.267), :rgba)
"f804"

julia> hex(ARGB(1,0.533,0,0.267), :S)
"4F80"
```

!!! compat "Colors v0.12"
    `style` requires at least Colors v0.12.
"""
hex(c::Colorant) = _hex(HexAuto, c) # there is no need to search the dictionary
hex(c::Colorant, style::Symbol) = _hex(get(_hex_styles, style, HexAuto), c)

const _hex_styles = Dict{Symbol, Type}(
    :AUTO => HexAuto,
    :S => HexShort{:upper}, :s => HexShort{:lower},
    :RGB => HexNotation{RGB,:upper,3}, :rgb => HexNotation{RGB,:lower,3},
    :ARGB => HexNotation{ARGB,:upper,4}, :argb => HexNotation{ARGB,:lower,4},
    :RGBA => HexNotation{RGBA,:upper,4}, :rgba => HexNotation{RGBA,:lower,4},
    :RRGGBB => HexNotation{RGB,:upper,6}, :rrggbb => HexNotation{RGB,:lower,6},
    :AARRGGBB => HexNotation{ARGB,:upper,8}, :aarrggbb => HexNotation{ARGB,:lower,8},
    :RRGGBBAA => HexNotation{RGBA,:upper,8}, :rrggbbaa => HexNotation{RGBA,:lower,8},
)
@inline function _hexstring(::Type{T}, u::U, itr) where {C, T <: HexNotation{C,:upper}, U <: Unsigned}
    s = UInt8(8sizeof(u) - 4)
    @inbounds String([b"0123456789ABCDEF"[((u << i) >> s) + 1] for i in itr])
end
@inline function _hexstring(::Type{T}, u::U, itr) where {C, T <: HexNotation{C,:lower}, U <: Unsigned}
    s = UInt8(8sizeof(u) - 4)
    @inbounds String([b"0123456789abcdef"[((u << i) >> s) + 1] for i in itr])
end

_to_uint32(c::Colorant) = reinterpret(UInt32, ARGB32(c))
_to_uint32(c::TransparentColor) = reinterpret(UInt32, alphacolor(RGB24(c), clamp01(alpha(c))))
_to_uint32(c::C) where C <: Union{AbstractRGB, TransparentRGB} =
    reinterpret(UInt32, ARGB32(correct_gamut(c)))
_to_uint32(c::C) where C <: Union{AbstractGray, TransparentGray} =
    reinterpret(UInt32, AGray32(clamp01(gray(c)), clamp01(alpha(c))))

_hex(t::Type, c::Colorant) = _hex(t, _to_uint32(c))

_hex(::Type{HexAuto}, c::Color) = _hex(HexNotation{RGB,:upper,6}, c)
_hex(::Type{HexAuto}, c::AlphaColor) = _hex(HexNotation{ARGB,:upper,8}, c)
_hex(::Type{HexAuto}, c::ColorAlpha) = _hex(HexNotation{RGBA,:upper,8}, c)

function _hex(::Type{HexShort{A}}, c::Colorant) where A
    u = _to_uint32(c)
    s = u == (u & 0x0F0F0F0F) * 0x11
    c isa AlphaColor && return _hex(HexNotation{ARGB, A, s ? 4 : 8}, u)
    c isa ColorAlpha && return _hex(HexNotation{RGBA, A, s ? 4 : 8}, u)
    _hex(HexNotation{RGB, A, s ? 3 : 6}, u)
end

# for 3-digit or 4-digit notations
function _hex(t::Type{T}, u::UInt32) where {C <:Union{RGB, ARGB, RGBA}, A, T <: HexNotation{C,A}}
    # To double the number of digits, we multiply each element by 17 (= 0x11).
    # Thus, we divide each element by 17 here, to halve the number of digits.
    u64 = UInt64(u)
    # TODO: use SIMD `move` with zero extension (e.g. vpmovzxbw)
    unpacked = ((u64 & 0xFF00FF00)<<24) | (u64 & 0x00FF00FF) # 0x00AA00GG00RR00BB
    # `all(x -> round(x / 17) == (x * 15 + 135) >> 8, 0:255) == true`
    q = muladd(unpacked, 0xF,  0x0087_0087_0087_0087) # 0x0Aaa0Ggg0Rrr0Bbb
    t <: HexNotation{ARGB} && return _hexstring(t, q, (0x04, 0x24, 0x14, 0x34))
    t <: HexNotation{RGBA} && return _hexstring(t, q, (0x24, 0x14, 0x34, 0x04))
    _hexstring(t, q, (0x24, 0x14, 0x34))
end

# for 6-digit or 8-digit notations
_hex(t::Type{HexNotation{ RGB,A,6}}, u::UInt32) where {A} = _hexstring(t, u, 0x8:0x4:0x1C)
_hex(t::Type{HexNotation{ARGB,A,8}}, u::UInt32) where {A} = _hexstring(t, u, 0x0:0x4:0x1C)
_hex(t::Type{HexNotation{RGBA,A,8}}, u::UInt32) where {A} =
    _hexstring(t, u, (0x8, 0xC, 0x10, 0x14, 0x18, 0x1C, 0x0, 0x4))

"""
    normalize_hue(h::Real)
    normalize_hue(c::Colorant)

Returns a normalized (wrapped-around) hue angle, or a color with the normalized
hue, in degrees, in [0, 360]. The normalization is essentially equivalent to
`mod(h, 360)`, but is faster at the expense of some accuracy.
"""
@fastmath normalize_hue(h::Real) = max(muladd(floor(h / 360), -360, h), zero(h))
@fastmath normalize_hue(h::Float16) = Float16(normalize_hue(Float32(h)))
normalize_hue(c::C) where {C <: Union{HSV, HSL, HSI}} = C(normalize_hue(c.h), c.s, comp3(c))
normalize_hue(c::C) where {Cb <: Union{HSV, HSL, HSI}, C <: Union{AlphaColor{Cb}, ColorAlpha{Cb}}} =
    C(normalize_hue(c.h), c.s, comp3(c), alpha(c))
normalize_hue(c::C) where C <: Union{LCHab, LCHuv} = C(c.l, c.c, normalize_hue(c.h))
normalize_hue(c::C) where {Cb <: Union{LCHab, LCHuv}, C <: Union{AlphaColor{Cb}, ColorAlpha{Cb}}} =
    C(c.l, c.c, normalize_hue(c.h), c.alpha)

"""
    mean_hue(h1::Real, h2::Real)
    mean_hue(a::C, b::C) where {C <: Colorant}

Compute the mean of two hue angles in degrees.

If the inputs are HSV-like or Lab-like color objects, this will also return a
hue, not a color. If one of the colors is achromatic, i.e. has zero saturation
or chroma, the hue of the other color is returned instead of the mean.

When the hue difference is exactly 180 degrees, which of the two mean hues is
returned depends on the implementation. In other words, it may vary by the color
type and version.
"""
function mean_hue(h1::T, h2::T) where {T <: Real}
    @fastmath hmin, hmax = minmax(h1, h2)
    d = 180 - normalize_hue(hmin - hmax - 180)
    F = typeof(zero(T) / 2)
    mh = muladd(F(0.5), d, hmin)
    return mh < 0 ? mh + 360 : mh
end
@inline function mean_hue(a::C, b::C) where {Cb <: Union{Lab, Luv},
                                     C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    a1, b1, a2, b2 = comp2(a), comp3(a), comp2(b), comp3(b)
    c1, c2 = chroma(a), chroma(b)
    k1 = c1 == zero(c1) ? oneunit(c1) : c1
    k2 = c2 == zero(c2) ? oneunit(c2) : c2
    dp = muladd(a1, a2, b1 * b2) # dot product
    dt = b1 * a2 - a1 * b2
    if dp < zero(dp)
        # 90 deg rotations
        a1, b1 = @fastmath flipsign( b1, dt), flipsign(-a1, dt)
        a2, b2 = @fastmath flipsign(-b2, dt), flipsign( a2, dt)
        if dt == zero(dt) && a1 > zero(a1)
            k1, k2 = -k1, -k2
        end
    end
    ma = muladd(k2, a1, k1 * a2)
    mb = muladd(k2, b1, k1 * b2)
    hue(Cb(zero(ma), ma, mb))
end
function mean_hue(a::C, b::C) where {Cb <: Union{LCHab, LCHuv},
                                     C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    mean_hue(a.c == 0 ? b.h : a.h, b.c == 0 ? a.h : b.h)
end
function mean_hue(a::C, b::C) where {Cb <: Union{HSV, HSL, HSI},
                                     C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    mean_hue(a.s == 0 ? b.h : a.h, b.s == 0 ? a.h : b.h)
end
mean_hue(a, b) = mean_hue(promote(a, b)...)

_delta_h_th(T) = zero(T)
_delta_h_th(::Type{Float32}) = 0.1f0
_delta_h_th(::Type{Float64}) = 6.5e-3

function delta_h(a::C, b::C) where {Cb <: Union{Lab, Luv},
                                    C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    a1, b1, a2, b2 = comp2(a), comp3(a), comp2(b), comp3(b)
    c1, c2 = chroma(a), chroma(b)
    dp = muladd(a1, a2, b1 * b2) # dot product
    dt = b1 * a2 - a1 * b2
    if _delta_h_th(typeof(dp)) * dp <= abs(dt)
        if dp isa Float32 || dp < 0
            dh = @fastmath sqrt(2 * muladd(c1, c2, -dp))
        else
            da, db, dc = a1 - a2, b1 - b2, c1 - c2
            dh = @fastmath sqrt(max(muladd(dc, -dc, muladd(da, da, db^2)), 0))
        end
        return copysign(dh, dt)
    else
        dtf = muladd(b1, a2, -a1 * b2)
        tn = dtf / dp #    tan(Δh) = Δh +  (1/3)*Δh^3 +   (2/15)*Δh^5 + ...
                      # 2sin(Δh/2) = Δh - (1/24)*Δh^3 + (1/1920)*Δh^5 - ...
                      #            ≈ tan(Δh) - (3/8)*tan(Δh)^3 + (31/128)*tan(Δh)^5
        sn = tn * @evalpoly(tn^2, 1, oftype(tn, -3/8), oftype(tn, 31/128)) # 2sin(Δh/2)
        return @fastmath sqrt(c1 * c2) * sn
    end
end
function delta_h(a::C, b::C) where {Cb <: Union{LCHab, LCHuv},
                                    C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    dh0 = hue(a) - hue(b)
    sh = muladd(dh0, oftype(dh0, 1 / 360), oftype(dh0, 0.5))
    d =  @fastmath sh - floor(sh) - oftype(sh, 0.5)
    2 * sqrt(chroma(a) * chroma(b)) * sinpi(d)
end
delta_h(a, b) = delta_h(promote(a, b)...)

@inline function delta_c(a::C, b::C) where {Cb <: Union{Lab{Float32}, Luv{Float32}},
                                    C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    n1, m1 = @fastmath minmax(comp2(a)^2, comp3(a)^2)
    n2, m2 = @fastmath minmax(comp2(b)^2, comp3(b)^2)
    ((m1 - m2) + (n1 - n2)) / @fastmath(max(chroma(a) + chroma(b), floatmin(Float32)))
end
delta_c(a, b) = chroma(a) - chroma(b)

"""
    weighted_color_mean(w1, c1, c2)

Returns the color `w1*c1 + (1-w1)*c2` that is the weighted mean of `c1` and
`c2`, where `c1` has a weight 0 ≤ `w1` ≤ 1.
"""
weighted_color_mean(w1::Real, c1::Colorant, c2::Colorant) = _weighted_color_mean(w1, c1, c2)
function weighted_color_mean(w1::Real, c1::C, c2::C) where {Cb <: Union{HSV, HSL, HSI, LCHab, LCHuv},
                                                            C <: Union{Cb, AlphaColor{Cb}, ColorAlpha{Cb}}}
    normalize_hue(_weighted_color_mean(w1, c1, c2))
end
function _weighted_color_mean(w1::Real, c1::Colorant{T1}, c2::Colorant{T2}) where {T1,T2}
    @fastmath min(w1, oneunit(w1) - w1) >= zero(w1) || throw(DomainError(w1, "`w1` must be in [0, 1]"))
    w2 = oneunit(w1) - w1
    mapc((x, y) -> convert(promote_type(T1, T2), muladd(w1, x, w2 * y)), c1, c2)
end
function _weighted_color_mean(w1::Integer, c1::C, c2::C) where C <: Colorant
    (w1 & 0b1) === w1 || throw(DomainError(w1, "`w1` must be in [0, 1]"))
    w1 == zero(w1) ? c2 : c1
end

"""
    weighted_color_mean(weights, colors)

Returns the weighted mean of the given collection `colors` with `weights`.
This is semantically equivalent to the calculation of `sum(weights .* colors)`.

# Examples
```jldoctest; setup = :(using Colors, FixedPointNumbers)
julia> rgbs = (RGB(1, 0, 0), RGB(0, 1, 0), RGB(0, 0, 1));

julia> weighted_color_mean([0.2, 0.2, 0.6], rgbs)
RGB{N0f8}(0.2,0.2,0.6)

julia> weighted_color_mean(0.5:-0.25:0.0, RGB{Float64}.(rgbs))
RGB{Float64}(0.5,0.25,0.0)
```

!!! compat "Colors v0.13"
    `wighted_color_mean` with collection or iterator inputs requires Colors
    v0.13 or later.

!!! note
    In a cylindrical color space such as HSV, a weighted mean of more than three
    colors is generally not meaningful.
"""
function weighted_color_mean(weights, colors)
    if length(weights) != length(colors)
        throw(DimensionMismatch("`weights` and `colors` must be the same length."))
    end
    w1, c1 = first(weights), first(colors)
    acc = mapc(_ -> zero(promote_type(typeof(w1), Float32)), c1)
    for (w, c) in zip(weights, colors)
        acc = mapc((a, v) -> a + w * v, acc, c)
    end
    convert(typeof(c1), acc)
end

"""
    range(start::T; stop::T, length=100) where T<:Colorant
    range(start::T, stop::T; length=100) where T<:Colorant

Generates N (=`length`) >2 colors in a linearly interpolated ramp from `start` to`stop`,
inclusive, returning an `Array` of colors.

!!! compat "Julia 1.1"
    `stop` as a positional argument requires at least Julia 1.1.
"""
function range(start::T; stop::T, length::Integer=100) where T<:Colorant
    return T[weighted_color_mean(w1, start, stop) for w1 in range(1.0,stop=0.0,length=length)]
end

if VERSION >= v"1.1"
    range(start::T, stop::T; kwargs...) where T<:Colorant = range(start; stop=stop, kwargs...)
end

if VERSION < v"1.0.0-"
import Base: linspace
Base.@deprecate linspace(start::Colorant, stop::Colorant, n::Integer=100) range(start, stop=stop, length=n)
end

# Double quadratic Bezier curve
function bezier(t::T, p0::T, p2::T, q0::T, q1::T, q2::T) where T<:Real
    B(t,a,b,c)=a*(1.0-t)^2 + 2.0b*(1.0-t)*t + c*t^2
    if t <= 0.5
        return B(2.0t, p0, q0, q1)
    else #t > 0.5
        return B(2.0(t-0.5), q1, q2, p2)
    end
end

# Inverse double quadratic Bezier curve
function inv_bezier(t::T, p0::T, p2::T, q0::T, q1::T, q2::T) where T<:Real
    invB(t,a,b,c)=(a-b+sqrt(b^2-a*c+(a-2.0b+c)*t))/(a-2.0b+c)
    if t < q1
        return 0.5*invB(t,p0,q0,q1)
    else #t >= q1
        return 0.5*invB(t,q1,q2,p2)+0.5
    end
end


text:
using Test
using MicroMacro

include("error.jl")

begin

    epsilons = [10.0^(-i) for i in 0:6]

    xmin     = 0
    xmax     = 2π
    T        = 2π
    size_x   = [64]
    size_tau = [64]
    Tfinal   = 0.25

    nb_dt  = 5

    for nx in size_x, ntau in size_tau

        println(rpad("dt", 15),  rpad("epsilon",15), rpad("error",15))

        for (kk, epsilon) in enumerate(epsilons)

            data = CosSin(xmin, xmax, nx, epsilon, T, Tfinal)

            for hh in 1:nb_dt

                dt = 2.0^(-hh) * Tfinal / 16

                micmac = MicMac(data, ntau)

                u, v = solve(micmac, dt)

                err = compute_error(u,v,data)
                println(rpad(dt, 15), rpad(epsilon,15), rpad(err, 15))
                @test isapprox(err, 0.0, atol=1e-5)

            end

        end

    end

end


text:
using Test, PolynomialBases

if !haskey(ENV, "JULIA_PKGEVAL") # sympy is not install on https://github.com/JuliaComputing/NewPkgEval.jl
  import SymPy
  x, α = SymPy.symbols("x, alpha")

  @test 0 == SymPy.simplify( gegenbauer(x, 0, α) - 1 )
  @test 0 == SymPy.simplify( gegenbauer(x, 1, α) - ( 2α*x ) )
  @test 0 == SymPy.simplify( gegenbauer(x, 2, α) - ( -α + 2α*(1+α)*x^2 ) )
  @test 0 == SymPy.simplify( gegenbauer(x, 3, α) - ( -2α*(1+α)*x + 4*α*(1+α)*(2+α)*x^3/3 ) )

  @test_skip gegenbauer(x, 4, α)
end
@inferred gegenbauer(10., 4, 3)


text:
mutable struct Game
    snake::Snake
    food::Point
    direction::Direction
    score::Int

    function Game()
        snake = Snake()
        food = place_food(snake)

        new(snake, food, RIGHT, 0)
    end
end

function Base.show(io::IO, game::Game)
    println("Snake Game")
    println("==========")
    println("Snake -> ", game.snake)
    println("Food -> ", game.food)
    println("Direction -> ", game.direction)
    println("Score -> ", game.score)
end

function place_food(s::Snake)
    x = rand(0:(WIDTH-BLOCK_SIZE)÷BLOCK_SIZE) * BLOCK_SIZE
    y = rand(0:(HEIGHT-BLOCK_SIZE)÷BLOCK_SIZE) * BLOCK_SIZE

    # Create new food
    food = Point(x, y)

    # Look if food is inside the snake
    if food in s.body
        return place_food(s)
    end

    return food
end

function reset!(game::Game)
    snake = Snake()
    food = place_food(snake)

    game.snake = snake
    game.food = food
    game.direction = RIGHT
    game.score = 0
end

function get_state(game::Game)
    head = game.snake.head
    point_l = Point(head.x - BLOCK_SIZE, head.y)
    point_r = Point(head.x + BLOCK_SIZE, head.y)
    point_u = Point(head.x, head.y - BLOCK_SIZE)
    point_d = Point(head.x, head.y + BLOCK_SIZE)

    dir_l = game.direction == LEFT
    dir_r = game.direction == RIGHT
    dir_u = game.direction == UP
    dir_d = game.direction == DOWN

    # Create the state vector
    state = [
        # Danger straight
        (dir_r && is_collision(game.snake, point_r)) ||
        (dir_l && is_collision(game.snake, point_l)) ||
        (dir_u && is_collision(game.snake, point_u)) ||
        (dir_d && is_collision(game.snake, point_d)),

        # Danger right
        (dir_u && is_collision(game.snake, point_r)) ||
        (dir_d && is_collision(game.snake, point_l)) ||
        (dir_l && is_collision(game.snake, point_u)) ||
        (dir_r && is_collision(game.snake, point_d)),

        # Danger left
        (dir_d && is_collision(game.snake, point_r)) ||
        (dir_u && is_collision(game.snake, point_l)) ||
        (dir_r && is_collision(game.snake, point_u)) ||
        (dir_l && is_collision(game.snake, point_d)),

        # Move direction
        dir_l,
        dir_r,
        dir_u,
        dir_d,

        # Food location
        game.food.x < game.snake.head.x,  # food left
        game.food.x > game.snake.head.x,  # food right
        game.food.y < game.snake.head.y,  # food up
        game.food.y > game.snake.head.y   # food down
    ]

    return convert.(Int, state)
end

function send_inputs!(game::Game, move::AbstractArray{<:Integer})
    if move == [1, 0, 0]
        new_dir = game.direction    # No changes in direction
    elseif move == [0, 1, 0]
        idx = Int(game.direction) + 1
        idx > 4 ? idx = 1 : nothing
        new_dir = Direction(idx)    # Turning clockwise
    elseif move == [0, 0, 1]
        idx = Int(game.direction) - 1
        idx < 1 ? idx = 4 : nothing
        new_dir = Direction(idx)    # Turning anticlockwise
    end

    # Perform the move
    game.direction = new_dir
    return
end

function play_step!(g::Game)
    # Move the snake position
    tail = move!(g.snake, g.direction)

    # Check if game is over
    reward = 0
    game_over = false
    # Look for collisions
    if is_collision(g.snake)
        game_over = true
        reward = -10

    elseif g.snake.head == g.food
        # Snake's length grows by one block
        push!(g.snake.body, tail)

        # Adjust the score and place new food
        g.score += 1
        reward = 10
        g.food = place_food(g.snake)
    end

    return reward, game_over, g.score
end


text:
# ---
# title: 754. Reach a Number
# id: problem754
# author: Indigo
# date: 2021-07-11
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/reach-a-number/description/>
# hidden: true
# ---
#
# You are standing at position `0` on an infinite number line. There is a goal
# at position `target`.
#
# On each move, you can either go left or right. During the _n_ -th move
# (starting from 1), you take _n_ steps.
#
# Return the minimum number of steps required to reach the destination.
#
# **Example 1:**
#
#
#
#     Input: target = 3
#     Output: 2
#     Explanation:
#     On the first move we step from 0 to 1.
#     On the second step we step from 1 to 3.
#
#
# **Example 2:**
#
#
#
#     Input: target = 2
#     Output: 3
#     Explanation:
#     On the first move we step from 0 to 1.
#     On the second move we step  from 1 to -1.
#     On the third move we step from -1 to 2.
#
#
# **Note:**
#
# * `target` will be a non-zero integer in the range `[-10^9, 10^9]`.
#
#
## @lc code=start
using LeetCode

function reach_number(target::Integer)
    target = abs(target)
    k = 0
    while target > 0
        k += 1
        target = target - k
    end
    return iseven(target) ? k : k + 1 + k % 2
end
## @lc code=end


text:
include("preprocessing/preprocessing.jl")
include("postprocessing/postprocessing.jl")
include("inprocessing/inprocessing.jl")


text:
function construct_icd10(icd10::AbstractString)
    system = "ICD10"
    value = icd10_add_dots(icd10)
    return DiagnosisClass(; system, value)
end

function construct_ccs10(ccs10::AbstractString)
    system = "CCS10"
    value = ccs10
    return DiagnosisClass(; system, value)
end


text:
using st16SseqJuliaTools
using ArgParse
using CSV
using DataFrames
using DataStructures
using BioSequences
using CodecZlib
using FASTX

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--16sr", "-r"
            help = "16S expected region in a read, <start>:<end> "
            arg_type = String
            required = true
        "--16st", "-t"
            help = "16S expected region in a hmm/alignment, <start>:<end> "
            arg_type = String
            required = true
        "--min", "-m"
            help = "minimum match length "
            arg_type = Int
            required = true
        "--input","-i"
            help = "Input file"
            required = true
        "--log", "-l"
            help = "Log file to write reads on target"
            arg_type = String
            required = true
        "--names", "-n"
            help = "Write out names of reads "
            arg_type = String
            required = true
        "--fastqin", "-q"
            help = "FASTQ.GZi input to filter out reads"
            arg_type = String
            required = true
        "--fastqout", "-o"
            help = "FATQ.GZ output to filter out reads"
            arg_type = String
            required = true
    end
    return parse_args(s)
end

function write_fasta(log_file::String, fastq_file::String, fastq_fileo::String, chosen_reads::Set{String})
    logf = open(log_file, "w")
    ct_t = length(chosen_reads)
    ct_p = 0
    ct_all = 0
    #read in fastq and filter
    reader = FASTQ.Reader(GzipDecompressorStream(open(fastq_file)))
    writer = FASTQ.Writer(GzipCompressorStream(open(fastq_fileo,"w")))
    println("Choosing reads with proper 16S fragments in file: $fastq_file")
    record = FASTQ.Record()

    while !eof(reader)
        read!(reader, record)
        ct_all += 1
        #println(record)
        id = FASTQ.identifier(record)

        if id in chosen_reads
            ct_p += 1
            write(writer, record)
            delete!(chosen_reads, id)
        end
        ## Do something.
    end

    if ct_p != ct_t
        println(stderr, "Read counts are not equal: $ct_t in csv and $ct_p in the fastq filr")
    end
    frac = ct_p/ct_all
    println(stderr, "On 'target' $frac reads")
    println(logf, "$frac")
    close(writer)
    close(logf)
end

if abspath(PROGRAM_FILE) == @__FILE__
    parsed_args = parse_commandline()
    chosen_reads = extract_rRNA_names(parsed_args["names"], parsed_args["16sr"], parsed_args["16st"], parsed_args["min"], parsed_args["input"])
    write_fasta(parsed_args["log"], parsed_args["fastqin"], parsed_args["fastqout"], chosen_reads)
end


text:
using StatsPlots

include("simulation_helpers.jl")
include("../helper_functions.jl")

begin
    begin
        total_simulation_length = 1000
        interval_length = 1000
        job_size = 1
        p1_arrival_rate = 1
        p2_arrival_rate = 1
        private_service_rate = 1 / 1.5
        public_service_rate = 1 / 0.7
        breakin = 3
        pst, qst = equilibrium_probability(
        p1_arrival_rate,
        p2_arrival_rate,
        job_size,
        α -> 1 / (1.5 * (α)),
        α -> 1 / (0.7 * (α)))
    end

event_data, interval_data = simulate_queue(
    total_simulation_length,
    interval_length,
    Exponential(1 / p1_arrival_rate),
    Exponential(1 / p2_arrival_rate),
    Exponential(job_size / public_service_rate),
    Exponential(job_size / private_service_rate),
    pst,
    qst,
)

transform!(event_data, :eventtype => ByRow(x -> split(x, "_")) => [:source, :event])
select!(event_data, Not(:eventtype))

event_data_wide = vcat([match_events(gdf) for gdf in groupby(event_data, :source)]...)
select!(event_data_wide, [:source, :arrival_time, :departure_time])
transform!(event_data_wide,
    [:arrival_time, :departure_time] => ((x, y) -> (y - x)) => :response_time,
    :arrival_time => (x -> cut(x, 0:interval_length:total_simulation_length + interval_length, extend=true)) => :interval)
performance_data = combine(groupby(event_data_wide, [:interval, :source]), :response_time => (x -> mean(skipmissing(x))) => :mean_response_time)

transform!(performance_data, :source => ByRow(x -> split(x, "-")) => [:player, :dummy])
select!(performance_data, Not(:dummy))
groupby(performance_data, [:player, :interval])
interval_data.cut = cut((interval_data.interval .- 1) .* interval_length, 0:interval_length:total_simulation_length + interval_length)
performance_data = leftjoin(performance_data, select(interval_data, Not(:interval)), on = [:interval => :cut])

cost_data = combine(groupby(performance_data, [:player, :interval])) do sdf
    evaluate_cost(sdf)
end
        rename!(cost_data, :x1 => :cost)

begin
    p2 = plot(title = "Cost by Communication Interval", xlabel = "Interval", ylabel = "Cost")
    cost_data = leftjoin(sort!(cost_data, :player), interval_data[:, [:interval, :cut]], on = (:interval => :cut), makeunique=true)
    for g in groupby(cost_data, :player)
        plot!(p2, g.interval_1, g.cost, label = first(g.player), lw = 2)
    end
    p2
end
p2
end

println(median(skipmissing(event_data_wide.response_time)))
println(mean(skipmissing(event_data_wide.response_time)))
histogram(log10.(collect(skipmissing(event_data_wide.response_time))),
    nbins = 50,
    label = false,
    xlabel = "log10(Response Time)",
    ylabel = "Frequency",
    title = "Histogram of Response Times")

ARs = 0.1:0.1:2
total_volume = 2
JSs = 1 / 2 * total_volume ./ ARs
n_reps = 50
dot_every = Int(floor(n_reps / 5))

begin
mean_resval = Float64[]
var_resval = Float64[]

for (AR, JS) in zip(ARs, JSs)
    mean_temp = Float64[]
    var_temp = Float64[]
    private_service_rate = 1 / (1.5 * JS)
    public_service_rate = 1 / (0.7 * JS)

    print("AR = ", AR, " JS = ", JS, " ")
    for i in 1:n_reps
        if mod(i, dot_every) == 0
            print(".")
        end
        event_data, interval_data = simulate_queue(
            total_simulation_length,
            total_simulation_length,
            AR,
            AR,
            1 / 0.7,
            1 / 1.5,
            JS
        )
        transform!(event_data, :eventtype => ByRow(x -> split(x, "_")) => [:source, :event])
        select!(event_data, Not(:eventtype))
        event_data_wide = vcat([match_events(gdf) for gdf in groupby(event_data, :source)]...)
        select!(event_data_wide, [:source, :arrival_time, :departure_time])
            transform!(event_data_wide,
                    [:arrival_time, :departure_time] => ((x, y) -> (y - x)) => :response_time)
        mean_and_var = combine(event_data_wide,
        :response_time => (
            x -> (mean_RT = mean(skipmissing(x)), var_RT = var(skipmissing(x)))
            ) => AsTable
        )

        append!(mean_temp, mean_and_var[1, 1])
        append!(var_temp, mean_and_var[1, 2])
    end
    println()
    append!(mean_resval, mean(mean_temp))
    append!(var_resval, mean(var_temp))
end
end

plot(ARs, JSs, xlabel = "Arrival Rate", ylabel = "Job Size")
mean_plot = plot(JSs, mean_resval, xlabel = "Job Size", yaxis = :log10, ylabel = "Mean Response Time", legend = :none)
sd_plot = plot(JSs, sqrt.(var_resval), xlabel = "Job Size", yaxis = :log10, ylabel = "SD Response Time", legend = :none)

begin
    total_simulation_length = 10000
    interval_length = total_simulation_length / 10
    job_size = 0.5
    p1_arrival_rate = 2
    p2_arrival_rate = 2
    private_service_rate = 1 / (1.5 * job_size)
    public_service_rate = 1 / (0.7 * job_size)
    breakin = 3
    pst, qst = equilibrium_probability(
        p1_arrival_rate,
        p2_arrival_rate,
        job_size,
        α -> 1 / (1.5 * (α)),
        α -> 1 / (0.7 * (α)))
end

# to ensure correctness: calculate_distributions

event_data, interval_data = simulate_queue(
    1000,
    1000,
    1,
    1,
    1 / 0.7,
    1 / 1.5,
    1
)

function evaluate_response_time(event_data)
    transform!(event_data, :eventtype => ByRow(x -> split(x, "_")) => [:source, :event])
    select!(event_data, Not(:eventtype))
    event_data_wide = vcat([match_events(gdf) for gdf in groupby(event_data, :source)]...)
    select!(event_data_wide, [:source, :arrival_time, :departure_time])
    transform!(event_data_wide,
                [:arrival_time, :departure_time] => ((x, y) -> (y - x)) => :response_time)
    df = combine(event_data_wide,
        :response_time => (
            x -> (mean_RT = mean(skipmissing(x)), sd_RT = sqrt(var(skipmissing(x))))
        ) => AsTable
    )
    return df[1, 1], df[1, 2]
end

evaluate_response_time(event_data)

prepared_data = prepare_event_data_for_plot(event_data)
@df prepared_data plot(
    :time,
    :effect_cumsum,
    group = :groupcolour,
)


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export shortest_distance

"""
  `shortest_distance(fst[, s=get_initial(fst); filter=arc->true, reverse=false)`

Computes the shortest distance between the state `s` to every other state.
The shortest distance between `s` and `i` is defined as the sum of the weights of all paths between these states.

Weights are required to be rigth distributive and k-closed.

If `fst` has `N` states a `N`-long vector is returned where the `i`-the element contains the shortest distance between the states `s` and `i`.

If `reversed==true` the shortest distance from the final states is computed.
In this case `s` has no effect since a new initial superstate is added to the reversed WFST.
Here weights are required to be left distributive and k-closed.

See Mohri "Semiring framework and algorithms for shortest-distance problems", Journal of Automata, Languages and Combinatorics 7(3): 321-350, 2002.

"""
function shortest_distance(fst::WFST{W}, s=get_initial(fst); filter=arc->true, reversed=false) where {W}

  if reversed
    if !isright(W)
      throw(ErrorException("Weight $W is not right distributive"))
    end
    if !ispath(W)
      throw(ErrorException("Weight $W must satisfy the \"path\" property"))
    end

    rfst = reverse(fst)
    d = no_check_shortest_distance(rfst,get_initial(rfst);filter=filter)
    return d[2:end]
  else
    if !isleft(W)
      throw(ErrorException("Weight $W is not left distributive"))
    end
    if !ispath(W)
      throw(ErrorException("Weight $W must satisfy the \"path\" property"))
    end
    return no_check_shortest_distance(fst,s;filter=filter)
  end

end

function no_check_shortest_distance(fst::WFST{W,A,D},
                                    s::D=get_initial(fst);
                                    filter=arc->true) where {W,A,D}
  Ns = size(fst,1)
  d = zeros(W,Ns) # d[q] estimate of shortest distance between q and s
  r = zeros(W,Ns) # r[q] total weight added to d[q] since last time q was extracted from S
  d[s], r[s] = one(W), one(W)
  S = Queue{D}()
  enqueue!(S,s)
  while !isempty(S)
    q = dequeue!(S)
    R = r[q]
    r[q] =  zero(W)
    for arc in fst[q]
      if filter(arc)
        n, w = get_nextstate(arc), get_weight(arc)
        dn = d[n]
        if dn != (dn+(R*w))
          d[n] = dn+(R*w)
          r[n] = r[n]+(R*w)
          if !(n in S)
            enqueue!(S,n)
          end
        end
      end
    end
  end
  d[s] = one(W)
  return d
end


text:
module test_05_ArraySignalsWithUnit

using ModiaResult
using Unitful
using OrderedCollections
ModiaResult.@usingModiaPlot

t = range(0.0, stop=1.0, length=100)

result = OrderedDict{String,Any}()

Ibase  = [1.1  1.2  1.3;
          2.1  2.2  2.3;
          3.1  3.2  3.3]u"kg*m^2"

result["time"]     = t*u"s"
result["Inertia"]  = [Ibase*t[i] for i in eachindex(t)]

println("\n... test_05_ArraySignalsWithUnit:")
ModiaResult.printResultInfo(result)

plot(result, ["Inertia[2,2]", "Inertia[2:3,3]"], heading="Array signals")

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
using Test     #src
# # Discrete-time Switched Invariant Set
#
#md # [![Binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/discrete_switched.ipynb)
#md # [![nbviewer](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/discrete_switched.ipynb)
#
# ## Introduction
#
# In this notebook, we compute the maximal (resp. minimal) invariant set contained in the square with vertices $(\pm 1, \pm 1)$ for the system
# ```math
# \begin{aligned}
# x_{k+1} & = -y_k\\
# y_{k+1} & =  x_k.
# \end{aligned}
# ```
#
# The system is $x_{k+1} = Ax_k$ where
# ```math
# A = \begin{bmatrix}
# 0 & -1\\
# 1 & 0
# \end{bmatrix}.
# ```
# A set $S$ is controlled invariant if
# $$AS \subseteq S.$$

# We need to pick an LP and an SDP solver, see [here](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) for a list of available ones. Run one of the following two cells to choose choose the solver.

using SetProg
import Clp
lp_solver = optimizer_with_attributes(Clp.Optimizer, MOI.Silent() => true)
import CSDP
sdp_solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
using Polyhedra
lib = DefaultLibrary{Float64}(lp_solver)
h = HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1)
□ = polyhedron(h, lib)
◇ = polar(□)
A1 = [-1 -1
      -4  0] / 4
A2 = [ 3  3
      -2  1] / 4
A = [0.0 -1.0
     1.0  0.0]

function maximal_invariant(template, heuristic::Function)
    solver = if template isa Polytope
        lp_solver
    else
        sdp_solver
    end
    model = Model(solver)
    @variable(model, S, template)
    @constraint(model, S ⊆ □)
    @constraint(model, A1 * S ⊆ S)
    @constraint(model, A2 * S ⊆ S)
    @objective(model, Max, heuristic(volume(S)))
    optimize!(model)
    @show solve_time(model)
    @show termination_status(model)
    @show objective_value(model)
    return value(S)
end
function minimal_invariant(template, heuristic::Function)
    solver = if template isa Polytope
        lp_solver
    else
        sdp_solver
    end
    model = Model(solver)
    @variable(model, S, template)
    @constraint(model, □ ⊆ S)
    @constraint(model, A1 * S ⊆ S)
    @constraint(model, A2 * S ⊆ S)
    @objective(model, Min, heuristic(volume(S)))
    optimize!(model)
    @show solve_time(model)
    @show termination_status(model)
    @show objective_value(model)
    @show JuMP.objective_sense(model)
    return value(S)
end

using Plots
function hexcolor(rgb::UInt32)
    r = ((0xff0000 & rgb) >> 16) / 255
    g = ((0x00ff00 & rgb) >>  8) / 255
    b = ((0x0000ff & rgb)      ) / 255
    Plots.RGBA(r, g, b)
end # Values taken from http://www.toutes-les-couleurs.com/code-couleur-rvb.php
lichen = hexcolor(0x85c17e)
canard = hexcolor(0x048b9a)
aurore = hexcolor(0xffcb60)
frambo = hexcolor(0xc72c48)

function primal_plot(min_sol, max_sol; npoints=256, args...)
    plot(ratio=:equal, tickfont=Plots.font(12); args...)
    plot!(min_sol, color=canard, npoints=npoints)
    plot!(□, color=lichen)
    plot!(max_sol, color=aurore, npoints=npoints)
end
function polar_plot(min_sol, max_sol; npoints=256, args...)
    plot(ratio=:equal, tickfont=Plots.font(12); args...)
    plot!(polar(max_sol), color=aurore, npoints=npoints)
    plot!(◇, color=lichen)
    plot!(polar(min_sol), color=canard, npoints=npoints)
end

# ## Polyhedral template
#
# ### Fixed point approach
#
# This section shows fixed point approach for computing a polyhedral maximal invariant set.
# We implement the fixed point iteration of the standard viability kernel algorithm with the following function.
# A more sophisticated approach specialized for discrete-time linear switched systems
# is available in [SwitchOnSafety](https://github.com/blegat/SwitchOnSafety.jl).

function backward_fixed_point_iteration(set::Polyhedron)
    new_set = set ∩ (A1 \ set) ∩ (A2 \ set)
    removehredundancy!(new_set)
    return new_set
end

# We start with $[-1, 1]^2$ and obtain the maximal invariant set after one iteration.

max_polytope_1 = backward_fixed_point_iteration(□)

# We can see that it has already converged.

backward_fixed_point_iteration(max_polytope_1)

# We now turn to the minimal invariant set with the following iteration:

function forward_fixed_point_iteration(set::Polyhedron)
    new_set = convexhull(set, A1 * set, A2 * set)
    removevredundancy!(new_set)
    return new_set
end

# We start with $[-1, 1]^2$ and obtain the minimal invariant set after three iteration.

min_polytope_1 = forward_fixed_point_iteration(□)

# We see below that the second iteration was indeed needed.

min_polytope_2 = forward_fixed_point_iteration(min_polytope_1)

# And the third iteration as well.

min_polytope_3 = forward_fixed_point_iteration(min_polytope_2)

# We see below that it has converged

forward_fixed_point_iteration(min_polytope_3)

# Let's see this visually

plot(ratio=:equal, tickfont=Plots.font(12))
plot!(min_polytope_3, color=frambo)
plot!(min_polytope_2, color=canard)
plot!(min_polytope_1, color=aurore)
plot!(□, color=lichen)
plot!(max_polytope_1, color=frambo)

# ### Linear programming approach
#
# As introduced in [R21], we can search over all polyhedra with a given face fan.
# The partition can be defined by the face fan of a given polytope as described in [Eq. (4.6), R21].
# For `maximal_invariant`, because the sets are represented with the support function,
# the fan will be used for the face fan of the polar set
# the face fan of the primal set will depend on the value of
# the support function on each piece.
# For `minimal_invariant`, the sets are represented with the gauge function,
# the fac will therefore be used for the face fan of the primal set.
#
# [R21] Raković, S. V.
# *Control Minkowski–Lyapunov functions*
# Automatica, Elsevier BV, 2021, 128, 109598

# We use the face fan of the polytope `pieces8` (which has the same face fan as the octogon `convexhull(□, √2 * ◇)`).

pieces8 = convexhull(□, 1.25 * ◇)
max_polytope = maximal_invariant(Polytope(symmetric=true, piecewise=pieces8), L1_heuristic)
min_polytope = minimal_invariant(Polytope(symmetric=true, piecewise=polar(pieces8)), L1_heuristic)
primal_plot(min_polytope, max_polytope)

# The polar plot is as follows:

polar_plot(min_polytope, max_polytope)

# ## Ellipsoidal template

# We now consider the ellipsoidal template.
# The ellipsoids of maximal volume are given as follows:

max_ell_vol = maximal_invariant(Ellipsoid(symmetric=true), nth_root)
min_ell_vol = minimal_invariant(Ellipsoid(symmetric=true), nth_root)
primal_plot(min_ell_vol, max_ell_vol)

# All ellipsoids are convex so we can also obtain the polar plot:

polar_plot(min_ell_vol, max_ell_vol)

# The ellipsoids of maximal sum of squared length of its axis.

max_ell_L1 = maximal_invariant(Ellipsoid(symmetric=true), vol -> L1_heuristic(vol, ones(2)))
min_ell_L1 = minimal_invariant(Ellipsoid(symmetric=true), vol -> L1_heuristic(vol, ones(2)))
primal_plot(min_ell_L1, max_ell_L1)

# And the polar plot is:

polar_plot(min_ell_vol, max_ell_vol)

# ## Polyset template

# We generalize the previous template to homogeneous polynomials of degree $2d$.
# As objective, we use the integral of the polynomial used to represent
# the set (either in the support function for `maximal_invariant`
# or in the gauge function for `minimal_invariant`) over the square `□`.
# We start with quartic polynomials.

max_4 = maximal_invariant(PolySet(symmetric=true, convex=true, degree=4), vol -> L1_heuristic(vol, ones(2)))
min_4 = minimal_invariant(PolySet(symmetric=true, degree=4), vol -> L1_heuristic(vol, ones(2)))
primal_plot(min_4, max_4)

# We now obtain the results for sextic polynomials, note that we do not
# impose the sets to be convex for `minimal_invariant`. This is why we don't
# show the polar plot for polysets.

max_6 = maximal_invariant(PolySet(symmetric=true, convex=true, degree=6), vol -> L1_heuristic(vol, ones(2)))
min_6 = minimal_invariant(PolySet(symmetric=true, degree=6), vol -> L1_heuristic(vol, ones(2)))
primal_plot(min_6, max_6)

# The nonconvexity is even more apparent for octic polynomials.

max_8 = maximal_invariant(PolySet(symmetric=true, convex=true, degree=8), vol -> L1_heuristic(vol, □))
min_8 = minimal_invariant(PolySet(symmetric=true, degree=8), vol -> L1_heuristic(vol, □))
primal_plot(min_8, max_8, npoints=1024)


text:
using PhaseRetrieval
using SampledDomains
using Test

@testset "utils.jl" begin
    # Write your tests here.
    q1d = PhaseRetrieval.quadratic(5)
    q2d = PhaseRetrieval.quadratic((5, 5))
    @test q1d == q2d ==  [
    8.0  5.0  4.0  5.0  8.0
    5.0  2.0  1.0  2.0  5.0
    4.0  1.0  0.0  1.0  4.0
    5.0  2.0  1.0  2.0  5.0
    8.0  5.0  4.0  5.0  8.0]

    dm = PhaseRetrieval.diskmatrix(10, .8)
    @test dm[1] == [
        0  0  0  0  0  0  0  0  0  0
        0  0  0  0  1  1  0  0  0  0
        0  0  1  1  1  1  1  1  0  0
        0  0  1  1  1  1  1  1  0  0
        0  1  1  1  1  1  1  1  1  0
        0  1  1  1  1  1  1  1  1  0
        0  0  1  1  1  1  1  1  0  0
        0  0  1  1  1  1  1  1  0  0
        0  0  0  0  1  1  0  0  0  0
        0  0  0  0  0  0  0  0  0  0
            ]
    # @test dm[2] == [
    #     NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
    #     NaN  NaN  NaN  NaN    1    1  NaN  NaN  NaN  NaN
    #     NaN  NaN    1    1    1    1    1    1  NaN  NaN
    #     NaN  NaN    1    1    1    1    1    1  NaN  NaN
    #     NaN    1    1    1    1    1    1    1    1  NaN
    #     NaN    1    1    1    1    1    1    1    1  NaN
    #     NaN  NaN    1    1    1    1    1    1  NaN  NaN
    #     NaN  NaN    1    1    1    1    1    1  NaN  NaN
    #     NaN  NaN  NaN  NaN    1    1  NaN  NaN  NaN  NaN
    #     NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN  NaN
    # ]

    @test PhaseRetrieval.lineararray(-1:.5:1, -1:.25:1, 3, 2, 1) == [
        -4.0  -2.5  -1.0  0.5  2.0
        -3.5  -2.0  -0.5  1.0  2.5
        -3.0  -1.5   0.0  1.5  3.0
        -2.5  -1.0   0.5  2.0  3.5
        -2.0  -0.5   1.0  2.5  4.0
        -1.5   0.0   1.5  3.0  4.5
        -1.0   0.5   2.0  3.5  5.0
        -0.5   1.0   2.5  4.0  5.5
         0.0   1.5   3.0  4.5  6.0
    ]

    @test PhaseRetrieval.lineararray(5, .5, 1, -1.5) == [
        0.0  0.5  1.0  1.5  2.0
        1.0  1.5  2.0  2.5  3.0
        2.0  2.5  3.0  3.5  4.0
        3.0  3.5  4.0  4.5  5.0
        4.0  4.5  5.0  5.5  6.0
    ]

    @test PhaseRetrieval.subdivide_sum(ones(10, 10), 5) == [
        25.0  25.0
        25.0  25.0
    ]

    a = SampledDomains.CartesianDomain2D(1:5, -2.5:.5:2.5)
    ap = PhaseRetrieval.aperture(a, 2.5, (3, -0.1))
    b = phwrap.(ap[2])
    # b  == [
    #     NaN  NaN    NaN    NaN    NaN
    #     NaN  NaN    NaN    NaN    NaN
    #     NaN  NaN    NaN    NaN    NaN
    #     NaN  NaN      1.0  NaN    NaN
    #     NaN    1.0    1.0    1.0  NaN
    #     NaN    1.0    1.0    1.0  NaN
    #     NaN    1.0    1.0    1.0  NaN
    #     NaN  NaN      1.0  NaN    NaN
    #     NaN  NaN    NaN    NaN    NaN
    #     NaN  NaN    NaN    NaN    NaN
    #     NaN  NaN    NaN    NaN    NaN
    #  ]

    @test isnan.(b) == [
        1  1  1  1  1
        1  1  1  1  1
        1  1  1  1  1
        1  1  0  1  1
        1  0  0  0  1
        1  0  0  0  1
        1  0  0  0  1
        1  1  0  1  1
        1  1  1  1  1
        1  1  1  1  1
        1  1  1  1  1
    ]

   a =  PhaseRetrieval.tile(reshape(1:60, (6, 10)), 2)
   @test a[:,:,1] == [ 1 7;2 8] && a[:,:,15] == [53 59; 54 60]

   a = 1:(4 * 3 * 3 * 3);
   b = reshape(a, 4 * 3, 3 * 3);
   c = PhaseRetrieval.subdivide(b, 3);
   d = PhaseRetrieval.subdivide_sum(b,3);

   @test  c[:,:,4] == [
    37  49  61
    38  50  62
    39  51  63
    40  52  64
   ]

   @test d == [
    369  477  585
    378  486  594
    387  495  603
    396  504  612
   ]

end


text:
module GatherShot

include("patches.jl")
include("observe.jl")
include("coverage_matrix.jl")

export generate_reports
export read_reports
export select_tests
export empty_return
export flip_return

end