text:
export deepmap_parser

"""
    deepmap_parser(f::Function,x,a...;kw...)

Perform a deep transformation of a `x`.
Used for [`log_names`](@ref).

Calls `deepmap_parser(f,IdDict(),x,a...)`.
"""
deepmap_parser(f::Function,x::CombinedParser, a...;kw...) =
    deepmap_parser(f,IdDict(),x,a...;kw...)

"""
    deepmap_parser(f::Function,mem::AbstractDict,x::LeafParser,a...;kw...)

return
```julia
    get!(mem,x) do
        f(x,a...;kw...)
    end
```
"""
deepmap_parser(f::Function,mem::AbstractDict,x::LeafParser,a...;kw...) =
    get!(mem,x) do
        f(x,a...;kw...)
    end

deepmap_parser(f::Function,mem::AbstractDict,x::SideeffectParser,a...;kw...) =
    get!(mem,x) do
        SideeffectParser(x.effect,
                         deepmap_parser(f,mem,x.parser,a...;kw...),
                         x.args...)
    end

deepmap_parser(f::Function,mem::AbstractDict,x::FlatMap,a...;kw...) =
    get!(mem,x) do
        FlatMap{result_type(x)}(x.right,# v -> deepmap_parser(f,mem,x.right(v),a...;kw...),
                                deepmap_parser(f,mem,x.left,a...;kw...))
    end

deepmap_parser(f::Function,mem::AbstractDict,x::Sequence,a...;kw...) =
    get!(mem,x) do
        Sequence( ( deepmap_parser(f,mem,p,a...;kw...)
                    for p in x.parts)... )
    end

deepmap_parser(f::Function,mem::AbstractDict,x::Repeat,a...;kw...) =
    get!(mem,x) do
        f(Repeat(x.range,
                 deepmap_parser(f,mem,x.parser,a...)),a...;kw...)
    end

deepmap_parser(f::Function,mem::AbstractDict,x::Optional,a...;kw...) =
    get!(mem,x) do
        Optional(deepmap_parser(f,mem,x.parser,a...;kw...);
                 default=x.default)
    end

function deepmap_parser(f::Function,mem::AbstractDict,x::Either{<:Vector},a...;kw...)
    if haskey(mem,x)
        mem[x]
    else
        mem[x] = r = Either{result_type(x)}(Any[])
        ## f(x,a...)
        for p in x.options
            push!(r,deepmap_parser(f,mem,p,a...;kw...))
        end
        r
    end
end

function deepmap_parser(f::Function,mem::AbstractDict,x::Either{<:Tuple},a...;kw...)
    get!(mem,x) do
        Either((deepmap_parser(f,mem,p,a...;kw...) for p in x.options)... )
    end
end

deepmap_parser(f::Function,mem::AbstractDict,x::Atomic,a...;kw...) =
    get!(mem,x) do
        Atomic(
            deepmap_parser(f,mem,x.parser,a...;kw...))
    end

deepmap_parser(f::Function,mem::AbstractDict,x::Lazy,a...;kw...) =
    get!(mem,x) do
        Lazy(deepmap_parser(f,mem,x.parser,a...;kw...))
    end

"""
    deepmap_parser(f::Function,mem::AbstractDict,x,a...;kw...)

Perform a deep transformation of a CombinedParser.

!!! note
    For a custom parser `P<:CombinedParser` with sub-parsers, provide a method
    ```julia
    deepmap_parser(f::Function,mem::AbstractDict,x::P,a...;kw...) =
        get!(mem,x) do
            ## construct replacement, e.g. if P <: WrappedParser
            P(deepmap_parser(f,mem,x.parser,a...;kw...))
        end
    ```
"""
deepmap_parser(f::Function,mem::AbstractDict,x::CombinedParser,a...;kw...) =
    error("""
    For a custom parser `$(typeof(x))` with sub-parsers, provide a method
    ```julia
    deepmap_parser(f::$(typeof(f)),mem::AbstractDict,x::$(typeof(x)),a...;kw...) =
        get!(mem,x) do
            ## construct replacement, e.g. if P <: WrappedParser
            $(typeof(x))(deepmap_parser(f,mem,x.parser,a...;kw...))
        end
    ```
""")

deepmap_parser(f::Function,mem::AbstractDict,x::NamedParser,a...;kw...) =
    get!(mem,x) do
        NamedParser(x.name,deepmap_parser(f,mem,x.parser,a...;kw...))
    end


text:
# ============================================================================ #
mutable struct Board
    hf::PyPlot.Figure
    ax::PyCall.PyObject
    squares::PatchDict
    coins::PatchDict
    player::Player
    counter::Int
end

function Board(hf, ax)
    if !isassigned(BOARD, 1)
        p = Player(1, 1, East)
        b = Board(hf, ax, PatchDict(), PatchDict(), p, 0)
        BOARD[1] = b
    else
        b = BOARD[1]
        b.player = Player(1, 1, East)
        b.hf = hf
        b.ax = ax
    end
    return b
end

increment!(b::Board) = begin b.counter += 1 end

function set_color(b::Board, x::Integer, y::Integer, col::AbstractString)
    if haskey(b.squares, (x, y))
        b.squares[(x,y)].set_facecolor(col)
    end
    return nothing
end

function get_color(b::Board, x::Integer, y::Integer)
    !haskey(b.squares, (x, y)) && error("Invalid patch corrdinates!")
    return colorname(b.squares[(x,y)].get_facecolor())
end

get_color(b::Board) = get_color(b, location(b.player)...)

function set_player(b::Board, x::Integer, y::Integer, facing::Bearing)
    set_to_face!(b.player, facing)
    moveto!(b.player, x, y)
    b.ax.add_patch(b.player.h)
    return nothing
end

function set_square(b::Board, x::Integer, y::Integer, color::AbstractString="gray")
    p = matplotlib.patches.Rectangle((x,y), 1, 1, color=color)
    b.squares[(x,y)] = p
    b.ax.add_patch(p)
    return nothing
end

function set_coin(b::Board, x::Integer, y::Integer)
    p = matplotlib.patches.CirclePolygon((x+.5,y+.5), .2, color="gold")
    b.coins[(x,y)] = p
    b.ax.add_patch(p)
    return nothing
end

function turn(b::Board, d::Direction)
    game_over() && return nothing
    increment!(b)
    turn!(b.player, d)
    sleep(1/SPEED)
    return nothing
end

@inline left(b::Board) = turn(b, Left)
@inline right(b::Board) = turn(b, Right)

function forward(b::Board)
    game_over() && return nothing
    increment!(b)
    dst = Tuple{Int,Int}(destination(b.player))
    if haskey(b.squares, dst)
        moveto!(b.player, dst...)
        b.hf.canvas.draw()

        if haskey(b.coins, dst)
            sleep(1/(SPEED*2))
            b.coins[dst].remove()
            delete!(b.coins, dst)
            if isempty(b.coins)
                message_box(b, "Final score: $(b.counter)!")
            end
        end
    end
    sleep(1/SPEED)
    return nothing
end

coins_remaining(b::Board) = length(b.coins)

@inline left() = isassigned(BOARD, 1) && left(BOARD[1])
@inline right() = isassigned(BOARD, 1) && right(BOARD[1])
@inline forward() = isassigned(BOARD, 1) && forward(BOARD[1])
@inline get_color() = isassigned(BOARD, 1) && get_color(BOARD[1])
@inline coins_remaining() = isassigned(BOARD, 1) && coins_remaining(BOARD[1])

@inline red() = get_color() == "red"
@inline green() = get_color() == "green"
@inline blue() = get_color() == "blue"
@inline gray() = get_color() == "gray"

function game_over()
    if isassigned(BOARD, 1)
        return isempty(BOARD[1].coins) || BOARD[1].counter > ceil(Int, length(BOARD[1].squares)* 1.5)
    else
        return true
    end
end

function refresh()
    if isassigned(RAW_STR, 1)
        make_board(RAW_STR[1])
    end
    return nothing
end

function clear!(b::Board)
    if plt.fignum_exists(b.hf.number)
        clear!(b.squares)
        clear!(b.coins)
        clear!(b.player)
    end
    b.counter = 0
    return b
end

# ============================================================================ #

ext: jl


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

name = "Zellij"
version = v"0.24.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/zellij-org/zellij/archive/refs/tags/v$(version).tar.gz",
                  "a7f2d1fa1dd9c55d37d1daebdf6af3c6666d144ee1e85ac7f805544ae03e3b1e"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/zellij*/
cargo build --release
install -Dvm 755 "target/${rust_target}/release/zellij${exeext}" "${bindir}/zellij${exeext}"
"""

platforms = [
    Platform("i686", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("aarch64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("aarch64", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("aarch64", "macos"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("zellij", :zellij),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", compilers=[:c, :rust],
               # Wasmer wants to use `-march`, sigh.
               lock_microarchitecture=false,
               )

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
using LinearAlgebra, Statistics
using ThreadTools

sample_eachcol_type = SubArray{Float64, 1, Matrix{Float64}, Tuple{Base.Slice{Base.OneTo{Int64}}, Int64}, true}

# GMM
@doc """
Function to initialize Gaussian mixture model with 'K' components and input vectors of 'dim' dimension.

# Examples
```jldoctest
Θ, gmm, gm = create_gmm(K, dim)
```
where the 'Θ' is dictionary of params, 'gmm' is Gaussian mixture model and 'gm' is a single Gaussian model.
""" ->
function create_gmm(K::Int64, dim::Int64)
    μ_coefs::Vector{Vector{Float64}} = [c[:] for c in eachcol(randn(dim, K))]
    Σ_coefs::Vector{Matrix{Float64}} = [randn(dim, dim) for _ in 1:K]
    α_coefs::Vector{Float64} = randn(K)
    coefs::Dict{Symbol, Vector} = Dict(:μ=>μ_coefs, :Σ=>Σ_coefs, :α=>α_coefs)

    function gmm(Θ::Dict{Symbol, Vector}, x::Vector{Float64})
        n = size(x)[1]
        return mapreduce(usa->usa[3]*(1/((2*π)^(n/2)*det(usa[2])^(1/2)))*
                            (ℯ^(-((x-usa[1])'*inv(usa[2])*(x-usa[1])/2))), +,
                                zip(Θ[:μ], Θ[:Σ], Θ[:α]))
    end
    gm(μ::Vector{Float64},
        Σ::Matrix{Float64},
        x::Vector{Float64}) = (1/((2*π)^(size(x)[1]/2)*det(Σ)^(1/2)))*(ℯ^(-((x-μ)'*inv(Σ)*(x-μ)/2)))
    return coefs, gmm, gm
end

# EM algorithm
@doc """
Function to learn params of Gaussian mixture model with the usage of Expectation Maximization Algorithm.

# Examples
```jldoctest
EM!(Θ, X, K, gmm, gm, steps)
```
where the 'Θ' is dictionary of params to be learned,
'X' is set of train data (size(dim, N) ~ (dimension of data, data count),
'gmm' is Gaussian mixture model, 'gm' is a single Gaussian model and 'steps' is count of algorithm steps.
""" ->
function EM!(Θ::Dict{Symbol, Vector}, X::Matrix{Float64}, K::Int64,
             gmm, gm, steps::Int64)
    # Initialize
    N::Int64 = size(X)[2]
    Θ[:Σ] = fill(cov(X, dims=2), K)
    Θ[:α] = fill(1/K, K)
    Θ[:μ] = [X[:, k] for k::Int64 in 1:K]

    saturate_cov(Σ) = det(Σ) < 10e-10 ? Σ+I*10e-5 : Σ
    Θ[:Σ] = tmap(x->saturate_cov(x), Θ[:Σ])

    E_step(k::Int64, x::Vector{Float64}) = Θ[:α][k]*gm(Θ[:μ][k], Θ[:Σ][k], x)/gmm(Θ, x)

    for _=1:steps
        # E-step
        γ = mapreduce(permutedims, vcat, [tmap(xᵢ::sample_eachcol_type->
                E_step(k, Vector{Float64}(xᵢ)), eachcol(X)) for k::Int64 in 1:K])'
        # M-step
        Θ[:α] = vec(mean(γ, dims=2))
        Γ = sum(γ, dims=1)[1, :]
        for k::Int64=1:K
            μₖ::Vector{Float64} = mapreduce(ix->γ[ix[1], k]*
                    Vector{Float64}(ix[2]), +, zip(1:N, eachcol(X)))/Γ[k]
            Σₖ::Matrix{Float64} = mapreduce(ix->γ[ix[1], k]*
                    (Vector{Float64}(ix[2])-μₖ)*(Vector{Float64}(ix[2])-μₖ)', +, zip(1:N, eachcol(X)))/Γ[k]
            Θ[:Σ][k] = saturate_cov(Σₖ)
            Θ[:μ][k] = μₖ
        end
    end
    Θ
end

max_stars_repo_name: kozvojtex/AnomalyDetection

max_issues_repo_name: kozvojtex/AnomalyDetection

max_forks_repo_name: kozvojtex/AnomalyDetection


text:
using DataFrames
using CSV

# Appending classification to the CSV data then merging into one file
# -----------------------------------------------------------------------------

# 1 (Sleep stage = W = wake)
dfW = CSV.File("wav/W/audio_recoded.csv") |> DataFrame
dfW[!,"Class"]  = vec(fill(1,(1,size(dfW)[1])))
#show(dfW)

# 2 (Sleep stage = 1)
df1 = CSV.File("wav/1/audio_recoded.csv") |> DataFrame
df1[!,"Class"]  = vec(fill(2,(1,size(df1)[1])))
#show(df1)

# 3 (Sleep stage = 2)
df2 = CSV.File("wav/2/audio_recoded.csv") |> DataFrame
df2[!,"Class"]  = vec(fill(3,(1,size(df2)[1])))
#show(df2)

# 4 (Sleep stage = 3 and 4)
df34 = CSV.File("wav/34/audio_recoded.csv") |> DataFrame
df34[!,"Class"] = vec(fill(4,(1,size(df34)[1])))
#show(df34)

# 5 (Sleep stage = R = REM)
dfR = CSV.File("wav/R/audio_recoded.csv") |> DataFrame
dfR[!,"Class"]  = vec(fill(5,(1,size(dfR)[1])))
#show(dfR)

dfCombined = vcat(dfW, df1, df2, df34, dfR)
#show(dfCombined)
CSV.write("audio_recoded_combined_classed.csv", dfCombined)

# Alternative with binary classification (0 = not REM, 1 = REM)
# -----------------------------------------------------------------------------
# 1 (Sleep stage = W = wake)
dfW = CSV.File("wav/W/audio_recoded.csv") |> DataFrame
dfW[!,"Class"]  = vec(fill(0,(1,size(dfW)[1])))
# 2 (Sleep stage = 1)
df1 = CSV.File("wav/1/audio_recoded.csv") |> DataFrame
df1[!,"Class"]  = vec(fill(0,(1,size(df1)[1])))
# 3 (Sleep stage = 2)
df2 = CSV.File("wav/2/audio_recoded.csv") |> DataFrame
df2[!,"Class"]  = vec(fill(0,(1,size(df2)[1])))
# 4 (Sleep stage = 3 and 4)
df34 = CSV.File("wav/34/audio_recoded.csv") |> DataFrame
df34[!,"Class"] = vec(fill(0,(1,size(df34)[1])))
# 5 (Sleep stage = R = REM)
dfR = CSV.File("wav/R/audio_recoded.csv") |> DataFrame
dfR[!,"Class"]  = vec(fill(1,(1,size(dfR)[1])))
# Combining
dfCombinedBinary = vcat(dfW, df1, df2, df34, dfR)
#show(dfCombined)
CSV.write("audio_recoded_combined_binary.csv", dfCombinedBinary)

# Alternative as 2 files (positives and negatives)
# -----------------------------------------------------------------------------
# 1 (Sleep stage = W = wake)
dfW = CSV.File("wav/W/audio_recoded.csv") |> DataFrame
# 2 (Sleep stage = 1)
df1 = CSV.File("wav/1/audio_recoded.csv") |> DataFrame
# 3 (Sleep stage = 2)
df2 = CSV.File("wav/2/audio_recoded.csv") |> DataFrame
# 4 (Sleep stage = 3 and 4)
df34 = CSV.File("wav/34/audio_recoded.csv") |> DataFrame
# 5 (Sleep stage = R = REM)
dfR = CSV.File("wav/R/audio_recoded.csv") |> DataFrame
# Combining
dfNegative = vcat(dfW, df1, df2, df34)
dfPositive = dfR
dfBareAll  = vcat(dfW, df1, df2, df34, dfR)
CSV.write("audio_recoded_combined_NonREM.csv", dfNegative)
CSV.write("audio_recoded_combined_REM.csv",    dfPositive)
CSV.write("audio_recoded_combined_bare.csv",   dfBareAll)


text:
using LazySets, MathematicalPredicates, SparseArrays

function beam_specification()
    # initial set: xᵢ = 0 if i ≤ 300 and xᵢ ∈ [0.0015, 0.002] otherwise
    X0 = Hyperrectangle(low=[zeros(300); fill(0.0015, 48)],
                        high=[zeros(300); fill(0.002, 48)])

    # safety property: x89 ≤ 2100
    property = is_contained_in(HalfSpace(sparsevec([89], [1.0], 348), 2100.0))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
"""
    integrate(f,x,w,j)
"""
function integrate(f,x,w,j)
  fx = evaluate_field(f,x)
  jx = evaluate_field(j,x)
  k = IntKernel()
  apply_kernel(k,fx,w,jx)
end

"""
    integrate(f::AbstractArray,x,w,j)
"""
function integrate(
  f::AbstractArray,x,w,j)
  fx = evaluate_field_array(f,x)
  jx = evaluate_field_array(j,x)
  k = IntKernel()
  apply(k,fx,w,jx)
end

struct IntKernel <: Kernel end

function kernel_cache(k::IntKernel,f::AbstractVector,w,j)
  T = _integrate_rt(f,w,j)
  zero(T)
end

function kernel_testitem!(c,k::IntKernel,f::AbstractVector,w,j)
  if _integrate_valid_sizes(f,w,j)
    apply_kernel!(c,k,f,w,j)
  else
    c
  end
end

@noinline function apply_kernel!(z,k::IntKernel,f::AbstractVector,w,j)
  _integrate_checks(f,w,j)
  r = z
  for p in eachindex(f)
    @inbounds r = r + f[p]*w[p]*meas(j[p])
  end
  r
end

function kernel_cache(k::IntKernel,f::AbstractArray,w,j)
  T = _integrate_rt(f,w,j)
  _, s = _split(size(f)...)
  r = zeros(T,s)
  c = CachedArray(r)
end

function kernel_testitem!(c,k::IntKernel,f::AbstractArray,w,j)
  if _integrate_valid_sizes(f,w,j)
    apply_kernel!(c,k,f,w,j)
  else
    c.array
  end
end

@inline function apply_kernel!(c,k::IntKernel,f::AbstractArray,w,j)
  _integrate_checks(f,w,j)
  np, s = _split(size(f)...)
  cis = CartesianIndices(s)
  setsize!(c,s)
  z = zero(eltype(c))
  r = c.array
  for i in cis
    @inbounds r[i] = z
  end
  for p in 1:np
    @inbounds dV = w[p]*meas(j[p])
    for i in cis
      @inbounds r[i] += f[p,i]*dV
    end
  end
  r
end

function _integrate_rt(f,w,j)
  Tf = eltype(f)
  Tw = eltype(w)
  Tj = eltype(j)
  return_type(*,Tf,Tw,return_type(meas,Tj))
end

function _integrate_checks(f,w,j)
  @assert _integrate_valid_sizes(f,w,j) "integrate: sizes  mismatch."
end

function _integrate_valid_sizes(f,w,j)
  nf, = size(f)
  nw = length(w)
  nj = length(j)
  (nf == nw) && (nw == nj)
end


text:
# Taken from https://github.com/mwarusz/Atum.jl
# Credit to mwarusz and lcw

function solve!(q, timeend, solver;
  after_step::Function=(x...) -> nothing,
  after_stage::Function=(x...) -> nothing,
  adjust_final=true)
  finalstep = false
  step = 0
  while true
    step += 1
    time = solver.time
    if time + solver.dt >= timeend
      adjust_final && (solver.dt = timeend - time)
      finalstep = true
    end
    dostep!(q, solver)
    after_step(step, time, q)
    finalstep && break
  end
end

mutable struct LSRK{FT,AT,NS,RHS}
  dt::FT
  time::FT
  rhs!::RHS
  dq::AT
  rka::NTuple{NS,FT}
  rkb::NTuple{NS,FT}
  rkc::NTuple{NS,FT}

  function LSRK(rhs!, rka, rkb, rkc, q, dt, t0)
    FT = eltype(eltype(q))
    dq = copy(q)
    fill!(dq, zero(eltype(q)))
    AT = typeof(dq)
    RHS = typeof(rhs!)
    new{FT,AT,length(rka),RHS}(FT(dt), FT(t0), rhs!, dq, rka, rkb, rkc)
  end
end

function dostep!(q, lsrk::LSRK)
  (; rhs!, dq, rka, rkb, rkc, dt, time) = lsrk
  for stage = 1:length(rka)
    stagetime = time + rkc[stage] * dt
    dq .*= rka[stage]
    rhs!(dq, q, stagetime)
    @. q += rkb[stage] * dt * dq
  end
  lsrk.time += dt
end

function LSRK144(rhs!, q, dt; t0=0)
  rka, rkb, rkc = coefficients_lsrk144()
  LSRK(rhs!, rka, rkb, rkc, q, dt, t0)
end

function coefficients_lsrk144()
  rka = (
    0,
    -0.7188012108672410,
    -0.7785331173421570,
    -0.0053282796654044,
    -0.8552979934029281,
    -3.9564138245774565,
    -1.5780575380587385,
    -2.0837094552574054,
    -0.7483334182761610,
    -0.7032861106563359,
    0.0013917096117681,
    -0.0932075369637460,
    -0.9514200470875948,
    -7.1151571693922548,
  )

  rkb = (
    0.0367762454319673,
    0.3136296607553959,
    0.1531848691869027,
    0.0030097086818182,
    0.3326293790646110,
    0.2440251405350864,
    0.3718879239592277,
    0.6204126221582444,
    0.1524043173028741,
    0.0760894927419266,
    0.0077604214040978,
    0.0024647284755382,
    0.0780348340049386,
    5.5059777270269628,
  )

  rkc = (
    0,
    0.0367762454319673,
    0.1249685262725025,
    0.2446177702277698,
    0.2476149531070420,
    0.2969311120382472,
    0.3978149645802642,
    0.5270854589440328,
    0.6981269994175695,
    0.8190890835352128,
    0.8527059887098624,
    0.8604711817462826,
    0.8627060376969976,
    0.8734213127600976,
  )

  rka, rkb, rkc
end

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
module EdgeTest

using Test
using ForneyLab: Interface, Edge, Variable, Interface, FactorNode, FactorGraph, currentGraph, addNode!, disconnect!, generateId

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    function MockNode(; id=generateId(MockNode))
        self = new(id, Interface[], Dict{Symbol,Interface}())
        addNode!(currentGraph(), self)

        return self
    end
end

@testset "Edge" begin
    # Edge should couple interfaces
    FactorGraph()

    a = Interface(MockNode())
    b = Interface(MockNode())
    edge = Edge(Variable(), a, b)

    @test ===(edge.a, a)
    @test ===(edge.b, b)
    @test ===(a.edge, edge)
    @test ===(b.edge, edge)
    @test ===(a.partner, b)
    @test ===(b.partner, a)
end

@testset "disconnect!" begin
    FactorGraph()
    a = Interface(MockNode())
    b = Interface(MockNode())
    edge = Edge(Variable(), a, b)

    # disconnect! from 'a' side should decouple interfaces
    disconnect!(edge, a)
    @test ===(edge.a, b)
    @test edge.b == nothing
    @test a.edge == nothing
    @test ===(b.edge, edge)
    @test a.partner == nothing
    @test b.partner == nothing

    FactorGraph()
    a = Interface(MockNode())
    b = Interface(MockNode())
    edge = Edge(Variable(), a, b)

    # disconnect! from 'b' side should decouple interfaces
    disconnect!(edge, b)
    @test ===(edge.a, a)
    @test edge.b == nothing
    @test ===(a.edge, edge)
    @test b.edge == nothing
    @test a.partner == nothing
    @test b.partner == nothing
end

end # module


text:
# Nondivergent velocity field tangent to the sphere
# Reference:
#     Lauritzen et. al. GMD 2012
#     Eq. (18) and (19)
function W(θϕ)
  θ,ϕ = θϕ
  u = 2*sin(θ)^2*sin(2*ϕ) + 2π/5*cos(ϕ)
  v = 2*sin(2*θ)*cos(ϕ)
  spherical_to_cartesian_matrix(VectorValue(θ,ϕ,1.0))⋅VectorValue(u,v,0)
end

"""
Spherical divergence on unit sphere applied to perp(W)
"""
function divM_perpW(θϕ)
  n_times_W=(θϕ)->(W(θϕ)×normal_unit_sphere(θϕ))
  f=divergence_unit_sphere(n_times_W)
  f(θϕ)
end

function compute_error_weak_div_perp(model,order,degree,ls=BackslashSolver())
  # Setup geometry
  Ω=Triangulation(model)
  dΩ=Measure(Ω,degree)
  n=get_normal_vector(Ω)

  # Setup H(div) spaces
  reffe_rt = ReferenceFE(raviart_thomas, Float64, order)
  V = FESpace(model, reffe_rt ; conformity=:HDiv)
  U = TrialFESpace(V)

  # H1 spaces
  reffe_lgn = ReferenceFE(lagrangian, Float64, order+1)
  S = FESpace(model, reffe_lgn; conformity=:H1)
  R = TrialFESpace(S)

  # Project the analytical velocity field W onto the H(div) space
  a1(u,v) = ∫(v⋅u)dΩ
  b1(v)   = ∫(v⋅(W∘xyz2θϕ))dΩ
  op      = AffineFEOperator(a1,b1,U,V)
  fels    = LinearFESolver(ls)
  wh      = solve(fels,op)

  e = (W∘xyz2θϕ)-wh
  println(sqrt(sum(∫(e⋅e)dΩ)))

  # Compute weak div_perp operator
  a2(u,v) = ∫(v*u)dΩ
  b2(v)   = ∫(⟂(∇(v),n)⋅(-wh) )dΩ
  op      = AffineFEOperator(a2,b2,R,S)
  fels    = LinearFESolver(ls)
  divwh   = solve(fels,op)
  e = divwh-divM_perpW∘xyz2θϕ

  sqrt(sum(∫(e*e)dΩ))
end

#  plotd=plot([ahs0,bihs0,biqhs0],[ak0errors,bik0errors,biqk0errors],
#  xaxis=:log, yaxis=:log,
#  label=["k=0 analytical map" "k=0 bilinear map" "k=0 biquadratic map"],
#  shape=:auto,
#  xlabel="h",ylabel="L2 error norm", legend=:bottomright)

#  savefig(plotd,"L2_error_weak_div_perp.png")

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
using Catalyst
using DiffEqJump

rn = @reaction_network begin
    (γA, θA),   Da + A ↔ Dpa
    (γR, θR),   Dr + A ↔ Dpr
    αA,         Da → Da + Ma
    αR,         Dr → Dr + Mr
    αpA,        Dpa → Dpa + Ma
    αpR,        Dpr → Dpr + Mr
    βA,         Ma → Ma + A
    βR,         Mr → Mr + R
    γC,         A + R → C
    δA,         C → R
    δMA,        Ma → ∅
    δMR,        Mr → ∅
    δA,         A → ∅
    δR,         R → ∅
end αA αpA αR αpR βA βR δMA δMR δA δR γA γR γC θA θR

p = [50., 500., 0.01, 50., 50., 5., 10., 0.5, 1., 0.2, 1., 1., 2., 50., 100.]
u0 = [1, 0, 0, 1, 0, 0, 0, 0, 0]
tspan = (0., 200.)
prob_discrete = DiscreteProblem(rn, u0, tspan, p)
prob_jump = JumpProblem(rn, prob_discrete, Direct(), save_positions=(false, false))
for i = 1:100
    sol = solve(prob_jump, SSAStepper(), saveat=1.)
end

max_stars_repo_name: Armavica/rebop

max_issues_repo_name: Armavica/bebop

max_forks_repo_name: Armavica/rebop


text:
using Base.Test
using GLM

dobson = DataFrame(counts=[18.,17,15,20,10,20,25,13,12], outcome=gl(3,1,9), treatment=gl(3,3))
gm1 = glm(:(counts ~ outcome + treatment), dobson, Poisson());
@test_approx_eq deviance(gm1) 5.12914107700115

## Example from http://www.ats.ucla.edu/stat/r/dae/logit.htm
nm = download("http://www.ats.ucla.edu/stat/data/binary.csv", tempname())
df = within(readtable(nm),:(rank=compact(PooledDataArray(rank))))
rm(nm)

rteps = sqrt(eps())

gm2 = glm(:(admit ~ gre + gpa + rank), df, Binomial());
@test_approx_eq deviance(gm2) 458.5174924758994
@test sumsqdiff(coef(gm2), [-3.9899786606380734,0.0022644256521549043,0.8040374535155766,-0.6754428594116577,-1.3402038117481079,-1.5514636444657492]) < rteps

gm3 = glm(:(admit ~ gre + gpa + rank), df, Binomial(), ProbitLink());
@test_approx_eq deviance(gm3) 458.4131713833386
@test sumsqdiff(coef(gm3), [-2.3867922998680786,0.0013755394922972369,0.47772908362647015,-0.4154125854823675,-0.8121458010130356,-0.9359047862425298]) < rteps

gm4 = glm(:(admit ~ gre + gpa + rank), df, Binomial(), CauchitLink());
@test_approx_eq deviance(gm4) 459.3401112751141

gm5 = glm(:(admit ~ gre + gpa + rank), df, Binomial(), CloglogLink());
@test_approx_eq deviance(gm5) 458.89439629612616


text:
@cdef :_PyObject_New PyPtr (PyPtr,)
@cdef :PyObject_ClearWeakRefs Cvoid (PyPtr,)
@cdef :PyObject_HasAttrString Cint (PyPtr, Cstring)
@cdef :PyObject_HasAttr Cint (PyPtr, PyPtr)
@cdef :PyObject_GetAttrString PyPtr (PyPtr, Cstring)
@cdef :PyObject_GetAttr PyPtr (PyPtr, PyPtr)
@cdef :PyObject_GenericGetAttr PyPtr (PyPtr, PyPtr)
@cdef :PyObject_SetAttrString Cint (PyPtr, Cstring, PyPtr)
@cdef :PyObject_SetAttr Cint (PyPtr, PyPtr, PyPtr)
@cdef :PyObject_GenericSetAttr Cint (PyPtr, PyPtr, PyPtr)
@cdef :PyObject_DelAttrString Cint (PyPtr, Cstring)
@cdef :PyObject_DelAttr Cint (PyPtr, PyPtr)
@cdef :PyObject_RichCompare PyPtr (PyPtr, PyPtr, Cint)
@cdef :PyObject_RichCompareBool Cint (PyPtr, PyPtr, Cint)
@cdef :PyObject_Repr PyPtr (PyPtr,)
@cdef :PyObject_ASCII PyPtr (PyPtr,)
@cdef :PyObject_Str PyPtr (PyPtr,)
@cdef :PyObject_Bytes PyPtr (PyPtr,)
@cdef :PyObject_IsSubclass Cint (PyPtr, PyPtr)
@cdef :PyObject_IsInstance Cint (PyPtr, PyPtr)
@cdef :PyObject_Hash Py_hash_t (PyPtr,)
@cdef :PyObject_IsTrue Cint (PyPtr,)
@cdef :PyObject_Length Py_ssize_t (PyPtr,)
@cdef :PyObject_GetItem PyPtr (PyPtr, PyPtr)
@cdef :PyObject_SetItem Cint (PyPtr, PyPtr, PyPtr)
@cdef :PyObject_DelItem Cint (PyPtr, PyPtr)
@cdef :PyObject_Dir PyPtr (PyPtr,)
@cdef :PyObject_GetIter PyPtr (PyPtr,)
@cdef :PyObject_Call PyPtr (PyPtr, PyPtr, PyPtr)
@cdef :PyObject_CallObject PyPtr (PyPtr, PyPtr)

const PyObject_Type__ref = Ref(PyNULL)
PyObject_Type() = pyglobal(PyObject_Type__ref, :PyBaseObject_Type)

PyObject_From(x::PyObjectRef) = (Py_IncRef(x.ptr); x.ptr)
PyObject_From(x::Nothing) = PyNone_New()
PyObject_From(x::Missing) = PyNone_New()
PyObject_From(x::Bool) = PyBool_From(x)
PyObject_From(
    x::Union{Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128,BigInt},
) = PyLong_From(x)
PyObject_From(
    x::Rational{
        <:Union{Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128,BigInt},
    },
) = PyFraction_From(x)
PyObject_From(x::Union{Float16,Float32,Float64}) = PyFloat_From(x)
PyObject_From(x::Complex{<:Union{Float16,Float32,Float64}}) = PyComplex_From(x)
PyObject_From(x::Union{String,SubString{String}}) = PyUnicode_From(x)
PyObject_From(x::Char) = PyUnicode_From(string(x))
PyObject_From(x::Tuple) = PyTuple_From(x)
PyObject_From(
    x::AbstractRange{
        <:Union{
            Bool,
            Int8,
            Int16,
            Int32,
            Int64,
            Int128,
            UInt8,
            UInt16,
            UInt32,
            UInt64,
            UInt128,
            BigInt,
        },
    },
) = PyRange_From(x)
PyObject_From(x::DateTime) = PyDateTime_From(x)
PyObject_From(x::Date) = PyDate_From(x)
PyObject_From(x::Time) = PyTime_From(x)
PyObject_From(x::Union{Second,Millisecond,Microsecond,Nanosecond}) = PyTimeDelta_From(x)
PyObject_From(x) =
    if ispyref(x)
        GC.@preserve x begin
            ptr = pyptr(x)
            Py_IncRef(ptr)
            ptr
        end
    else
        PyJuliaValue_From(x)
        # PyErr_SetString(PyExc_TypeError(), "Cannot convert this Julia '$(typeof(x))' to a Python object.")
        # PyNULL
    end

function PyObject_CallArgs(
    f,
    args::Tuple,
    kwargs::Union{
        Nothing,
        NamedTuple,
        Base.Iterators.Pairs{Symbol},
        Base.Iterators.Pairs{Union{}},
    } = nothing,
)
    if kwargs !== nothing && !isempty(kwargs)
        argso = PyTuple_From(args)
        isnull(argso) && return PyNULL
        kwargso = PyDict_From(kwargs)
        isnull(kwargso) && (Py_DecRef(argso); return PyNULL)
        r = PyObject_Call(f, argso, kwargso)
        Py_DecRef(argso)
        Py_DecRef(kwargso)
        r
    elseif !isempty(args)
        argso = PyTuple_From(Tuple(args))
        isnull(argso) && return PyNULL
        r = PyObject_CallObject(f, argso)
        Py_DecRef(argso)
        return r
    else
        PyObject_CallObject(f, C_NULL)
    end
end

PyObject_CallNice(f, args...; kwargs...) = PyObject_CallArgs(f, args, kwargs)

"""
Mapping of Julia types to mappings of Python types to vectors of compiled functions implementing the conversion.
"""
const TRYCONVERT_COMPILED_RULES = Dict{Type,Dict{PyPtr,Vector{Ptr{Cvoid}}}}()

const TRYCONVERT_COMPILED_RULES_CACHE = Dict{Any,PyPtr}()

"""
Mapping of type names to lists of rules.

A rule is a triple `(priority, type, impl)`.

Rules are applied in priority order, then in MRO order, then in the order in this list. You should currently only use the following values:
- `0` (the default) for most rules.
- `100` for "canonical" rules which are the preferred way to convert a type, e.g. `float` to `Float64`.
- `-100` for "last ditch" rules once all other options have been exhausted (currently used for anything to `PyObject` or `PyRef`).

The rule is applied only when the `type` has non-trivial intersection with the target type `T`.

Finally `impl` is the function implementing the conversion.
Its signature is `(o, T, S)` where `o` is the `PyPtr` object to convert, `T` is the desired type, and `S = typeintersect(T, type)`.
On success it returns `putresult(x)` where `x::T` is the converted value (this stores the result and returns `1`).
If conversion was not possible, returns `0` (indicating we move on to the next rule in the list).
On error, returns `-1`.
"""
const TRYCONVERT_RULES = Dict{String,Vector{Tuple{Int,Type,Any}}}()

"""
List of niladic functions returning a pointer to a type. We always check for subtypes of these types in TryConvert.

Can also return NULL. Without an error set, this indicates the type is not loaded (e.g. its containing module is not loaded) and therefore is skipped over.
"""
const TRYCONVERT_EXTRATYPES = Vector{Any}()

@generated PyObject_TryConvert_CompiledRules(::Type{T}) where {T} =
    get!(Dict{PyPtr,Vector{Ptr{Cvoid}}}, TRYCONVERT_COMPILED_RULES, T)

PyObject_TryConvert_Rules(n::String) =
    get!(Vector{Tuple{Int,Type,Function}}, TRYCONVERT_RULES, n)
PyObject_TryConvert_AddRule(n::String, T, rule, priority = 0) =
    push!(PyObject_TryConvert_Rules(n), (priority, T, rule))
PyObject_TryConvert_AddRules(n::String, xs) =
    for x in xs
        PyObject_TryConvert_AddRule(n, x...)
    end
PyObject_TryConvert_AddExtraType(tfunc) = push!(TRYCONVERT_EXTRATYPES, tfunc)
PyObject_TryConvert_AddExtraTypes(xs) =
    for x in xs
        PyObject_TryConvert_AddExtraType(x)
    end

PyObject_TryConvert_CompileRule(::Type{T}, t::PyPtr) where {T} = begin

    ### STAGE 1: Get a list of supertypes to consider.
    #
    # This includes the MRO of t, but also the MRO of any extra type (such as an ABC)
    # that we explicitly check for because it may not appear the the MRO of t.

    # MRO of t
    tmro = PyType_MROAsVector(t)
    # find the "MROs" of all base types we are considering
    basetypes = PyPtr[t]
    basemros = Vector{PyPtr}[tmro]
    for xtf in TRYCONVERT_EXTRATYPES
        xt = xtf()
        isnull(xt) && (PyErr_IsSet() ? (return PYERR()) : continue)
        xb = PyNULL
        for b in tmro
            r = PyObject_IsSubclass(b, xt)
            ism1(r) && return PYERR()
            r != 0 && (xb = b)
        end
        if xb != PyNULL
            push!(basetypes, xt)
            xmro = PyType_MROAsVector(xt)
            push!(basemros, xb in xmro ? xmro : [xb; xmro])
        end
    end
    for b in basetypes[2:end]
        push!(basemros, [b])
    end
    # merge the MROs
    # this is a port of the merge() function at the bottom of:
    # https://www.python.org/download/releases/2.3/mro/
    alltypes = PyPtr[]
    while !isempty(basemros)
        # find the first head not contained in any tail
        ok = false
        b = PyNULL
        for bmro in basemros
            b = bmro[1]
            if all(bmro -> b ∉ bmro[2:end], basemros)
                ok = true
                break
            end
        end
        ok || error(
            "Fatal inheritence error: could not merge MROs (alltypes=$alltypes, basemros=$basemros)",
        )
        # add it to the list
        push!(alltypes, b)
        # remove it from consideration
        for bmro in basemros
            filter!(x -> x != b, bmro)
        end
        # remove empty lists
        filter!(x -> !isempty(x), basemros)
    end
    # check the original MRO is preserved
    @assert filter(x -> x in tmro, alltypes) == tmro
    # some special cases
    extranames = Dict()
    for (i, b) in reverse(collect(enumerate(alltypes)))
        if PyObject_HasAttrString(b, "__array_struct__") != 0
            push!(get!(Vector, extranames, i), "<arraystruct>")
            break
        end
    end
    for (i, b) in reverse(collect(enumerate(alltypes)))
        if PyObject_HasAttrString(b, "__array_interface__") != 0
            push!(get!(Vector, extranames, i), "<arrayinterface>")
            break
        end
    end
    for (i, b) in reverse(collect(enumerate(alltypes)))
        if PyObject_HasAttrString(b, "__array__") != 0
            push!(get!(Vector, extranames, i), "<array>")
            break
        end
    end
    for (i, b) in reverse(collect(enumerate(alltypes)))
        if PyType_CheckBuffer(b)
            push!(get!(Vector, extranames, i), "<buffer>")
            break
        end
    end
    allnames = String[]
    for (i, t) in enumerate(alltypes)
        n = PyType_FullName(t)
        if n === PYERR()
            # if we cannot get the fully qualified name of the type, just skip it (it was probably dynamically generated)
            PyErr_Clear()
        else
            push!(allnames, n)
        end
        append!(allnames, get(Vector, extranames, i))
    end

    ### STAGE 2: Get a list of applicable conversion rules.
    #
    # These are the conversion rules of the types found above

    # gather rules of the form (priority, order, S, rule) from these types
    rules = Tuple{Int,Int,Type,Any}[]
    for n in allnames
        for (p, S, r) in PyObject_TryConvert_Rules(n)
            push!(rules, (p, length(rules) + 1, S, r))
        end
    end
    # sort by priority
    sort!(rules, by = x -> (-x[1], x[2]))
    # intersect S with T
    rules = [typeintersect(S, T) => r for (p, i, S, r) in rules]
    # discard rules where S is a subtype of the union of all previous S for rules with the same implementation
    # in particular this removes rules with S==Union{} and removes duplicates
    rules = [S => rule for (i, (S, rule)) in enumerate(rules) if
     !(S <: Union{[S′ for (S′, rule′) in rules[1:i-1] if rule == rule′]...})]

    @debug "PYTHON CONVERSION FOR '$(PyType_Name(t))' to '$T'" basetypes =
        map(PyType_Name, basetypes) alltypes = allnames rules = rules

    ### STAGE 3: Define and compile functions implementing these rules.

    map(rules) do x
        get!(TRYCONVERT_COMPILED_RULES_CACHE, x) do
            S, rule = x
            # make the function implementing the rule
            rulefunc = @eval (o::PyPtr) -> $rule(o, $S)
            # compile it
            crulefunc = @cfunction($rulefunc, Int, (PyPtr,))
            # cache it
            push!(CACHE, crulefunc)
            # get the function pointer
            Base.unsafe_convert(Ptr{Cvoid}, crulefunc)
        end
    end
end

const TRYCONVERT_ERR_ARRAY = Ptr{Cvoid}[]

PyObject_TryConvert(o::PyPtr, ::Type{T}) where {T} = begin
    # First try based only on the type T
    # Used mainly by wrapper types for immediate conversion.
    r = PyObject_TryConvert__initial(o, T)::Int
    r == 0 || return r

    # Try to find an appropriate conversion based on `T` and the type of `o`.
    rules = PyObject_TryConvert_CompiledRules(T)
    t = Py_Type(o)
    crules = get(rules, t, TRYCONVERT_ERR_ARRAY)
    if crules === TRYCONVERT_ERR_ARRAY
        crules = PyObject_TryConvert_CompileRule(T, t)
        crules === PYERR() && return -1
        rules[t] = crules
    end
    for crule in crules
        r = ccall(crule, Int, (PyPtr,), o)
        r == 0 || return r
    end

    # Failed to convert
    return 0
end
PyObject_TryConvert(o, ::Type{T}) where {T} =
    GC.@preserve o PyObject_TryConvert(Base.unsafe_convert(PyPtr, o), T)

PyObject_TryConvert__initial(o, ::Type{T}) where {T} = 0

PyObject_Convert(o::PyPtr, ::Type{T}) where {T} = begin
    r = PyObject_TryConvert(o, T)
    if r == 0
        PyErr_SetString(
            PyExc_TypeError(),
            "Cannot convert this '$(PyType_Name(Py_Type(o)))' to a Julia '$T'",
        )
        -1
    elseif r == -1
        -1
    else
        0
    end
end
PyObject_Convert(o, ::Type{T}) where {T} =
    GC.@preserve o PyObject_Convert(Base.unsafe_convert(PyPtr, o), T)


text:
using SparseArrays
import StatsBase

export TPNState, TPN, isenabled, isready, isvalidstate, update_clock!, max_delay, elapse!, fire!, fire_all!, fire_necessary!, run!

"""
    TPNState

Class to hold the state (=place and transition markings) of the Time Petri Net.
"""
struct TPNState{M,H}
    m::Vector{M}                # p-marking
    h::Vector{Union{Nothing,H}} # t-marking
end

# forward declare function so it can be used in constructor
function isenabled end
function update! end

"""
    TPN{M,H}

Class to hold the entire Time Petri Net.
This implementation uses an unusual parameterization inspired by [*Time and Petri Nets*, 2013][1]:

# Attributes:
- `P::Vector{Symbol}`: names of the places
- `T::Vector{Symbol}`: names of the transitions
- `R::SparseMatrixCSC{UInt}`: additional read arcs of the transitions
- `ΔF::SparseMatrixCSC{Int}`: token-changes caused by the transitions (one column per transition)
- `eft::Vector{H}`: earliest firing time of transition
- `lft::Vector{H}`: last firing time of transition
- `x₀::TPNState{M,H}`: initial state vector

[1]: https://link.springer.com/book/10.1007/978-3-642-41115-1
"""
struct TPN{M,H}
    P::Vector{Symbol}           # names of the places
    T::Vector{Symbol}           # names of the transitions
    C::SparseMatrixCSC{M}       # all dependencies of the transitions (read arcs + regular arcs)
    ΔF::SparseMatrixCSC{M}      # token-changes caused by the transitions (one column per transition)
    eft::Vector{H}              # earliest firing time of transition
    lft::Vector{H}              # last firing time of transition
    trans_probs::Vector{Float64}# probability to fire (if multiple options are available simultaneously)
    x₀::TPNState{M,H}           # initial state vector
    # internal helper variables
    may_disable::SparseMatrixCSC{Bool}
    may_enable::SparseMatrixCSC{Bool}
    is_alternative::SparseMatrixCSC{Bool}
    P_index::Dict{Symbol,Int}
    T_index::Dict{Symbol,Int}

    """
        TPN(P,T,C,ΔF,eft,lft,m₀::Union{TPNState,Vector{<:Integer}})

    Construct the `TPN`.

    # Arguments:
    - `P`: names of the places
    - `T`: names of the transitions
    - `R`: additional read arcs of the transitions
    - `ΔF`: token-changes caused by the transitions (one column per transition)
    - `eft`: earliest firing time of transition
    - `lft`: last firing time of transition
    - `m₀::Union{TPNState,Vector{<:Integer}}`: initial state vector (if type `TPNState`) or initial marking
    """
    function TPN(P,T,R,ΔF,eft,lft,trans_probs,m₀::Union{TPNState,Vector{<:Integer}})
        M = eltype(ΔF)
        H = eltype(eft)

        C = sparse(R - ΔF.*(ΔF .< 0))
        nT = length(T)
        # initial state or initial marking provided
        x₀ = isa(m₀, TPNState) ? m₀ : TPNState(Vector{M}(m₀), Vector{Union{Nothing,H}}(nothing, nT))

        may_disable = spzeros(Bool, nT, nT)
        may_enable = spzeros(Bool, nT, nT)
        is_alternative = spzeros(Bool, nT, nT)

        P_index = Dict(value=>key for (key,value) in enumerate(P))
        T_index = Dict(value=>key for (key,value) in enumerate(T))

        pn = new{M,H}(P,T,C,ΔF,eft,lft,trans_probs,x₀, may_disable, may_enable, is_alternative, P_index, T_index)

        update!(pn)
        pn
    end
end

"""
    isenabled(pn::TPN, x, t_id)

Returns `true` if the transition with ID `t_id` of the TPN `pn` is enabled in state `x` and false otherwise.

!!! info
    The notion of a transition being `enabled` in a Time Petri Net only refers to
    the place marking (i.e. tokens) and does **not** check the transition markings (i.e. clocks)!

    For checking if the transition can be actually fired, see [`isready`](@ref)
"""
Base.@propagate_inbounds function isenabled(pn::TPN, x, t_id)
    @boundscheck @assert 1 ≤ t_id ≤ length(pn.T) "Transition $(t_id) out of bounds!";
    rows = rowvals(pn.C)
    vals = nonzeros(pn.C)
    @inbounds for i in nzrange(pn.C, t_id)
        if x.m[rows[i]] < vals[i]
            return false
        end
    end
    return true
end

"""
    isready(pn::TPN, x, t_id)

Returns `true` if the transition with ID `t_id` of the TPN `pn` is ready in state `x` and false otherwise.

!!! info
    The notion of a transition being `ready` in a Time Petri Net requires both
    the place marking (i.e. tokens) and the transition markings (i.e. clocks) to allow firing!

    For only checking the place markings, see [`isenabled`](@ref)
"""
Base.@propagate_inbounds function isready(pn::TPN, x, t_id)
    return !isnothing(x.h[t_id]) && pn.eft[t_id] ≤ x.h[t_id] ≤ pn.lft[t_id]
end

"""
    isvalidstate(pn::TPN, x)

Returns `true` if the place and transition markings of the state `x` in TPN `pn` are consistent, and `false` otherwise.
"""
function isvalidstate(pn::TPN, x)
    for (t_id,h_t) in enumerate(x.h)
        en = isenabled(pn, x, t_id)
        @inbounds if (en && h_t > lft[t_id]) || (!en && !isnothing(h_t))
            return false
        end
    end
    return true
end

"""
    update_clock!(pn::TPN, x::TPNState{M,H}, t_id) where {M,H}

Sets the clock of the TPN `pn`'s transition `t_id` in state `x` to
`zero(H)` if the transition was just enabled
or to `nothing` if the transition was just disabled

!!! warning
    This modifies the state `x`!
"""
Base.@propagate_inbounds function update_clock!(pn::TPN, x::TPNState{M,H}, t_id) where {M,H}
    en = isenabled(pn, x, t_id)
    @inbounds val = x.h[t_id]
    @inbounds x.h[t_id] = if en && isnothing(val)
        zero(H)
    elseif en
        val
    else
        nothing
    end
end

"""
    update!(pn::TPN)

Updates the TPN after structural changes have been made (e.g. adding arcs).
"""
function update!(pn::TPN{M,H}) where {M,H}
    # check which transitions can disable each other
    ΔF_rows = rowvals(pn.ΔF)
    ΔF_vals = nonzeros(pn.ΔF)
    C_rows = rowvals(pn.C)
    @inbounds for t_id in eachindex(pn.T)
        # if the number of tokens is decreased, may disable other transition
        decreased = [ΔF_rows[i] for i in nzrange(pn.ΔF, t_id) if ΔF_vals[i] < 0]
        # if the number of tokens is increased, may enable other transition
        increased = [ΔF_rows[i] for i in nzrange(pn.ΔF, t_id) if ΔF_vals[i] > 0]

        for t_id′ in eachindex(pn.T)
            checked = C_rows[nzrange(pn.C, t_id′)]
            if !isdisjoint(decreased, checked)
                pn.may_disable[t_id′,t_id] = true
                # if t_id′ disables t_id or vice versa, both transactions are alternatives to each other
                pn.is_alternative[t_id′,t_id] = true
                pn.is_alternative[t_id,t_id′] = true
            end
            if !isdisjoint(increased, checked)
                pn.may_enable[t_id′,t_id] = true
            end
        end
    end

    # set clocks of all enabled transitions in initial state to zero
    @inbounds for t_id in eachindex(pn.T)
        if isenabled(pn, pn.x₀, t_id)
            pn.x₀.h[t_id] = zero(H)
        end
    end
    return nothing
end

"""
    fire!(pn::TPN, x::TPNState, t_id; check_ready=true)

Fires the transition `x` of TPN `pn` from state `x`.

If `check_ready` is `true` (the default), then ensure that the transition is ready before firing.

!!! warning
    This modifies the state `x`!
"""
Base.@propagate_inbounds function fire!(pn::TPN, x::TPNState, t_id; check_ready=true)
    if !check_ready || isready(pn, x, t_id)
        rows = rowvals(pn.ΔF)
        vals = nonzeros(pn.ΔF)
        @inbounds for i in nzrange(pn.ΔF, t_id)
            @inbounds x.m[rows[i]] += vals[i]
        end

        @inbounds for t_id′ in rowvals(pn.may_disable)[nzrange(pn.may_disable,t_id)] ∪ rowvals(pn.may_enable)[nzrange(pn.may_enable,t_id)]
            update_clock!(pn,x,t_id′)
        end
    else
        error("Transition #$(t_id) not ready!")
    end
    nothing
end

"""
    fire_all!(pn::TPN, x; must_fire=isready)

Fires all transitions `t` of TPN `pn` in state `x` for which `must_fire(pn, x, t)` returns `true`.

Returns a vector of the fired transition labels.

!!! warning
    This modifies the state `x`!
"""
function fire_all!(pn::TPN, x; must_fire=isready)
    options = Set([t_id for t_id in eachindex(pn.T) if must_fire(pn, x, t_id)])
    transitions = Symbol[]
    while !isempty(options)
        # pick a random option
        opt = rand(options)
        # find the alternative options that are currently viable (including 'opt'!)
        # and might conflict with this one
        alternatives = [val for val in rowvals(pn.is_alternative)[nzrange(pn.is_alternative,opt)] if val ∈ options]
        # find the probabilities of these alternatives
        alternatives_probs = [pn.trans_probs[t] for t in alternatives]
        # select one alternative according to its probability
        t_id = pop!(options, StatsBase.sample(alternatives, StatsBase.AnalyticWeights(alternatives_probs)))
        if length(alternatives)>1
            println("Picked $(pn.T[t_id]) out of these options: $([(pn.T[t], pn.trans_probs[t]) for t in alternatives])")
        end
        fire!(pn, x, t_id; check_ready=false)
        push!(transitions, pn.T[t_id])

        # check if transition must fire again
        if must_fire(pn, x, t_id)
            push!(options, t_id)
        end

        # check transitions that might have been disabled
        @inbounds for t_id′ in rowvals(pn.may_disable)[nzrange(pn.may_disable,t_id)]
            if !isready(pn, x, t_id′)
                delete!(options, t_id′)
            end
        end

        # check transitions that might have been enabled
        @inbounds for t_id′ in rowvals(pn.may_enable)[nzrange(pn.may_enable,t_id)]
            if must_fire(pn, x, t_id′)
                push!(options, t_id′)
            end
        end
    end
    return transitions
end

"""
    fire_necessary!(pn::TPN, x::TPNState)

Fires all transitions `t` of TPN `pn` in state `x` that are ready and cannot be delayed any further.

Returns a vector of the fired transition labels.

!!! warning
    This modifies the state `x`!
"""
function fire_necessary!(pn::TPN, x::TPNState)
    isnecessary(pn, x, t_id) =
        @inbounds isready(pn, x, t_id) && x.h[t_id] == pn.lft[t_id]
    fire_all!(pn, x; must_fire=isnecessary)
end

"""
    max_delay(pn, x::TPNState{TT}) where TT

Computes the maximum delay possible of TPN `pn` in state `x` that does not invalidate the state.
"""
function max_delay(pn, x::TPNState{TT}) where TT
    m = typemax(TT)
    for (tᵢ,lftᵢ) in zip(x.h,pn.lft)
        if !isnothing(tᵢ)
            m = min(m, lftᵢ-tᵢ)
        end
    end
    return m
end

"""
    elapse!(pn::TPN, x::TPNState, τ; check_valid=true)

Elapses the time `τ`, i.e. advances all clocks of TPN `pn` in state `x` by `τ`.
If `check_valid` is `true` (the default), then ensure that `τ` does not exceed the maximum allowed delay.

!!! warning
    This modifies the state `x`!
"""
Base.@propagate_inbounds function elapse!(pn::TPN, x::TPNState, τ; check_valid=true)
    if !check_valid || τ ≤ max_delay(pn, x)
        for (t_id,tᵢ) in enumerate(x.h)
            if !isnothing(tᵢ)
                @inbounds x.h[t_id] += τ
            end
        end
    else
        error("Time delay exceeds maximum delay ($(τ) > $(max_delay(pn,x)))!")
    end
end

"""

Runs the TPN `pn` from state `x`

!!! warning
    This modifies the state `x`!
"""
function run!(pn::TPN, x, events; tmax=first(last(events)))
    @assert issorted(events, by=first) "Events must be a sequence of (time,transition) pairs with increasing time."

    t_total = 0
    for (time, ext_t_id) in events
        while t_total < time
            delay = max_delay(pn, x)
            # elapse time until the next (external or internal event)
            if time ≤ t_total + delay
                # squeeze in our externally triggered transition
                t_total = time
                elapse!(pn, x, time-t_total)
                fire!(pn, x, ext_t_id)

                # fire all other events (to be sure)
                fire_necessary!(pn, x, ext_t_id)
            else
                # proceed to next regular event
                elapse!(pn, x, delay)
                fire_necessary!(pn, x, ext_t_id)
            end

            t_new = min(time, t_total+_max_delay)
            elapse!(pn, x, t_new-t_total)
            t_total = t_new
        end
    end

    x=copy(x)
    trace = []
    states = [copy(x)]
    for t_id in t_ids
        fire!(pn, x, t_id)
        push!(trace, pn.T[t_id])
        push!(states, copy(x))
    end
    return (;trace, states)
end

"""
    |(tpns::TPN{M,H}...) where {M,H}

Combines multiple TPNs together.

This procedure:
- merges transitions by adding arcs,
- merges earliest and last firing times by taking the more restrictive bounds,
- merges initial tokens by adding them together

"""
function Base.:|(tpns::TPN{M,H}...) where {M,H}
    all_places = union((tpn.P for tpn in tpns)...)
    all_transitions = union((tpn.T for tpn in tpns)...)
    all_R = spzeros(M, length(all_places), length(all_transitions))
    all_ΔF = spzeros(M, length(all_places), length(all_transitions))
    all_eft = zeros(H, length(all_transitions))
    all_lft = fill(typemax(H), length(all_transitions))
    all_p_trans = ones(Float64, length(all_transitions))
    all_m₀ = zeros(M, length(all_places))

    # merge transitions by adding arcs,
    # merge earliest and last firing times by taking the more restrictive bounds,
    # merge initial tokens by adding them together
    for tpn in tpns
        # map the TPN's places and transitions to the new places and transitions
        p_idx = indexin(tpn.P, all_places)
        t_idx = indexin(tpn.T, all_transitions)
        R = tpn.C + tpn.ΔF .* (tpn.ΔF .< 0)
        all_R[p_idx,t_idx] += R
        all_ΔF[p_idx,t_idx] += tpn.ΔF
        all_eft[t_idx] .= max(all_eft[t_idx], tpn.eft)
        all_lft[t_idx] .= min(all_lft[t_idx], tpn.lft)
        all_p_trans[t_idx] .= min(all_p_trans[t_idx], tpn.trans_probs)
        all_m₀[p_idx] .+= tpn.x₀.m
    end

    TPN(all_places, all_transitions, all_R, all_ΔF, all_eft, all_lft, all_p_trans, all_m₀)
end


text:
@inline function ExaTronKernels.daxpy(n::Int,da::Float64,
                       dx::oneDeviceArray{Float64,1},incx::Int,
                       dy::oneDeviceArray{Float64,1},incy::Int)
    tx = get_local_id()

    if tx <= n
        @inbounds dy[tx] = dy[tx] + da*dx[tx]
    end
    barrier()

    return
end


text:
# This file holds the definition of the functions pertaining to the Attrition
#   type.

# The functions of the Attrition type requires the State type.
requiredTypes = [ "state", "attrition" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export setAttritionRate
"""
setAttritionRate( attrScheme::Attrition,
                  rate::T )
    where T <: Real

This function sets the attrition curve of the attrition scheme `attrScheme` to a
constant value of `rate` for all terms. If the given attrition rate is 0, this
function will also set the attrition period to 1.

This function returns `nothing`.
"""
function setAttritionRate( attrScheme::Attrition, rate::T ) where T <: Real

    if ( rate < 0.0 ) || ( rate >= 1.0 )
        warn( "Attrition rate must be a percentage between 0.0 and 1.0, not making any changes." )
        return
    end  # if ( rate < 0.0 ) || ( rate > 1.0 )

    attrScheme.attrCurvePoints = [ 0.0 ]
    attrScheme.attrRates = [ rate ]

    # If the attrition rate becomes 0.0, reset the period to 1.0.
    # if rate == 0.0
    #     attrScheme.attrPeriod = 1.0
    # end  # if rate == 0.0

    computeDistPars( attrScheme )

end  # setAttritionRate( attrScheme, rate )

export setAttritionPeriod
"""
setAttritionPeriod( attrScheme::Attrition,
                    period::T )
    where T <: Real

This function sets the attrition period of the attrition scheme `attrScheme` to
`period`. If the attrition curve is a flat 0, the function warns that it has no
practical effect.

This function returns `nothing`.
"""
function setAttritionPeriod( attrScheme::Attrition, period::T ) where T <: Real

    if period <= 0.0
        warn( "Attrition period must be > 0.0, not making any changes." )
        return
    end  # if period <= 0.0

    # if attrScheme.attrRates == [ 0.0 ]
    #     warn( "Attrition rate is a flat 0.0, changing the period has no effect." )
    #     return
    # end  # if attrScheme.attrRates == [ 0.0 ]

    attrScheme.attrPeriod = period
    computeDistPars( attrScheme )

end  # setAttritionPeriod( attrScheme, period )

export setAttritionParameters
"""
setAttritionParameters( attrScheme::Attrition,
                        rate::T1,
                        period::T2 )
    where T1 <: Real where T2 <: Real

This function sets the attrition rate of attrition scheme `attrScheme` to a
fixed attrition rate `rate` per period of length `period`.

This function returns `nothing`.
"""
function setAttritionParameters( attrScheme::Attrition, rate::T1, period::T2 ) where T1 <: Real where T2 <: Real

    setAttritionRate( attrScheme, rate )

    # if rate != 0.0
        setAttritionPeriod( attrScheme, period )
    # end  # if rate != 0.0

end  # setAttritionParameters( attrScheme, rate, period )

"""
setAttritionParameters( mpSim::ManpowerSimulation,
                        rate::T1,
                        period::T2 )
    where T1 <: Real where T2 <: Real

This function sets the attrition rate of attrition scheme in the manpower
simulation `mpSim` to a fixed attrition rate `rate` per period of length
`period`.

This function returns `nothing`.
"""
function setAttritionParameters( mpSim::ManpowerSimulation, rate::T1,
    period::T2 ) where T1 <: Real where T2 <: Real

    if rate == 0.0
        mpSim.defaultAttritionScheme = nothing
    elseif mpSim.defaultAttritionScheme === nothing
        mpSim.defaultAttritionScheme = Attrition( rate, period )
    else
        setAttritionParameters( mpSim.defaultAttritionScheme, rate, period )
    end  # if rate == 0.0

end  # setAttritionParameters( mpSim, rate, period )

export setAttritionCurve
"""
setAttritionCurve( attrScheme::Attrition,
                   curve::Array{Float64, 2} )

This function sets the attrition curve of attrition scheme `attrScheme` to
`curve`. The attrition curve is passed as a 2-dimensional `Array{Float64, 2}`
with 2 columns. The first column has the time a person exists in the simulation,
and the second column has the attrition rates per period from the time specified
by that key to the next.

This function returns `nothing`. If the given array does not have exactly 2
columns, or if there are duplicate entries in the first column, this function
does not make any changes.

This function ignores negative terms and non sensical attrition rates. If this
results in an empty set of eligible curve points, the attrition rate is set to a
constant 0 instead. The first eligible node has its term set to 0.
"""
function setAttritionCurve( attrScheme::Attrition, curve::Array{Float64, 2} )

    # Check for correct dimensions.
    if size( curve )[ 2 ] != 2
        warn( "Invalid array given to define attrition curve, not making any changes." )
        return
    end  # if size( curve )[ 2 ] != 2

    # Check for duplicate terms.
    if length( curve[ :, 1 ] ) != length( unique( curve[ :, 1 ] ) )
        warn( "Duplicate entries in the terms of the attrition curve, not making any changes." )
        return
    end  # if length( curve[ :, 1 ] ) != length( unique( curve[ :, 1 ] ) )

    # Remove curve points with non sensical attrition rates.
    tmpCurve = curve[ map( ii -> 0 <= curve[ ii, 2 ] < 1,
        eachindex( curve[ :, 1 ] ) ), : ]

    # Flat 0 rate if no more nodes are eligible.
    if isempty( tmpCurve )
        setAttritionRate( attrScheme, 0 )
        return
    end  # if isempty( tmpNodes )

    # Remove curve points with negative terms.
    tmpCurve = sortrows( tmpCurve, by = x -> x[ 1 ] )
    lastNegIndex = findlast( tmpCurve[ :, 1 ] .<= 0 )
    tmpCurve = tmpCurve[ lastNegIndex:end, : ]

    # Filter out consecutive nodes with the same attrition rate.
    tmpCurve[ 1, 1 ] = 0.0
    isDiffFromPrevious = vcat( true,
        tmpCurve[ 2:end, 2 ] .!= tmpCurve[ 1:( end - 1 ), 2 ] )

    # Update the attrition scheme.
    attrScheme.attrCurvePoints = tmpCurve[ isDiffFromPrevious, 1 ]
    attrScheme.attrRates = tmpCurve[ isDiffFromPrevious, 2 ]
    computeDistPars( attrScheme )

end  # setAttritionCurve( attrScheme, curve )

"""
setAttritionCurve( attrScheme::Attrition,
                   curve::Dict{Float64, Float64} )

This function sets the attrition curve of attrition scheme `attrScheme` to
`curve`. The attrition curve is passed as a `Dict{Float64, Float64}` with the
keys as the time a person exists in the simulation, and the values as the
attrition rates per period from the time specified by that key to the next.

This function returns `nothing`.

This function ignores negative keys and non sensical attrition rates. If this
results in an empty set of eligible curve points, the attrition rate is set to a
constant 0 instead. The first eligible node has its term set to 0.
"""
function setAttritionCurve( attrScheme::Attrition,
    curve::Dict{Float64, Float64} )

    terms = collect( keys( curve ) )
    rates = map( term -> curve[ term ], terms )
    setAttritionCurve( attrScheme, hcat( terms, collect( rates ) ) )

end  # setAttritionCurve( attrScheme, curve )

export addAttritionScheme!
"""
```
addAttritionScheme!( attrScheme::Attrition,
                     mpSim::ManpowerSimulation )
```
This function adds the attrition scheme `attrScheme` to the manpower simulation
`mpSim`. If the attrition scheme has the name `default`, the function overwrites
the default scheme of the manpower simulation (and warns about this), and if the
attrition scheme has a name that is already used, the other scheme with this
name gets overwritten.

This function returns `nothing`.
"""
function addAttritionScheme!( attrScheme::Attrition,
    mpSim::ManpowerSimulation )::Void

    # If the scheme's name is "default", overwrite the default scheme.
    if attrScheme.name == "default"
        warn( "Overwriting default attrition scheme." )
        mpSim.defaultAttritionScheme = attrScheme
        return
    end  # if attrScheme.name == "default"

    # If an attrition scheme with that name exists, overwrite it.
    if haskey( mpSim.attritionSchemes, attrScheme.name )
        warn( "Overwriting attrition scheme '$(attrScheme.name)'.")
    end  # if haskey( mpSim.attritionSchemes, attrScheme.name )

    mpSim.attritionSchemes[ attrScheme.name ] = attrScheme
    return

end  # addAttritionScheme!( attrScheme, mpSim )

export clearAttritionSchemes!
"""
```
clearAttritionSchemes!( mpSim::ManpowerSimulation )
```
This function clears all attrition schemes, except the default, from the
manpower simulation `mpSim`.

This function returns `nothing`.
"""
function clearAttritionSchemes!( mpSim::ManpowerSimulation )::Void

    empty!( mpSim.attritionSchemes )
    return

end  # clearAttritionSchemes!( mpSim )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

"""
```
readAttrition( attrCat::XLSX.Worksheet,
               sLine::Int )
```
This function reads the Excel sheet `attrCat`, starting in the row with number
`sLine`, to read the parameters of the attritition scheme.

This function returns an object of type `Attrition`, the attrition scheme.
"""
function readAttrition( attrCat::XLSX.Worksheet, sLine::Int )::Attrition

    attrName = string( attrCat[ "A$sLine" ] )

    # Check if it's the default attrition scheme.
    if lowercase( attrName ) == "default"
        attrName = "default"
    end  # if lowercase( attrName ) == "default"

    newAttrScheme = Attrition( attrName )
    attrPeriod = attrCat[ "B$sLine" ]
    nPoints =  attrCat[ "C$sLine" ]
    attrCurve = zeros( nPoints, 2 )

    for ii in 1:nPoints
        attrCurve[ ii, : ] = [ attrCat[ XLSX.CellRef( sLine, 2 + 2 * ii ) ],
            attrCat[ XLSX.CellRef( sLine, 3 + 2 * ii ) ] ]
    end  # for ii in 1:nPoints

    attrCurve[ :, 1 ] *= attrPeriod
    setAttritionCurve( newAttrScheme, attrCurve )
    setAttritionPeriod( newAttrScheme, attrPeriod )
    return newAttrScheme

end  # readAttrition( attrCat, sLine )

"""
```
computeDistPars( attrScheme::Attrition )
```
This function pre-computes some reoccurring parameters of the distribution of
the time to attrition of the attrition scheme `attrScheme`.

This function returns `nothing`.
"""
function computeDistPars( attrScheme::Attrition )::Void

    attrScheme.lambdas = - log.( 1 - attrScheme.attrRates ) /
        attrScheme.attrPeriod
    alpha = vcat( 1.0, exp.( - ( attrScheme.lambdas[ 1:(end - 1) ] -
        attrScheme.lambdas[ 2:end ] ) .* attrScheme.attrCurvePoints[ 2:end ] ) )
    attrScheme.betas = cumprod( alpha )
    attrScheme.gammas = attrScheme.betas .*
        exp.( - attrScheme.lambdas .* attrScheme.attrCurvePoints )
    return

end  # computeDistPars( attrScheme )

"""
```
determineAttritionScheme( stateList::Vector{String},
                          mpSim::ManpowerSimulation )
```
This function determines the appropriate attrition scheme to use for a personnel
member belonging to the states in `statList`. The available attrition schemes
are found in the manpower simulation `mpSim`.

This function returns an object of type `Attrition` (an attrition scheme). If
multiple schemes qualify, the first is returned; if none qualify, the default
attrition scheme is returned.
"""
function determineAttritionScheme( stateList::Vector{State},
    mpSim::ManpowerSimulation )::Attrition

    # XXX Check states and determine appropriate schemes.
    attrList = map( state -> state.attrScheme, stateList )
    filter!( attrScheme -> attrScheme !== mpSim.defaultAttritionScheme,
        attrList )

    return isempty( attrList ) ? mpSim.defaultAttritionScheme : attrList[ 1 ]

end  # determineAttritionScheme( stateList, mpSim )

"""
```
generateTimeOfAttrition( attrScheme::Attrition,
                         nowTime::Float64 )
```
This function generates a time to attrition for the attrition scheme
`attrScheme`, starting from the current time.

This function returns a `Float64`, the time to attrition, or a `String` equal to
`"NULL"` if the time to attrition is infinite.
"""
function generateTimeOfAttrition( attrScheme::Attrition, nowTime::Float64 )

    u = rand()
    ii = findlast( 1 - attrScheme.gammas .<= u )

    if ( ii == length( attrScheme.attrCurvePoints ) ) &&
        ( attrScheme.lambdas[ end ] == 0 )
        return "NULL"
    end  # if ( ii == length( attrScheme.attrCurvePoints ) ) && ...

    return nowTime -
        log( ( 1 - u ) / attrScheme.betas[ ii ] ) / attrScheme.lambdas[ ii ]

end  # generateTimeOfAttrition( attrScheme, nowTime )

"""
```
updateTimeToAttrition( stateList::Dict{String, Vector{State}},
                       mpSim::ManpowerSimulation )
```
This function updates the time to attrition for all the personnel members with
ids in the key list of `stateList` (the argument is a list of the states the
personnel member is in), using the attrition schemes defined in the manpower
simulation `mpSim`.

This function returns `nothing`.
"""
function updateTimeToAttrition( stateList::Dict{String, Vector{State}},
    mpSim::ManpowerSimulation )::Void

    if isempty( stateList )
        return
    end  # if isempty( stateList )

    queryCmd = "SELECT $(mpSim.idKey), attritionScheme
        FROM $(mpSim.personnelDBname)
        WHERE $(mpSim.idKey) IN ('$(join( collect( keys( stateList ) ), "', '" ))')"
    attrList = SQLite.query( mpSim.simDB, queryCmd )
    idList = Vector{String}()
    newAttrTimeList = Vector{Float64}()

    for ii in eachindex( attrList[ Symbol( mpSim.idKey ) ] )
        id = attrList[ Symbol( mpSim.idKey ) ][ ii ]
        attrScheme = determineAttritionScheme( stateList[ id ], mpSim )

        # Generate a new time to attrition
        if attrList[ :attritionScheme ][ ii ] != attrScheme.name
            timeToAttr = generateTimeOfAttrition( attrScheme, now( mpSim ) )
            push!( idList, id )
            push!( newAttrTimeList,
                isa( timeToAttr, String ) ? +Inf : timeToAttr )

            # Update database
            updateCmd = "UPDATE $(mpSim.personnelDBname)
                SET attritionScheme = '$(attrScheme.name)',
                    expectedAttritionTime = $timeToAttr
                WHERE $(mpSim.idKey) IS '$id'"
            SQLite.execute!( mpSim.simDB, updateCmd )
        end  # if attrSchemeList[ id ] != attrScheme.name
    end  # for ii in length( attrList[ Symbol( mpSim.idKey ) ] )

    if isempty( idList )
        return
    end  # if isempty( idList )

    # Check if any of the new attrition times occur before the next attrition
    #   process check.
    newProcessCheckTime = ceil( now( mpSim ) / mpSim.attritionTimeSkip ) *
        mpSim.attritionTimeSkip

    for ii in eachindex( idList )
        if newAttrTimeList[ ii ] <= newProcessCheckTime
            executeAttritionProcess( mpSim.sim, idList[ ii ],
                newAttrTimeList[ ii ], mpSim )
        end  # if newAttrTimeList[ ii ] <= newProcessCheckTime
    end

end  # updateTimeToAttrition( stateList, mpSim )

#=
@resumable
attritionProcess( sim::SimJulia.Simulation,
                  id::String,
                  timeOfRetirement::T,
                  retProc::ResumableFunctions.Process,
                  mpSim::ManpowerSimulation )
    where T <: Real

This SimJulia process handles the career attrition for the person with ID `id`
in the manpower simulation `mpSim`.

The expected time of retirement for this person is `timeOfRetirement`, and the
associated retirement process is `retProc`. The argument `sim` is the
SimJulia.Simulation component of the manpower simulation, and is required by the
@process macro.
=#
@resumable function attritionProcess( sim::Simulation, id::String,
    timeOfRetirement::T, #retProc::Process, mpSim::ManpowerSimulation ) where T <: Real
    mpSim::ManpowerSimulation ) where T <: Real

    attrScheme = mpSim.defaultAttritionScheme
    timeOfEntry = now( sim )
    checkForAttrition = true
    rateIndex = 1
    attrRate = -1
    lambda = 0
    timeOfRateChange = 0
    queryCmd = "SELECT $(mpSim.idKey) FROM $(mpSim.personnelDBname)
        WHERE ( $(mpSim.idKey) IS '$id' ) AND
            ( status NOT IN ('fired', 'retirement') )"

    while checkForAttrition
        # Update attrition rate and time of next rate change if necessary.
        if attrScheme.attrRates[ rateIndex ] != attrRate
            attrRate = attrScheme.attrRates[ rateIndex ]
            lambda = - log( 1 - attrRate ) / attrScheme.attrPeriod
            timeOfRateChange = length( attrScheme.attrRates ) == rateIndex ?
                +Inf : attrScheme.attrCurvePoints[ rateIndex + 1 ] + timeOfEntry
        end  # if attrScheme.attrRates[ rateIndex ] != attrRate

        # If the attrition rate is 0, skip right ahead or terminate the
        #   attrition process, whichever is required.
        if attrRate == 0
            # Terminate if the rest of the curve is a flat 0 attrition rate.
            if length( attrScheme.attrRates ) == rateIndex
                checkForAttrition = false
            # Terminate if the time of the next rate change is
            # 1. Outside the simulation time frame, or
            # 2. Past the expected retirement time.
            elseif ( timeOfRateChange > mpSim.simLength ) ||
                ( timeOfRateChange > timeOfRetirement )
                checkForAttrition = false
            # Otherwise, skip ahead to the next attrition rate change.
            else
                @yield timeout( sim, timeOfRateChange - now( sim ),
                    priority = mpSim.phasePriorities[ :attrition ] )
                rateIndex += 1
            end  # if length( attrScheme.attrRates ) == rateIndex
        # Otherwise, perform attrition checks.
        else
            attrTime = rand( Exponential( 1 / lambda ) )
            timeOfAttr = now( sim ) + attrTime
            nextAttrPeriodStart = now( sim ) + attrScheme.attrPeriod

            # Check if attrition happens in the current attrition period.
            # 1. Attrition happens in simulation timeframe
            # 2. If yes, attrition happens in attrition period.
            # 3. If yes, attrition happens before next rate change.
            # 4. If yes, attrition happens before retirement.
            # 5. If yes, interrupt retirement process.
            if ( timeOfAttr <= mpSim.simLength ) &&
                ( attrTime <= attrScheme.attrPeriod ) &&
                ( timeOfAttr <= timeOfRateChange ) &&
                ( timeOfAttr <= timeOfRetirement )
                @yield timeout( sim, attrTime,
                    priority = mpSim.phasePriorities[ :attrition ] )

                # Make sure the person hasn't been fired.
                result = SQLite.query( mpSim.simDB, queryCmd )

                if size( result ) == ( 1, 1 )
                    retirePerson( mpSim, id, "attrition" )
                end  # if size( result ) == ( 1, 1 )

                checkForAttrition = false

                # The interrupt is only sensible if the retirement process (still)
                #   exists. This can be detected by testing for finite exptected
                #   retirement time.
                # if timeOfRetirement < +Inf
                #     interrupt( retProc )
                # end  # if timeOfRetirement < +Inf
            end  # if ( timeOfAttr <= mpSim.simLength ) && ...
        end  # if attrRate == 0

        # Check if it's sensible to start a next attrition period.
        # 1. Attrition hasn't happened yet.
        # 2. If yes, next attrition period starts in simulation timeframe.
        # 3. If yes, next attrition period starts before retirement.
        # 4. If yes, wait for next attrition period.
        if checkForAttrition && ( nextAttrPeriodStart <= mpSim.simLength ) &&
            ( nextAttrPeriodStart < timeOfRetirement )
            tmpWaitTime = min( attrScheme.attrPeriod,
                timeOfRateChange - now( sim ) )

            @yield timeout( sim, tmpWaitTime,
                priority = mpSim.phasePriorities[ :attrition ] )

            # Make sure the person hasn't been fired.
            result = SQLite.query( mpSim.simDB, queryCmd )

            if size( result ) != ( 1, 1 )
                checkForAttrition = false
            end  # if size( result ) != ( 1, 1 )

            # Prepare to update attrition rate if necessary.
            if now( sim ) == timeOfRateChange
                rateIndex += 1
            end
        else
            checkForAttrition = false
        end  # if checkForAttrition && ...
    end  # while checkForAttrition

end  # attritionProcess( sim, id, timeOfRetirement, mpSim )

@resumable function checkAttritionProcess( sim::Simulation, mpSim::ManpowerSimulation )

    processTime = Dates.Millisecond( 0 )
    tStart = now()
    timeOfNextCheck = 0.0
    priority = mpSim.phasePriorities[ :attrCheck ]
    queryCmd = "SELECT $(mpSim.idKey), expectedAttritionTime
        FROM $(mpSim.personnelDBname)
        WHERE status NOT IN ('retirement', 'attrition', 'fired')
            AND expectedAttritionTime <= "

    while now( sim ) < mpSim.simLength
        processTime += now() - tStart

        @yield( timeout( sim, timeOfNextCheck - now( sim ),
            priority = priority ) )

        tStart = now()
        timeOfNextCheck += mpSim.attritionTimeSkip
        timeOfNextCheck = timeOfNextCheck > mpSim.simLength ? mpSim.simLength :
            timeOfNextCheck
        queryResult = SQLite.query( mpSim.simDB,
            queryCmd * string( timeOfNextCheck ) )
        processTime += now() - tStart

        for ii in 1:length( queryResult[ 1 ] )
            @process executeAttritionProcess( sim, queryResult[ 1 ][ ii ],
                queryResult[ 2 ][ ii ], mpSim )
        end  # for ii in 1:length( queryResult[ 1 ] )

        tStart = now()
    end  # while now( sim ) < mpSim.simLength

    processTime += now() - tStart
    println( "Attrition check process took $(processTime.value / 1000) seconds." )

end  # checkAttritionProcess( sim, mpSim )

# This process executes the attrition if it occurs in the next period.
@resumable function executeAttritionProcess( sim::Simulation, id::String,
    timeOfAttrition::Float64, mpSim::ManpowerSimulation )

    # Wait until the time attrition should take place.
    @yield( timeout( sim, timeOfAttrition - now( sim ),
        priority = mpSim.phasePriorities[ :attrition ] ) )

    tStart = now()
    # Get the personnel record.
    queryCmd = "SELECT expectedAttritionTime, status
        FROM $(mpSim.personnelDBname)
        WHERE $(mpSim.idKey) IS '$id'"
    persRecord = SQLite.query( mpSim.simDB, queryCmd )

    # Only perform the attrition if the person is still active and the time of
    #   attrition hasn't changed.
    if ( persRecord[ :status ][ 1 ] ∉
        [ "attrition", "retirement", "fired" ] ) &&
        ( persRecord[ :expectedAttritionTime ][ 1 ] == timeOfAttrition )
        retirePerson( mpSim, id, "attrition" )
    end  # if ( persRecord[ :status ][ 1 ] ∉ ...

    mpSim.attrExecTimeElapsed += now() - tStart

end  # executeAttritionProcess( sim, id, timeOfAttrition, mpSim )

function Base.show( io::IO, attrScheme::Attrition )

    print( io, "  Attrition scheme '$(attrScheme.name)'" )

    if attrScheme.attrRates == [ 0.0 ]
        print( io, "\n    No attrition." )
        return
    end  # if attrScheme.attrRates == [ 0.0 ]

    print( io, "\n    Attrition period: $(attrScheme.attrPeriod)" )

    if length( attrScheme.attrRates ) == 1
        print( io, "\n    Attrition rate: $(attrScheme.attrRates[ 1 ] * 100)%" )
    else
        print( io, "\n    Attrition curve" )

        for ii in eachindex( attrScheme.attrRates )
            print( io, "\n      term: $(attrScheme.attrCurvePoints[ ii ]);" )
            print( io, "   rate: $(attrScheme.attrRates[ ii ] * 100)%" )
        end  # for ii in eachindex( attrScheme.attrRates )
    end  # if length( attrScheme.attrRates ) == 1

end  # Base.show( io, retScheme )


text:
using Games: pure_nash
using QuantEcon: MVNSampler
using Games.Generators

include("test_bimatrix_generators.jl")


text:
using OpenStreetMapIO, Test

@testset "Testing `pbf` file" begin
    @time osmdata = OpenStreetMapIO.readpbf("data/map.pbf")

    @testset "Testing Node" begin
        node = osmdata.nodes[1675598406]

        @test typeof(node) === Node
        @test node.latlon === LatLon(54.2619665, 9.9854149)
        @test length(node.tags) === 5
        @test node.tags[Symbol("addr:country")] === "DE"
    end

     @testset "Testing Way" begin
        way = osmdata.ways[889648159]

        @test typeof(way) === Way
        @test length(way.refs) === 56
        @test way.refs[23] === 1276389426
        @test length(way.tags) === 2
        @test way.tags[Symbol("wetland")] === "wet_meadow"
    end

    @testset "Testing Relation" begin
        relation = osmdata.relations[12475101]

        @test typeof(relation) === Relation
        @test length(relation.refs) === 136
        @test relation.refs[23] === 324374700
        @test length(relation.types) === 136
        @test relation.types[23] === :node
        @test length(relation.roles) === 136
        @test relation.roles[23] === "platform"
        @test length(relation.tags) === 8
        @test relation.tags[Symbol("from")] === "Bordesholm Bahnhof"
    end

end

@testset "Testing `osm` file" begin
    @time osmdata = OpenStreetMapIO.readosm("data/map.osm")

    @testset "Testing Node" begin
        node = osmdata.nodes[1675598406]

        @test typeof(node) === Node
        @test node.latlon === LatLon(54.2619665, 9.9854149)
        @test length(node.tags) === 5
        @test node.tags[Symbol("addr:country")] === "DE"
    end

    @testset "Testing Way" begin
        way = osmdata.ways[889648159]

        @test typeof(way) === Way
        @test length(way.refs) === 56
        @test way.refs[23] === 1276389426
        @test length(way.tags) === 2
        @test way.tags[Symbol(:wetland)] === "wet_meadow"
    end

    @testset "Testing Relation" begin
        relation = osmdata.relations[12475101]

        @test typeof(relation) === Relation
        @test length(relation.refs) === 136
        @test relation.refs[23] === 324374700
        @test length(relation.types) === 136
        @test relation.types[23] === :node
        @test length(relation.roles) === 136
        @test relation.roles[23] === "platform"
        @test length(relation.tags) === 8
        @test relation.tags[Symbol("from")] === "Bordesholm Bahnhof"
    end
end


text:
module EC2

using ...AWSTools: timeout
using HTTP: HTTP
using Mocking

export instance_metadata, instance_region, instance_availability_zone

"""
    instance_metadata(path) -> Union{String,Nothing}

Retrieve AWS EC2 instance metadata as a string from the provided `path`. If no instance
metadata is available (typically due to not running within an EC2 instance) then `nothing`
will be returned. See the AWS documentation for details on what metadata is available.

https://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ec2-instance-metadata.html
"""
function instance_metadata(path::AbstractString)
    # Retrieve details about the instance:
    # https://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ec2-instance-metadata.html
    #
    # Note: When running outside of EC2 the connection to the local-link address will fail
    # with a connection timeout (ETIMEDOUT) after a 60 seconds (tested on HTTP.jl v0.8.2)
    # See: https://github.com/JuliaWeb/HTTP.jl/issues/114
    uri = HTTP.URI(scheme="http", host="169.254.169.254", path=path)

    r = timeout(5) do
        @mock HTTP.get(uri, status_exception=false)
    end

    return r !== nothing ? String(something(r).body) : r
end

"""
    instance_availability_zone() -> Union{String,Nothing}

Get the availability zone of the host if running inside of an EC2 instance. If not running
within an EC2 instance `nothing` is returned.
"""
function instance_availability_zone()
    # Get the availability zone information from the EC2 instance metadata.
    return instance_metadata("/latest/meta-data/placement/availability-zone")
end

"""
    instance_availability_zone() -> Union{String,Nothing}

Get the region of the host if executed inside of an EC2 instance. If not running within an
EC2 instance `nothing` is returned.
"""
instance_region() = chop(instance_availability_zone())

end


text:
using Clang
using Clang.LibClang.Clang_jll

const CHIPMUNK_INCLUDE = "Chipmunk2D/include/chipmunk"

const CHIPMUNK_HEADERS = [
    "Chipmunk2D/include/chipmunk/chipmunk.h",
    "Chipmunk2D/include/chipmunk/chipmunk_structs.h",
    "Chipmunk2D/include/chipmunk/chipmunk_types.h",
]

const EXCLUDED_CURSOR_NAMES = Set([
    "cpcalloc",
    "cpfree",
    "cprealloc",

    "cpfsqrt",
    "cpfsin",
    "cpfcos",
    "cpfacos",
    "cpfatan2",
    "cpfmod",
    "cpfexp",
    "cpfpow",
    "cpffloor",
    "cpfceil",
    "CPFLOAT_MIN",
])

wrap_context = init(
    headers=CHIPMUNK_HEADERS,
    output_file=joinpath(@__DIR__, "libchipmunk_api.jl"),
    common_file=joinpath(@__DIR__, "libchipmunk_common.jl"),
    clang_includes=vcat(CHIPMUNK_INCLUDE, CLANG_INCLUDE),
    clang_args=["-I", joinpath(CHIPMUNK_INCLUDE, "..")],

    cursor_wrapped=((cursor_name, cursor) -> begin
        !in(cursor_name, EXCLUDED_CURSOR_NAMES)
    end),

    header_wrapped=((root, current) -> begin
        # current == root
        startswith(current, "Chipmunk2D/include/chipmunk/")
    end),

    header_library=(header_name -> "libchipmunk_name"),
    clang_diagnostics=true,
)

run(wrap_context)


text:
##########################################################################################
# Basic imports and dimension definitions
##########################################################################################
using DimensionalData
using DimensionalData: @dim, hasdim, Dimension, IndependentDim
using DimensionalData: basetypeof
using Dates

Time = DimensionalData.Ti

AbDimArray = DimensionalData.AbstractDimArray

dimindex(A, i::Int) = i
dimindex(A, dim) = DimensionalData.dimnum(A, dim)

export At, Between, Near # Selectors from DimensionalArrays.jl
export hasdim, dimnum, dims
export Time, Lon, Lat, dims, Coord, Hei, Pre, Ti
export UnstructuredGrid, LonLatGrid, spacestructure
export DimensionalData # for accessing its functions

@dim Lon IndependentDim "Longitude"
@dim Lat IndependentDim "Latitude"
@dim Hei IndependentDim "Height"
@dim Pre IndependentDim "Pressure"

STANDARD_DIMS = (Lon, Lat, Time, Hei, Pre, Coord)

"""
    COMMONNAMES
A dictionary of common names of dimensions (as strings) to actual dimension types.
"""
const COMMONNAMES = Dict(
    "lat" => Lat,
    "latitude" => Lat,
    "y" => Lat,
    "yc" => Lat,
    "lon" => Lon,
    "long" => Lon,
    "x" => Lon,
    "xc" => Lon,
    "longitude" => Lon,
    "time" => Time,
    "t" => Time,
    "height" => Hei,
    "altitude" => Hei,
    "pressure" => Pre,
    "level" => Pre,
)

##########################################################################################
# Space types
##########################################################################################

# the following traits for the the way space is configured. currently the options are
# UnstructuredGrid, which is for equal area "grids" (or better points/coordinates)
# while the LonLatGrid is for standard Longitude x Latitude dimensions.
# Dispatch on `UnstructuredGrid` or other types is type-stable
# because it comes from the underlying dimension

abstract type SpaceType end

"""
Space coordinates are represented by two orthogonal dimensions `Lon, Lat`,
one being longitude and the other being latitude.
"""
struct LonLatGrid <: SpaceType end

"""
Space coordinates are represented by a single dimension `Coord`, whose
elements are coordinate locations, i.e. 2-element `SVector(longitude, latitude)`.
Each coordinate represents an **equal area polygon** corresponding to the point in space.
The actual limits of each polygon are not included in the dimension for performance reasons.

This dimension allows indexing according to the underlying `Lon, Lat` representation,
e.g. you can do
```julia
A # some `ClimArray` with unstructured grid type.
A[Coord(Lon(Between(0, 30)), Lat(Between(-30, 30)))]
```

To use functions such as [`zonalmean`](@ref) or [`hemispheric_means`](@ref) with this grid,
you must first sort the `ClimArray` so that the latitudes
of its coordinates are sorted in ascending order. I.e.
```julia
A # some `ClimArray` with unstructured grid type.
coords = dims(A, Coord).val
si = sortperm(coords, by = reverse)
A = A[Coord(si)]
```
**This is done automatically by [`ncread`](@ref).**

!!! warn
    `UnstructuredGrid` functionality is currently in an **experimental phase**!
    Notice that non-equal area unstructured grids are not supported yet.
"""
struct UnstructuredGrid <: SpaceType end

spacestructure(a::AbDimArray) = spacestructure(dims(a))
function spacestructure(dims)
    if hasdim(dims, Coord)
        UnstructuredGrid()
    elseif hasdim(dims, Lon) || hasdim(dims, Lat)
        LonLatGrid()
    else
        error("Array does not have any spatial dimensions: `Lon`, `Lat`, or `Coord`.")
    end
end

##########################################################################################
# ClimArray definition and DimensionalData.jl extensions
##########################################################################################
export ClimArray
struct ClimArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Me} <: AbstractDimensionalArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Symbol
    attrib::Me
end
ClimArray(A::DimensionalArray) = ClimArray(A.data, A.dims, A.refdims, A.name, A.metadata)
ClimArray(A::ClimArray, dims::Tuple = A.dims; name = A.name, attrib = A.attrib) =
ClimArray(A.data, DimensionalData.formatdims(A.data, dims), A.refdims, Symbol(name), attrib)

"""
    ClimArray(A::Array, dims::Tuple; name = "", attrib = nothing)

`ClimArray` is a structure that contains numerical array data bundled with dimensional
information, a name and an `attrib` field (typically a dictionary) that holds general
attributes.
You can think of `ClimArray` as a in-memory representation of a CFVariable.

At the moment, a `ClimArray` is using `DimensionalArray` from DimensionalData.jl, and
all basic handling of `ClimArray` is offered by `DimensionalData` (see below).

`ClimArray` is created by passing in standard array data `A` and a
tuple of dimensions `dims`. See [`ncread`](@ref) to automatically create a `ClimArray`
from a .nc file.

## Example
```julia
using ClimateBase, Dates
Time = ClimateBase.Ti # more intuitive name for time dimension
lats = -90:5:90
lons = 0:10:359
t = Date(2000, 3, 15):Month(1):Date(2020, 3, 15)
# dimensional information:
dimensions = (Lon(lons), Lat(lats), Time(t))
data = rand(36, 37, 241) # numeric data
A = ClimArray(data, dimensions)
```
"""
ClimArray(A::AbstractArray, dims::Tuple; refdims=(), name="", attrib=nothing) =
    ClimArray(A, DimensionalData.formatdims(A, dims), refdims, Symbol(name), attrib)
ClimArray(A::AbstractArray, dims::Tuple, name; refdims=(), attrib=nothing) =
    ClimArray(A, DimensionalData.formatdims(A, dims), refdims, Symbol(name), attrib)

Base.parent(A::ClimArray) = A.data
Base.@propagate_inbounds Base.setindex!(A::ClimArray, x, I::Vararg{DimensionalData.StandardIndices}) =
    setindex!(A.data, x, I...)

DimensionalData.metadata(A::ClimArray) = A.attrib
DimensionalData.rebuild(A::ClimArray, data::Any, dims::Tuple=dims(A), refdims=DimensionalData.refdims(A),
name="", attrib=nothing) = ClimArray(data, dims, refdims, Symbol(name), attrib)
DimensionalData.basetypeof(::ClimArray) = ClimArray

DimensionalData.rebuild(
    A::ClimArray;
    data=data(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A)
) = ClimArray(data, dims, refdims, name, metadata)

# The following basic methods allow indexing with tuples, (Time(5), Lon(3))
Base.getindex(A::ClimArray, i::Tuple) = A[i...]
Base.setindex!(A::ClimArray, x, i::Tuple) = setindex!(A, x, i...)
Base.view(A::ClimArray, i::Tuple) = view(A, i...)

# Convenience
Base.ones(A::AbDimArray) = basetypeof(A)(ones(eltype(A), size(A)), dims(A))
Base.zeros(A::AbDimArray) = basetypeof(A)(zeros(eltype(A), size(A)), dims(A))


text:
"""
    Arrow.jl

A pure Julia implementation of the [apache arrow](https://arrow.apache.org/) memory format specification.

This implementation supports the 1.0 version of the specification, including support for:
  * All primitive data types
  * All nested data types
  * Dictionary encodings and messages
  * Extension types
  * Streaming, file, record batch, and replacement and isdelta dictionary messages

It currently doesn't include support for:
  * Tensors or sparse tensors
  * Flight RPC
  * C data interface

Third-party data formats:
  * csv and parquet support via the existing CSV.jl and Parquet.jl packages
  * Other Tables.jl-compatible packages automatically supported (DataFrames.jl, JSONTables.jl, JuliaDB.jl, SQLite.jl, MySQL.jl, JDBC.jl, ODBC.jl, XLSX.jl, etc.)
  * No current Julia packages support ORC or Avro data formats

See docs for official Arrow.jl API with `Arrow.Table`, `Arrow.write`, and `Arrow.Stream`.
"""
module Arrow

using Mmap
import Dates
using DataAPI, Tables, SentinelArrays, PooledArrays, CodecLz4, CodecZstd

using Base: @propagate_inbounds
import Base: ==

const DEBUG_LEVEL = Ref(0)

function setdebug!(level::Int)
    DEBUG_LEVEL[] = level
    return
end

function withdebug(f, level)
    lvl = DEBUG_LEVEL[]
    try
        setdebug!(level)
        f()
    finally
        setdebug!(lvl)
    end
end

macro debug(level, msg)
    esc(quote
        if DEBUG_LEVEL[] >= $level
            println(string("DEBUG: ", $(QuoteNode(__source__.file)), ":", $(QuoteNode(__source__.line)), " ", $msg))
        end
    end)
end

const ALIGNMENT = 8
const FILE_FORMAT_MAGIC_BYTES = b"ARROW1"
const CONTINUATION_INDICATOR_BYTES = 0xffffffff

# vendored flatbuffers code for now
include("FlatBuffers/FlatBuffers.jl")
using .FlatBuffers

include("metadata/Flatbuf.jl")
using .Flatbuf; const Meta = Flatbuf

include("arrowtypes.jl")
using .ArrowTypes
include("utils.jl")
include("arraytypes.jl")
include("eltypes.jl")
include("table.jl")
include("write.jl")

end  # module Arrow

max_stars_repo_name: StefanKarpinski/arrow

max_issues_repo_name: StefanKarpinski/arrow

max_forks_repo_name: StefanKarpinski/arrow


text:
using Documenter, PackageAnalyzer

DocMeta.setdocmeta!(PackageAnalyzer, :DocTestSetup, :(using PackageAnalyzer); recursive=true)

makedocs(
    format = Documenter.HTML(
        prettyurls = true
    ),
    modules = [PackageAnalyzer],
    sitename = "PackageAnalyzer.jl",
    pages = ["Home" => "index.md",
             "API Reference" => "api.md",
             "Saving results" => "serialization.md"]
)

deploydocs(
    repo = "github.com/JuliaEcosystem/PackageAnalyzer.jl.git",
    push_preview=true,
    devbranch = "main",
)


text:
module ElasticWave
    using DSP, Requires
    const spmatveclib = abspath(joinpath(splitdir(Base.source_path())[1],"..","deps","builds","spmatvec"))
    include("dataStructure/dataStructure.jl")
    include("forwardAndAdjoint/forwardAndAdjoint.jl")
    include("spmatvec/spmatvec.jl")
    include("MOD/MOD.jl")
    @require PyPlot include("plotting/plotting.jl")
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
using QPDAS, Random, Test, LinearAlgebra
import QPDAS: deleterowcol!, addrowcol!, CholeskySpecialShifted, ldiv2!

@testset "Projections and solutions" begin
@testset "Projections onto nullspace, $ϵ, $lambdamin, $closesol" for (ϵ, lambdamin, closesol) in [(1e-8, false, false), (1e-5, false, false), (1e-6, false, false), (1e-7, false, false),
                                 (1e-8, true , false), (1e-5, true , false), (1e-6, true , false), (1e-7, true , false),
                                 (1e-8, false, true ), (1e-5, false, true ), (1e-6, false, true ), (1e-7, false, true ),
                                 (1e-8, true , true ), (1e-5, true , true ), (1e-6, true , true ), (1e-7, true , true )]

Random.seed!(3141592)

m, n = 60, 40

A = randn(m,n)
M = A*A'
b = randn(m)

# Almost sinugar
F2 = svd(M, full=true)
F2.S[n+1] += 1e-1
F2.S[n+2] += 1e-3
F2.S[n+3] += 1e-5
if lambdamin
    M = F2.U*Diagonal(F2.S)*F2.U'
end

F = CholeskySpecialShifted(Symmetric(M), ϵ)

b = randn(m)

#  # Almost solution
x0 = randn(m)
# d contains nothing in range(M)
d = 1e-5*(I-M*pinv(M))*randn(m)
if closesol
    b = M*x0
    b .+= d
end

@testset "no sol" begin
    x = copy(b)
    #data1 = Data()
    xsol, projection = ldiv2!(F, x)#, data=data1)
    @test projection == true
    @test norm((I-M*pinv(M))*b - xsol) < 2/ϵ*1e-11
    @test norm(M*xsol) < 1e-8
    #println("xsol $closesol $lambdamin: $(dot(b, xsol)/norm(b)/norm(xsol))")
    #data2 = Data()
    x2sol, projection2 = ldiv2!(F, 0*x, x0=b)#, data=data2)
    @test norm((I-M*pinv(M))*b-x2sol) < 5/ϵ*1e-11 #This seems to work better
    #println("x2sol $closesol $lambdamin: $(dot(b, x2sol)/norm(b)/norm(x2sol))")
    @test norm(M*x2sol) < 1e-8
end

@testset "test sol zero residual?" begin
    b = M*x0
    x = copy(b)
    #data2 = Data()
    xsol, projection = ldiv2!(F, x)
    @test projection == false
    @test norm(M*xsol - b) < 1e-10
end

@testset "test sol ortogonal" begin
    # We should get x0 back
    b = M*x0
    x02 = x0 + pinv(M)*M*randn(m) # Add something orthogonal to Mx=0
    xsol, projection = ldiv2!(F, copy(b), x0=x02)
    @test projection == false
    @test norm(M*xsol - b) < 1e-12
    @test norm(xsol - x0) < 1e-4
end

@testset "test sol projection from zero" begin
    # We should get x0 back
    x0 = pinv(M)*M*randn(m) # Something orthogonal to Mx=0 (should be solution from 0)
    b = M*x0
    xsol, projection = ldiv2!(F, copy(b), x0=0)
    @test projection == false
    @test norm(M*xsol - b) < 1e-12
    @test norm(xsol - x0) < 1e-4
end

end

end

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
# tut2_widgets.jl copied from Dash tutorial

#=
The DashCoreComponents package includes a set of higher level components like dropdowns, graphs, markdown blocks, and more. Like all Dash components, they are described entirely declaratively. Every option that is configurable is available as a keyword argument to the component.
=#

using Dash, DashHtmlComponents, DashCoreComponents

# app = dash()
app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

dropdown_options = [
	Dict("label" => "New York City", "value" => "NYC"),
	Dict("label" => "Montreal", "value" => "MTL"),
	Dict("label" => "San Francisco", "value" => "SF"),
]
app.layout = html_div(style = Dict("columnCount" => 2)) do
	html_label("Dropdown"),
	dcc_dropdown(options = dropdown_options, value = "MTL"),
	html_label("Multi-Select Dropdown"),
	dcc_dropdown(
		options = dropdown_options,
		value = ["MTL", "SF"],
		multi = true,
	),
	html_label("Radio Items"),
	dcc_radioitems(options = dropdown_options, value = "MTL"),
	html_label("Checkboxes"),
	dcc_checklist(options = dropdown_options, value = ["MTL", "SF"]),
	html_label("Text Input"),
	dcc_input(value = "MTL", type = "text"),
	html_label("Slider"),
	dcc_slider(
		min = 0,
		max = 9,
		marks = Dict([i => (i == 1 ? "Label $(i)" : "$(i)") for i = 1:6]),
		value = 5,
	)
end

run_server(app, "0.0.0.0", 8000, debug = true)

max_stars_repo_name: ProfessorUdGuru/treebard

max_issues_repo_name: ProfessorUdGuru/treebard

max_forks_repo_name: ProfessorUdGuru/treebard