text:
################################################################################
# api.jl
#
# The core API of the Air library. Abstract types are defined here as well as
# a number of functions that are common across many of the persistent
# collections. Functions that are defined for core Julia types (like Tuple and
# Array) are defined here as well.
#
# @author Noah C. Benson
#
# MIT License
# Copyright (c) 2020-2021 Noah C. Benson

import Base.setindex
import Base.delete!

################################################################################
# Abstract Types
"""
    AbstractPDict{K,V}

`AbstractPDict` is a subtype of `AbstractDict` that is extended only by
persistent dictionary types such as `PDict` and `LazyDict`.

See also [`PDict`](@ref), [`PWDict`](@ref), [`LazyDict`](@ref),
[`PIdDict`](@ref), [`LazyIdDict`](@ref), [`PWIdDict`](@ref).
"""
abstract type AbstractPDict{K,V} <: AbstractDict{K,V} end
"""
    AbstractPSet{T}

`AbstractPSet` is an abstract type extended by persistent set types such as
 `PSet` and `PWSet.

See also: [`PSet`](@ref), [`PIdSet`](@ref), [`PWSet`](@ref), [`PWIdSet`](@ref).
"""
abstract type AbstractPSet{T} <: AbstractSet{T} end
"""
    AbstractPArray{T,N}

The Abstract persistent array type represents any `Array`-like type that is
persistent. The Air library provides the reified type `PArray{T,N}`.

See also: [`PArray`](@ref), [`PVector`](@ref).
"""
abstract type AbstractPArray{T,N} <: AbstractArray{T,N} end
export AbstractPDict, AbstractPSet, AbstractPArray
"""
    ReentrantRef{T}

A `ReentrantRef` is a type of ref object that is thread-safe. There are a few
strategies for this, each of which are encoded in a different object type that
inheits from `ReentrantRef`. These types are as follows.

`Var`: `Var` objects act like `Ref` objects except that changes to them are
exclusively task-local and must be performed in specific scoped expressions.
However, the state of all `Var` objects for the current task can also be
saved and restored at a later point.

`Actor`: `Actor` objects obey the actor pattern; you can call
`send(fn, actor, args...)` where `fn` is a function that is, in another thread
at some point, called as `fn(actor[], args...)`. The new value of the actor
after running `fn` is the return value of the call.

`Volatile`: `Volatile` objects are `Ref`s that can be changed by any thread but
that must be changed only within a synchronized transaction that ensures that
all reads and writes of volatiles, as well as reads from and sends to actors,
are atomic: either they all happen successfully, or none of them do.
"""
abstract type ReentrantRef{T} <: Ref{T} end
export ReentrantRef
"""
    TransactionalRef{T}

A `TransactionalRef` is a reentrant reference that additional participates in
transactions. Transactional refs include `Actor`s and `Volatile`s.
"""
abstract type TransactionalRef{T} <: ReentrantRef{T} end
export TransactionalRef

# #equalfn and #hashfn #########################################################
"""
    equalfn(x)

If `x` is an object (such as a `PSet` or `Dict`) that has an opinion about
equality, `equalfn(x)` returns the function that it uses.

See also: [`hashfn`](@ref)

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> equalfn(Dict) == isequal
true

julia> equalfn(IdDict) == (===)
true

julia> equalfn(PSet) == isequal
true

julia> equalfn(PArray)
ERROR: ArgumentError: no equalfn for type PArray
```
"""
equalfn(x::T) where {T} = equalfn(T)
equalfn(x::Type{T}) where {T} = throw(ArgumentError("no equalfn for type $T"))
equalfn(x::Type{Dict}) = isequal
equalfn(x::Type{Set}) = isequal
equalfn(x::Type{IdDict}) = (===)
equalfn(x::Type{Base.IdSet}) = (===)
equalfn(x::Type{Dict{K,V}}) where {K,V} = isequal
equalfn(x::Type{Set{T}}) where {T} = isequal
equalfn(x::Type{IdDict{K,V}}) where {K,V} = (===)
equalfn(x::Type{Base.IdSet{T}}) where {T} = (===)
"""
    hashfn(x)

If `x` is an object (such as a `PSet` or `Dict`) that has an opinion about how
it hashes objects, `hashfn(x)` returns the function that it uses. It is
sufficient in almost all circumstances to define `equalfn(T)`; the `hashfn`
should always match the `equalfn` regardless.

See also: [`equalfn`](@ref)

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> hashfn(Dict) == hash
true

julia> hashfn(IdDict) == objectid
true

julia> hashfn(PSet) == hash
true

julia> hashfn(PArray)
ERROR: ArgumentError: no equalfn for type PArray
```
"""
hashfn(x::T) where {T} = hashfn(equalfn(T))
hashfn(::Type{T}) where {T} = hashfn(equalfn(T))
# equalfn and hashfn can also be used on their respective functions
equalfn(::typeof(objectid))  = (===)
equalfn(::typeof(hash))      = isequal
hashfn(::typeof(===))     = objectid
hashfn(::typeof(isequal)) = hash
export equalfn, hashfn

# #setindex ====================================================================
"""
    setindex(coll, val, index)

Yields a copy of the given collection `coll` with the value `val` set at the
given index. This is a persistent version of `setindex!` and works for most
collection types including `Array`s, `Dict`s, `IdDict`s, and Air's persistent
versions of these.

Note that setindex() *always* returns a copy of the argument `coll` or fails.
For Air's persistent collections these operations are efficient, but for the
mutable counterparts, they are typically `O(n)`.

See also: `Base.setindex!`, [`push`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> setindex((:a,:b,:e,:d), :c, 3)
(:a, :b, :c, :d)
```

```jldoctest; filter=r"4-element (Array{Symbol,1}|Vector{Symbol}):"
julia> u = [:a,:b,:e,:d]; v = setindex(u, :c, 3)
4-element Array{Symbol,1}:
 :a
 :b
 :c
 :d

julia> u == v
false

julia> u[3]
:e
```

```jldoctest; filter=[r"Dict{Any, ?Any} with 1 entry:", r"Dict{Any, ?Any}()"]
julia> d1 = Dict(); setindex(d1, 10, :a)
Dict{Any,Any} with 1 entry:
  :a => 10

julia> d1
Dict{Any,Any}()
```
"""
setindex(u::AbstractArray{T,N}, x::S, I...) where {T,N,S} =
    Base.setindex!(copy(u), x, I...)
setindex(d::AbstractDict{K,V}, v::U, k::J) where {K,V,U,J} =
    Base.setindex!(copy(d), v, k)
export setindex

# #push ========================================================================
"""
    push(coll, val)
    push(coll1, val1, val2...)

Yields a copy of the given collection `coll` with the given value `val`
appended. This function is essentially a persistent equivalent of the `push!`
function that never modifies the object `coll`. For persistent collections in
the Air library, this operation is efficient, but for most mutable objects, it
is `O(n)`.

See also: `Base.push!`, [`pop`](@ref), `Base.setindex!`.

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"5-element (Array{Symbol,1}|Vector{Symbol}):"
julia> push((:a,:b,:c,:d), :e)
(:a, :b, :c, :d, :e)
```

```jldoctest; filter=r"5-element (Array{Symbol,1}|Vector{Symbol}):"
julia> u = [:a,:b,:c,:d]; push(u, :e)
5-element Array{Symbol,1}:
 :a
 :b
 :c
 :d
 :e
```

```jldoctest
julia> s = Set([:a, :b, :c]); push(s, :d)
Set{Symbol} with 4 elements:
  :a
  :b
  :d
  :c
```

```jldoctest; filter=r"Dict{Symbol, ?Int64} with 4 entries:"
julia> d = Dict(:a => 1, :b => 2); push(d, :c => 3, :d => 4)
Dict{Symbol,Int64} with 4 entries:
  :a => 1
  :b => 2
  :d => 4
  :c => 3
```
"""
function push end
push(A, a, b, c...) = reduce(push, c, init=push(push(A, a), b))
push(u::NTuple{N,T}, val::S) where {T,N,S} = begin
    return NTuple{N+1,T}(T[u..., val])
end
push(u::AbstractVector{T}, val::S) where {T,S} = T[u..., val]
push(d::AbstractDict{K,V}, kv::Pair{J,U}) where {K,V,J,U} = begin
    return setindex(d, kv[2], kv[1])
end
push(s::AbstractSet{T}, u::S) where {T,S} = push!(copy(s), u)
export push

# #pop =========================================================================
"""
    pop(coll)

Yields a tuple `(last, most)` where `last` is the last element of the given
collection and most is a duplicate tuple of all but the last element of tup.
This is basically a persistent equivalent to the `pop!` function that never
modifies the given `coll`. For persistent collections in Air, these operations
are efficient, but for most mutable types, it is `O(n)`.

    pop(coll, key[, default])

Similar to `pop!`, pops the specific `key` from the given collection `coll` and
yields `(val, rest)` where `val` is the value associated with `key` in `coll`
and `rest` is the remainder of the collection without `key`. If the `key` is not
in `coll`, then `default` is yielded or an error is thrown.

See also: `Base.pop!`, [`push`](@ref), [`popat`](@ref),
[`delete`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> pop((:a,:b,:c,:d))
(:d, (:a, :b, :c))

julia> u = [:a,:b,:c,:d]; pop(u)
(:d, [:a, :b, :c])

julia> length(u)
4

julia> s = Set([:a, :b, :c]); pop(s, :c)
(:c, Set([:a, :b]))

julia> :c in s
true
```

```jldoctest; filter=r"\\(3, Dict\\(:a => 1, ?:b => 2\\)\\)"
julia> d = Dict(:a => 1, :b => 2, :c => 3); pop(d, :c)
(3, Dict(:a => 1,:b => 2))

julia> d[:c]
3
```
"""
function pop end
pop(u::NTuple{N,T}) where {N,T} = (u[end], u[1:end-1])
pop(u::Tuple{}) = throw(ArgumentError("n-tuple must be non-empty"))
pop(u::AbstractArray{T,1}) where {T} = (u[end], copy(u[1:end-1]))
pop(d::AbstractDict{K,V}) where {K,V} = begin
    rem = copy(d)
    v = pop!(rem)
    return (v, rem)
end
pop(d::AbstractDict{K,V}, k::J) where {K,V,J} = begin
    rem = copy(d)
    v = pop!(rem, k)
    return (v, rem)
end
pop(d::AbstractPDict{K,V}) where {K,V} = begin
    (length(d) == 0) && throw(
        ArgumentError("cannot pop from an empty collection"))
    kv = iterate(d)[1]
    return (kv, delete(d, kv[1]))
end
pop(d::AbstractPDict{K,V}, k::J) where {K,V,J} = begin
    kv = getpair(d, k)
    (kv === missing) && throw(KeyError(k))
    return (kv, delete(d, k))
end
pop(d::AbstractDict{K,V}, k::J, dv) where {K,V,J} = begin
    rem = copy(d)
    v = pop!(rem, k, dv)
    return (v, rem)
end
pop(s::AbstractSet{T}) where {T,S} = begin
    rem = copy(s)
    u = pop!(rem)
    return (u, rem)
end
pop(s::AbstractSet{T}, u::S) where {T,S} = begin
    rem = copy(s)
    u = pop!(rem, u)
    return (u, rem)
end
pop(s::AbstractPSet{T}) where {T,S} = begin
    (length(d) == 0) && throw(
        ArgumentError("cannot pop from an empty collection"))
    x = iterate(d)[1]
    return (x, delete(d, x))
end
pop(s::AbstractPSet{T}, u::S) where {T,S} = begin
    rem = copy(s)
    u = pop!(rem, u)
    return (u, rem)
end
export pop

# #pushfirst ===================================================================
"""
    pushfirst(coll, val)
    pushfirst(coll1, val1, val2...)

Yields a copy of the given collection `coll` with the given value `val`
prepended. This function is essentially a persistent equivalent of the
`pushfirst!` function that never modifies the object `coll`. For persistent
collections in the Air library, this operation is efficient, but for most
mutable objects, it is `O(n)`.

See also: `Base.pushfirst!`, [`popfirst`](@ref), `Base.setindex!`,
[`push`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest; filter=r"5-element (Array{Symbol,1}|Vector{Symbol}):"
julia> pushfirst((:a,:b,:c,:d), :e)
(:e, :a, :b, :c, :d)

julia> u = [:a,:b,:c,:d]; pushfirst(u, :e)
5-element Array{Symbol,1}:
 :e
 :a
 :b
 :c
 :d

julia> length(u)
4
```
"""
function pushfirst end
pushfirst(u::NTuple{N,T}, val::S) where {T,N,S} = begin
    return NTuple{N+1,T}(T[val, u...])
end
pushfirst(u::AbstractVector{T}, val::S) where {T,S} = T[val, u...]
pushfirst(A, a, b...) = push(push(A, a), b...)
export pushfirst

# #popfirst ====================================================================
"""
    popfirst(coll)

Yields a tuple `(first, rest)` where `first` is the first element of the given
collection and rest is a duplicate tuple of all but the first element of `tup`.
This is basically a persistent equivalent to the `popfirst!` function that never
modifies the given `coll`. For persistent collections in Air, these operations
are efficient, but for most mutable types, it is `O(n)`.

See also: `Base.popfirst!`, [`pushfirst`](@ref), [`popat`](@ref),
[`delete`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> popfirst((:a,:b,:c,:d))
(:a, (:b, :c, :d))

julia> u = [:a,:b,:c,:d]; popfirst(u)
(:a, [:b, :c, :d])

julia> length(u)
4
```
"""
function popfirst end
popfirst(u::NTuple{N,T}) where {N,T} = (u[1], u[2:end])
popfirst(u::Tuple{}) = throw(ArgumentError("n-tuple must be non-empty"))
popfirst(u::AbstractArray{T,1}) where {T} = (u[1], copy(u[2:end]))
export popfirst

# #popat =======================================================================
"""
    popat(coll, k)

Yields a tuple `(el, rest)` where `el` is the `k`th element of the given
collection `coll` and `rest` is a duplicate of all but the `k`th element of
`coll`. This is basically a persistent equivalent to the `popat!` function that
never modifies the given `coll`. For persistent collections in Air, these
operations are efficient, but for most mutable types, it is `O(n)`.

See also: `Base.popat!`, [`pop`](@ref), [`popfirst`](@ref),
[`delete`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> popat((:a,:b,:c,:d), 2)
(:b, (:a, :c, :d))

julia> u = [:a,:b,:c,:d]; popat(u, 3)
(:c, [:a, :b, :d])

julia> length(u)
4
```
"""
function popat end
popat(u::NTuple{N,T}, k::K) where {N,T,K<:Integer} =
    (u[k], (u[1:k-1]..., u[k+1:end]...))
popat(u::Tuple{}) = throw(ArgumentError("n-tuple must be non-empty"))
popat(u::AbstractArray{T,1}, k::K) where {T,K<:Integer} = begin
    c = copy(u)
    el = popat!(c, k)
    return (el, c)
end
export popat

# #insert ======================================================================
"""
    insert(coll, idx, val)

Yields a copy of the given collection `coll` with the given value `cal` inserted
at the given index. Roughly equivalent to `insert!(copy(arr), idx, va)`: the
`insert` function never modifies its arguments and always yields a copy. For the
persistent collections defined in Air, this operation is efficient, but for most
mutable objects, this is `O(n)`.

See also: `Base.insert!`, [`push`](@ref), [`delete`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> insert((:a,:b,:c,:d), 2, :x)
(:a, :x, :b, :c, :d)
```

```jldoctest; filter=r"5-element (Array{Symbol,1}|Vector{Symbol}):"
julia> u = [:a,:b,:c,:d]; insert(u, 1, :x)
5-element Array{Symbol,1}:
 :x
 :a
 :b
 :c
 :d

julia> length(u)
4
```

```jldoctest; filter=r"5-element Bit.+:"
julia> b = BitArray([0,0,0,0]); insert(b, 2, 1)
5-element BitArray{1}:
 0
 1
 0
 0
 0

julia> b[2]
false
```
"""
function insert end
insert(a::AbstractVector{T}, k::K, val::S) where {T,K<:Integer,S} =
    insert!(copy(a), k, val)
insert(a::Vector{T}, k::K, val::S) where {T,K<:Integer,S} = begin
    n = length(a) + 1
    (k > n) && throw(
        ArgumentError("insert: $k is out of range for Vector of length $(n-1)"))
    out = Vector{T}(undef, n)
    (k > 1) && copyto!(out, 1, a, 1, k - 1)
    @inbounds out[k] = val
    (k < n) && copyto!(out, k + 1, a, k, n - k)
    return out
end
insert(a::Tuple, idx::II, val::T) where {II<:Integer,T} =
    (a[1:idx]..., val, a[idx:end]...)
insert(a::Tuple{}, idx::II, val::T) where {II<:Integer,T} = begin
    if idx == 1
        return (val,)
    else
        throw(ArgumentError("invalid index $idx for tuple of size 0"))
    end
end
# For tuples, we want to generate versions of this for NTuples up to size 64;
# beyond that we can use a generic function.
macro _tuple_insert_gencode(N::Int)
    # We'll want to refer to the tuple elements:
    els = [:(tup[$k]) for k in 1:N]
    # Build up the if-elseif-else expression, starting with the else:
    ifexpr = :(throw(ArgumentError(
        "insert: $k if out of range for a Tuple of length $(length(tup))")))
    ifexpr = Expr(:elseif, :(k == $(N+1)), :(push(tup, el)), ifexpr)
    for k in N:-1:1
        ifexpr = Expr(k == 1 ? :if : :elseif,
                      :(k == $k),
                      :(($(els[1:k-1]...), el, $(els[k:end]...))),
                      ifexpr)
    end
    return quote
        insert(tup::NTuple{$N,T}, k::K, el::S) where {T,K<:Integer,S} = $ifexpr
    end |> esc
end
# Now generate functions for up to 64:
(@_tuple_insert_gencode  1)
(@_tuple_insert_gencode  2)
(@_tuple_insert_gencode  3)
(@_tuple_insert_gencode  4)
(@_tuple_insert_gencode  5)
(@_tuple_insert_gencode  6)
(@_tuple_insert_gencode  7)
(@_tuple_insert_gencode  8)
(@_tuple_insert_gencode  9)
(@_tuple_insert_gencode 10)
(@_tuple_insert_gencode 11)
(@_tuple_insert_gencode 12)
(@_tuple_insert_gencode 13)
(@_tuple_insert_gencode 14)
(@_tuple_insert_gencode 15)
(@_tuple_insert_gencode 16)
(@_tuple_insert_gencode 17)
(@_tuple_insert_gencode 18)
(@_tuple_insert_gencode 19)
(@_tuple_insert_gencode 20)
(@_tuple_insert_gencode 21)
(@_tuple_insert_gencode 22)
(@_tuple_insert_gencode 23)
(@_tuple_insert_gencode 24)
(@_tuple_insert_gencode 25)
(@_tuple_insert_gencode 26)
(@_tuple_insert_gencode 27)
(@_tuple_insert_gencode 28)
(@_tuple_insert_gencode 29)
(@_tuple_insert_gencode 30)
(@_tuple_insert_gencode 31)
(@_tuple_insert_gencode 32)
(@_tuple_insert_gencode 33)
(@_tuple_insert_gencode 34)
(@_tuple_insert_gencode 35)
(@_tuple_insert_gencode 36)
(@_tuple_insert_gencode 37)
(@_tuple_insert_gencode 38)
(@_tuple_insert_gencode 39)
(@_tuple_insert_gencode 40)
(@_tuple_insert_gencode 41)
(@_tuple_insert_gencode 42)
(@_tuple_insert_gencode 43)
(@_tuple_insert_gencode 44)
(@_tuple_insert_gencode 45)
(@_tuple_insert_gencode 46)
(@_tuple_insert_gencode 47)
(@_tuple_insert_gencode 48)
(@_tuple_insert_gencode 49)
(@_tuple_insert_gencode 50)
(@_tuple_insert_gencode 51)
(@_tuple_insert_gencode 52)
(@_tuple_insert_gencode 53)
(@_tuple_insert_gencode 54)
(@_tuple_insert_gencode 55)
(@_tuple_insert_gencode 56)
(@_tuple_insert_gencode 57)
(@_tuple_insert_gencode 58)
(@_tuple_insert_gencode 59)
(@_tuple_insert_gencode 60)
(@_tuple_insert_gencode 61)
(@_tuple_insert_gencode 62)
(@_tuple_insert_gencode 63)
(@_tuple_insert_gencode 64)
export insert

# #delete ======================================================================
"""
    delete(coll, k)

Yields a copy of the collection `coll` with the item associated with the given
index or key `k` deleted. This function works on dictionaries, vectors, tuples,
and sets, and it always yields a copy of the collection without modifying its
aguments. For the persistent collections defined in Air, this is efficient, but
for most mutable objects, this is `O(n)`.

See also: `Base.delete!`, [`insert`](@ref), [`pop`](@ref).

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> delete((:a,:b,:c,:d), 2)
(:a, :c, :d)
```

```jldoctest; filter=r"3-element (Array{Symbol,1}|Vector{Symbol}):"
julia> u = [:a,:b,:c,:d]; delete(u, 1)
3-element Array{Symbol,1}:
 :b
 :c
 :d
```

```jldoctest
julia> s1 = Set([:a,:b,:c,:d]); s2 = delete(s1, :d); (:d in s1, :d in s2)
(true, false)
```

```jldoctest; filter=r"4-element Bit.+:"
julia> b = BitArray([0,0,0,1,0]); delete(b, 4)
4-element BitArray{1}:
 0
 0
 0
 0

julia> b[2]
false
```
"""
delete(d::AbstractDict{K,V}, k::J) where {K,V,J} = delete!(copy(d), k)
delete(d::AbstractSet{V}, k::U) where {V,U} = delete!(copy(d), k)
delete(u::AbstractVector{T}, k::K) where {T,K<:Integer} = begin
    n = length(u)
    V = typeof(u)
    (n == 0 || k < 0 || k > n) && throw(
        ArgumentError("delete: $k is out of range for Vector of length $n"))
    out = V(undef, n - 1)
    (k > 1) && copyto!(out, 1, u, 1, k - 1)
    (k < n) && copyto!(out, k, u, k + 1, n - k)
    return out
end
delete(a::Tuple, idx::K) where {K<:Integer,T} = (a[1:idx-1]..., a[idx+1:end]...)
delete(a::Tuple{}, idx::II) where {II<:Integer} =
    throw(ArgumentError("delete: $idx is out of range for a Tuple{}"))
# For tuples, we want to generate versions of this for NTuples up to size 64;
# beyond that we can use a generic function.
macro _tuple_delete_gencode(N::Int)
    # We'll want to refer to the tuple elements:
    els = [:(tup[$k]) for k in 1:N]
    # Build up the if-elseif-else expression, starting with the else:
    ifexpr = :(throw(ArgumentError(
        "delete: $k if out of range for a Tuple of length $(length(tup))")))
    for k in N:-1:1
        ifexpr = Expr(k == 1 ? :if : :elseif,
                      :(k == $k),
                      :(($(els[1:k-1]...), $(els[k+1:end]...))),
                      ifexpr)
    end
    return quote
        delete(tup::NTuple{$N,T}, k::K) where {T,K<:Integer} = $ifexpr
    end |> esc
end
# Now generate functions for up to 64:
(@_tuple_delete_gencode  1)
(@_tuple_delete_gencode  2)
(@_tuple_delete_gencode  3)
(@_tuple_delete_gencode  4)
(@_tuple_delete_gencode  5)
(@_tuple_delete_gencode  6)
(@_tuple_delete_gencode  7)
(@_tuple_delete_gencode  8)
(@_tuple_delete_gencode  9)
(@_tuple_delete_gencode 10)
(@_tuple_delete_gencode 11)
(@_tuple_delete_gencode 12)
(@_tuple_delete_gencode 13)
(@_tuple_delete_gencode 14)
(@_tuple_delete_gencode 15)
(@_tuple_delete_gencode 16)
(@_tuple_delete_gencode 17)
(@_tuple_delete_gencode 18)
(@_tuple_delete_gencode 19)
(@_tuple_delete_gencode 20)
(@_tuple_delete_gencode 21)
(@_tuple_delete_gencode 22)
(@_tuple_delete_gencode 23)
(@_tuple_delete_gencode 24)
(@_tuple_delete_gencode 25)
(@_tuple_delete_gencode 26)
(@_tuple_delete_gencode 27)
(@_tuple_delete_gencode 28)
(@_tuple_delete_gencode 29)
(@_tuple_delete_gencode 30)
(@_tuple_delete_gencode 31)
(@_tuple_delete_gencode 32)
(@_tuple_delete_gencode 33)
(@_tuple_delete_gencode 34)
(@_tuple_delete_gencode 35)
(@_tuple_delete_gencode 36)
(@_tuple_delete_gencode 37)
(@_tuple_delete_gencode 38)
(@_tuple_delete_gencode 39)
(@_tuple_delete_gencode 40)
(@_tuple_delete_gencode 41)
(@_tuple_delete_gencode 42)
(@_tuple_delete_gencode 43)
(@_tuple_delete_gencode 44)
(@_tuple_delete_gencode 45)
(@_tuple_delete_gencode 46)
(@_tuple_delete_gencode 47)
(@_tuple_delete_gencode 48)
(@_tuple_delete_gencode 49)
(@_tuple_delete_gencode 50)
(@_tuple_delete_gencode 51)
(@_tuple_delete_gencode 52)
(@_tuple_delete_gencode 53)
(@_tuple_delete_gencode 54)
(@_tuple_delete_gencode 55)
(@_tuple_delete_gencode 56)
(@_tuple_delete_gencode 57)
(@_tuple_delete_gencode 58)
(@_tuple_delete_gencode 59)
(@_tuple_delete_gencode 60)
(@_tuple_delete_gencode 61)
(@_tuple_delete_gencode 62)
(@_tuple_delete_gencode 63)
(@_tuple_delete_gencode 64)
export delete

# #getpair =====================================================================
"""
    getpair(d, k)

If the key `k` is found in the dictionary `d`, yields the pair (`k => d[k]`);
otherwise yields `missing`.

See also: `Base.get`, `Base.getindex`.

# Examples

```@meta
DocTestSetup = quote
    using Air
end
```

```jldoctest
julia> d = Dict(:a => 1, :b => 2, :c => 3); getpair(d, :a)
:a => 1

julia> getpair(d, :x)
missing
```
"""
function getpair end
let _private_symbol = gensym("private_getpair_c61d3fee7deae4d3d2237de9044247b2")
    global getpair(d::AbstractDict, k) = begin
        v = get(d, k, _private_symbol)
        return (v === _private_symbol) ? missing : (k => v)
    end
end
export getpair

# #isequal #####################################################################
function _equalfn_isequal(
    l::DL, r::DR
) where {
    KL,VL, DL <: AbstractDict{KL,VL},
    KR,VR, DR <: AbstractDict{KR,VR}
}
    # Identical dictionaries are always equal.
    (l === r) && return true
    # Equality functions must match.
    (equalfn(DL) === equalfn(DR)) || return false
    # Dictionaries must be the same length to be equal.
    (length(l) == length(r)) || return false
    # All pairs from one must be in the other.
    for pair in l
        in(pair, r, isequal) || return false
    end
    return true
end
Base.isequal(l::AbstractDict, r::AbstractPDict) = _equalfn_isequal(l, r)
Base.isequal(l::AbstractPDict, r::AbstractDict) = _equalfn_isequal(l, r)
Base.isequal(l::AbstractPDict, r::AbstractPDict) = _equalfn_isequal(l, r)
function _equalfn_isequal(
    l::SL, r::SR
) where {
    TL, SL <: AbstractSet{TL},
    TR, SR <: AbstractSet{TR}
}
    # Sets must have the same equality base to be equal.
    (equalfn(SL) === equalfn(SR)) || return false
    # Identical sets are always equal.
    (l === r) && return true
    # Sets must be the same length to be equal.
    (length(l) == length(r)) || return false
    # And all elements must be found in both.
    for el in l
        in(el, r) || return false
    end
    return true
end
Base.isequal(l::AbstractSet, r::AbstractPSet) = _equalfn_isequal(l, r)
Base.isequal(l::AbstractPSet, r::AbstractSet) = _equalfn_isequal(l, r)
Base.isequal(l::AbstractPSet, r::AbstractPSet) = _equalfn_isequal(l, r)
function _equalfn_eq(
    l::DL, r::DR
) where {
    KL,VL, DL <: AbstractDict{KL,VL},
    KR,VR, DR <: AbstractDict{KR,VR}
}
    # Dictionaries must have the same equality base to be equal.
    (equalfn(DL) === equalfn(DR)) || return false
    # Identical dictionaries are always equal.
    (l === r) && return true
    # Dictionaries must be the same length to be equal.
    (length(l) == length(r)) || return false
    # And all pairs must be found in both.
    anymissing = false
    for pair in l
        isin = in(pair, r)
        if ismissing(isin)
            anymissing = true
        else
            isin || return false
        end
    end
    return anymissing ? missing : true
end
Base.:(==)(l::AbstractDict, r::AbstractPDict) = _equalfn_eq(l, r)
Base.:(==)(l::AbstractPDict, r::AbstractDict) = _equalfn_eq(l, r)
Base.:(==)(l::AbstractPDict, r::AbstractPDict) = _equalfn_eq(l, r)
function _equalfn_eq(
    l::SL, r::SR
) where {
    TL, SL <: AbstractSet{TL},
    TR, SR <: AbstractSet{TR}
}
    # Sets must have the same equality base to be equal.
    (equalfn(SL) === equalfn(SR)) || return false
    # Identical sets are always equal.
    (l === r) && return true
    # Sets must be the same length to be equal.
    (length(l) == length(r)) || return false
    # And all elements must be found in both.
    for el in l
        in(el, r) || return false
    end
    return true
end
Base.:(==)(l::AbstractSet, r::AbstractPSet) = _equalfn_eq(l, r)
Base.:(==)(l::AbstractPSet, r::AbstractSet) = _equalfn_eq(l, r)
Base.:(==)(l::AbstractPSet, r::AbstractPSet) = _equalfn_eq(l, r)


text:
using HTTP, MbedTLS

ssl = MbedTLS.SSLConfig()
MbedTLS.config_defaults!(ssl)
MbedTLS.rng!(ssl, MbedTLS.seed!(MbedTLS.CtrDrbg(), MbedTLS.Entropy()))
MbedTLS.ca_chain!(ssl, MbedTLS.crt_parse_file("zyte-smartproxy-ca.crt"))
MbedTLS.authmode!(ssl, MbedTLS.MBEDTLS_SSL_VERIFY_REQUIRED)

proxy = "http://" * ENV["KEY"] * ":@" * ENV["PROXY"]

r = HTTP.get(ENV["URL"]; sslconfig=ssl, proxy=proxy)

print(String(r.body))

ext: jl

max_stars_repo_name: zytedata/zyte-smartproxy-clients

max_issues_repo_name: scrapinghub/crawlera-clients

max_forks_repo_name: zytedata/zyte-smartproxy-clients


text:
# Abstract path function
abstract type PathFunction{type <: FunctionType} end

# PathFunction contructors
function PathFunction(type::FunctionType, func!::Function, nFuncs::Int, nStates::Int, nControls::Int, nStatic::Int)
    return ADPathFunction(type, func!, nFuncs, nStates, nControls, nStatic)
end
function PathFunction(type::FunctionType, func!::Function, stateJac!::Union{Function,Nothing},
    controlJac!::Union{Function,Nothing}, staticJac!::Union{Function,Nothing}, timeJac!::Union{Function,Nothing},
    nFuncs::Int, nStates::Int, nControls::Int, nStatic::Int,
    stateSP::AbstractVecOrMat, controlSP::AbstractVecOrMat, staticSP::AbstractVecOrMat, timeSP::AbstractVecOrMat)
    return AnalyticPathFunction(type, func!, stateJac!, controlJac!, staticJac!, timeJac!,
        nFuncs, nStates, nControls, nStatic, stateSP, controlSP, staticSP, timeSP)
end

# Get the type of path function
function GetFunctionType(pf::PathFunction{type}) where {type}
    return type
end

# Evaluate functions method
function EvaluateFunction(pf::PathFunction,
                          out::AbstractVector,
                          state::AbstractVector,
                          control::AbstractVector,
                          static::AbstractVector,
                          time::AbstractFloat)
    pf.func!(out,state,control,static,time)
    return nothing
end

# Get number of parameters methods
function GetNumberOfFunctions(pf::PathFunction)
    return pf.nFuncs
end
function GetNumberOfStates(pf::PathFunction)
    return pf.nStates
end
function GetNumberOfControls(pf::PathFunction)
    return pf.nControls
end
function GetNumberOfStatics(pf::PathFunction)
    return pf.nStatic
end

# Function name methods
function SetFunctionNames!(pf::PathFunction, names::Vector{String})
    for i in 1:length(names)
        push!(pf.funcNames, name[i])
    end
    return nothing
end
function GetFunctionNames(pf::PathFunction)
    return pf.funcNames
end

# Get Jacobian Sparsity Methods
function GetJacobianSparsity(jacType::State, pf::PathFunction)
    return pf.stateSP
end
function GetJacobianSparsity(jacType::Control, pf::PathFunction)
    return pf.controlSP
end
function GetJacobianSparsity(jacType::Static, pf::PathFunction)
    return pf.staticSP
end
function GetJacobianSparsity(jacType::Time, pf::PathFunction)
    return pf.timeSP
end


text:
# Common Colorspaces
# ------------------

# The base type
abstract ColorValue{T}
typealias ColourValue{T} ColorValue{T}
abstract AbstractRGB{T} <: ColorValue{T} # allow different memory layouts of RGB
abstract AbstractAlphaColorValue{C <: ColorValue, T <: Real}  # allow different memory layouts of AlphaColorValues

eltype{T}(::ColorValue{T}) = T
eltype{T}(::Type{ColorValue{T}}) = T
eltype{CV<:ColorValue}(::Type{CV}) = eltype(super(CV))

# colortype is defined at the end of this file

# Transparency support
immutable AlphaColorValue{C <: ColorValue, T <: Fractional} <: AbstractAlphaColorValue{C, T}
    c::C
    alpha::T

    function AlphaColorValue(x1::Real, x2::Real, x3::Real, alpha::Real = 1.0)
        new(C(x1, x2, x3), alpha)
    end
    AlphaColorValue(c::ColorValue, alpha::Real) = new(c, alpha)
end
AlphaColorValue{T<:Fractional}(c::ColorValue{T}, alpha::T = one(T)) = AlphaColorValue{typeof(c),T}(c, alpha)

eltype{C<:ColorValue,T}(::AbstractAlphaColorValue{C,T}) = T
eltype{CV<:AbstractAlphaColorValue}(::Type{CV}) = _eltype(CV, super(CV))
_eltype{CV1<:AbstractAlphaColorValue, CV2<:AbstractAlphaColorValue}(::Type{CV1}, ::Type{CV2}) = _eltype(CV2, super(CV2))
_eltype{CV<:ColorValue,T}(::Type{AbstractAlphaColorValue{CV,T}}, ::Type{Any}) =  T

colortype{C}(::AbstractAlphaColorValue{C}) = colortype(C)
colortype{AC<:AbstractAlphaColorValue}(::Type{AC}) = _colortype(AC, super(AC))
_colortype{CV1<:AbstractAlphaColorValue, CV2<:AbstractAlphaColorValue}(::Type{CV1}, ::Type{CV2}) = _colortype(CV2, super(CV2))
_colortype{CV<:ColorValue,T}(::Type{AbstractAlphaColorValue{CV,T}}, ::Type{Any}) =  colortype(CV)

# sRGB (standard Red-Green-Blue)
immutable RGB{T<:Fractional} <: AbstractRGB{T}
    r::T # Red [0,1]
    g::T # Green [0,1]
    b::T # Blue [0,1]

    function RGB(r::Real, g::Real, b::Real)
        new(r, g, b)
    end

end
RGB{T<:Fractional}(r::T, g::T, b::T) = RGB{T}(r, g, b)
RGB(r, g, b) = (T = promote_type(typeof(r), typeof(g), typeof(b)); RGB{T}(r, g, b))
RGB(r::Integer, g::Integer, b::Integer) = RGB{Float64}(r, g, b)
RGB() = RGB(0.0, 0.0, 0.0)

typemin{T}(::Type{RGB{T}}) = RGB{T}(zero(T), zero(T), zero(T))
typemax{T}(::Type{RGB{T}}) = RGB{T}(one(T),  one(T),  one(T))

# HSV (Hue-Saturation-Value)
immutable HSV{T<:FloatingPoint} <: ColorValue{T}
    h::T # Hue in [0,360]
    s::T # Saturation in [0,1]
    v::T # Value in [0,1]

    function HSV(h::Real, s::Real, v::Real)
        new(h, s, v)
    end

end
HSV{T<:FloatingPoint}(h::T, s::T, v::T) = HSV{T}(h, s, v)
HSV(h, s, v) = HSV{Float64}(h, s, v)
HSV() = HSV(0.0, 0.0, 0.0)

HSB(h, s, b) = HSV(h, s, b)

# HSL (Hue-Lightness-Saturation)
immutable HSL{T<:FloatingPoint} <: ColorValue{T}
    h::T # Hue in [0,360]
    s::T # Saturation in [0,1]
    l::T # Lightness in [0,1]

    function HSL(h::Real, s::Real, l::Real)
        new(h, s, l)
    end
end
HSL{T<:FloatingPoint}(h::T, s::T, l::T) = HSL{T}(h, s, l)
HSL(h, s, l) = HSL{Float64}(h, s, l)
HSL() = HSL(0.0, 0.0, 0.0)

HLS(h, l, s) = HSL(h, s, l)

# XYZ (CIE 1931)
immutable XYZ{T<:Fractional} <: ColorValue{T}
    x::T
    y::T
    z::T

    function XYZ(x::Real, y::Real, z::Real)
        new(x, y, z)
    end
end
XYZ{T<:Fractional}(x::T, y::T, z::T) = XYZ{T}(x, y, z)
XYZ(x, y, z) = (T = promote_type(typeof(x), typeof(y), typeof(z)); XYZ{T}(x, y, z))
XYZ(x::Integer, y::Integer, z::Integer) = XYZ{Float64}(x, y, z)
XYZ() = XYZ(0.0, 0.0, 0.0)

# CIE 1931 xyY (chromaticity + luminance) space
immutable xyY{T<:FloatingPoint} <: ColorValue{T}
    x::T
    y::T
    Y::T

    function xyY(x::Real, y::Real, Y::Real)
        new(x, y, Y)
    end
end
xyY{T<:FloatingPoint}(x::T, y::T, Y::T) = xyY{T}(x, y, Y)
xyY(x, y, Y) = xyY{Float64}(x, y, Y)
xyY() = xyY(0.0,0.0,0.0)

# Lab (CIELAB)
immutable Lab{T<:FloatingPoint} <: ColorValue{T}
    l::T # Luminance in approximately [0,100]
    a::T # Red/Green
    b::T # Blue/Yellow

    function Lab(l::Real, a::Real, b::Real)
        new(l, a, b)
    end
end
Lab{T<:FloatingPoint}(l::T, a::T, b::T) = Lab{T}(l, a, b)
Lab(l, a, b) = Lab{Float64}(l, a, b)
Lab() = Lab(0.0, 0.0, 0.0)

typealias LAB Lab

# LCHab (Luminance-Chroma-Hue, Polar-Lab)
immutable LCHab{T<:FloatingPoint} <: ColorValue{T}
    l::T # Luminance in [0,100]
    c::T # Chroma
    h::T # Hue in [0,360]

    function LCHab(l::Real, c::Real, h::Real)
        new(l, c, h)
    end
end
LCHab{T<:FloatingPoint}(l::T, c::T, h::T) = LCHab{T}(l, c, h)
LCHab(l, c, h) = LCHab{Float64}(l, c, h)
LCHab() = LCHab(0.0, 0.0, 0.0)

# Luv (CIELUV)
immutable Luv{T<:FloatingPoint} <: ColorValue{T}
    l::T # Luminance
    u::T # Red/Green
    v::T # Blue/Yellow

    function Luv(l::Real, u::Real, v::Real)
        new(l, u, v)
    end
end
Luv{T<:FloatingPoint}(l::T, u::T, v::T) = Luv{T}(l, u, v)
Luv(l, u, v) = Luv{Float64}(l, u, v)
Luv() = Luv(0.0, 0.0, 0.0)

typealias LUV Luv

# LCHuv (Luminance-Chroma-Hue, Polar-Luv)
immutable LCHuv{T<:FloatingPoint} <: ColorValue{T}
    l::T # Luminance
    c::T # Chroma
    h::T # Hue

    function LCHuv(l::Real, c::Real, h::Real)
        new(l, c, h)
    end
end
LCHuv{T<:FloatingPoint}(l::T, c::T, h::T) = LCHuv{T}(l, c, h)
LCHuv(l, c, h) = LCHuv{Float64}(l, c, h)
LCHuv() = LCHuv(0.0, 0.0, 0.0)

# DIN99 (L99, a99, b99) - adaptation of CIELAB
immutable DIN99{T<:FloatingPoint} <: ColorValue{T}
    l::T # L99
    a::T # a99
    b::T # b99

    function DIN99(l::Real, a::Real, b::Real)
        new(l, a, b)
    end
end
DIN99{T<:FloatingPoint}(l::T, a::T, b::T) = DIN99{T}(l, a, b)
DIN99(l, a, b) = DIN99{Float64}(l, a, b)
DIN99() = DIN99(0.0, 0.0, 0.0)

# DIN99d (L99d, a99d, b99d) - Improvement on DIN99
immutable DIN99d{T<:FloatingPoint} <: ColorValue{T}
    l::T # L99d
    a::T # a99d
    b::T # b99d

    function DIN99d(l::Real, a::Real, b::Real)
        new(l, a, b)
    end
end
DIN99d{T<:FloatingPoint}(l::T, a::T, b::T) = DIN99d{T}(l, a, b)
DIN99d(l, a, b) = DIN99d{Float64}(l, a, b)
DIN99d() = DIN99d(0.0, 0.0, 0.0)

# DIN99o (L99o, a99o, b99o) - adaptation of CIELAB
immutable DIN99o{T<:FloatingPoint} <: ColorValue{T}
    l::T # L99o
    a::T # a99o
    b::T # b99o

    function DIN99o(l::Real, a::Real, b::Real)
        new(l, a, b)
    end
end
DIN99o{T<:FloatingPoint}(l::T, a::T, b::T) = DIN99o{T}(l, a, b)
DIN99o(l, a, b) = DIN99o{Float64}(l, a, b)
DIN99o() = DIN99o(0.0, 0.0, 0.0)

# LMS (Long Medium Short)
immutable LMS{T<:FloatingPoint} <: ColorValue{T}
    l::T # Long
    m::T # Medium
    s::T # Short

    function LMS(l::Real, m::Real, s::Real)
        new(l, m, s)
    end
end
LMS{T<:FloatingPoint}(l::T, m::T, s::T) = LMS{T}(l, m, s)
LMS(l, m, s) = LMS{Float64}(l, m, s)
LMS() = LMS(0.0, 0.0, 0.0)

# 24 bit RGB and 32 bit ARGB (used by Cairo)
# It would be nice to make this a subtype of AbstractRGB, but it doesn't have operations like c.r defined.
immutable RGB24 <: ColorValue{Uint8}
    color::Uint32
end
RGB24() = RGB24(0)
RGB24(r::Uint8, g::Uint8, b::Uint8) = RGB24(convert(UInt32, r)<<16 | convert(UInt32, g)<<8 | convert(UInt32, b))
RGB24(r::Ufixed8, g::Ufixed8, b::Ufixed8) = RGB24(reinterpret(r), reinterpret(g), reinterpret(b))

immutable ARGB32 <: AbstractAlphaColorValue{RGB24, Uint8}
    color::Uint32
end
ARGB32() = ARGB32(0)
ARGB32(r::Uint8, g::Uint8, b::Uint8, alpha::Uint8) = ARGB32(convert(UInt32, alpha)<<24 | convert(UInt32, r)<<16 | convert(UInt32, g)<<8 | convert(UInt32, b))
ARGB32(r::Ufixed8, g::Ufixed8, b::Ufixed8, alpha::Ufixed8) = ARGB32(reinterpret(r), reinterpret(g), reinterpret(b), reinterpret(alpha))

AlphaColorValue(c::RGB24, alpha::Uint8 = 0xff) = AlphaColorValue{typeof(c),Uint8}(c, alpha)

eltype(::RGB24) = Uint8
eltype(::Type{RGB24}) = Uint8
eltype(::ARGB32) = Uint8
eltype(::Type{ARGB32}) = Uint8

# Versions with transparency
typealias RGBA{T} AlphaColorValue{RGB{T},T}
typealias HSVA{T} AlphaColorValue{HSV{T},T}
typealias HSLA{T} AlphaColorValue{HSL{T},T}
typealias XYZA{T} AlphaColorValue{XYZ{T},T}
typealias xyYA{T} AlphaColorValue{xyY{T},T}
typealias LabA{T} AlphaColorValue{Lab{T},T}
typealias LCHabA{T} AlphaColorValue{LCHab{T},T}
typealias LuvA{T} AlphaColorValue{Luv{T},T}
typealias LCHuvA{T} AlphaColorValue{LCHuv{T},T}
typealias DIN99A{T} AlphaColorValue{DIN99{T},T}
typealias DIN99dA{T} AlphaColorValue{DIN99d{T},T}
typealias DIN99oA{T} AlphaColorValue{DIN99o{T},T}
typealias LMSA{T} AlphaColorValue{LMS{T},T}

rgba{T}(c::ColorValue{T}) = AlphaColorValue(convert(RGB{T},c))
hsva{T}(c::ColorValue{T}) = AlphaColorValue(convert(HSV{T},c))
hsla{T}(c::ColorValue{T}) = AlphaColorValue(convert(HSL{T},c))
xyza{T}(c::ColorValue{T}) = AlphaColorValue(convert(XYZ{T},c))
xyYa{T}(c::ColorValue{T}) = AlphaColorValue(convert(xyY{T},c))
laba{T}(c::ColorValue{T}) = AlphaColorValue(convert(Lab{T},c))
lchaba{T}(c::ColorValue{T}) = AlphaColorValue(convert(LCHab{T},c))
luva{T}(c::ColorValue{T}) = AlphaColorValue(convert(Luv{T},c))
lchuva{T}(c::ColorValue{T}) = AlphaColorValue(convert(LCHuv{T},c))
din99a{T}(c::ColorValue{T}) = AlphaColorValue(convert(DIN99{T},c))
din99da{T}(c::ColorValue{T}) = AlphaColorValue(convert(DIN99d{T},c))
din99oa{T}(c::ColorValue{T}) = AlphaColorValue(convert(DIN99o{T},c))
lmsa{T}(c::ColorValue{T}) = AlphaColorValue(convert(LMS{T},c))
argb32{T}(c::ColorValue{T}) = ARGB32(convert(RGB24,c).color | 0xff000000)

const CVconcrete = (HSV, HSL, XYZ, xyY, Lab, Luv, LCHab, LCHuv, DIN99, DIN99d, DIN99o, LMS)
const CVparametric = tuple(RGB, CVconcrete...)
const CVfractional = (RGB, XYZ)
const CVfloatingpoint = (HSV, HSL, xyY, Lab, Luv, LCHab, LCHuv, DIN99, DIN99d, DIN99o, LMS)
const CVAlpha = (RGBA, HSVA, HSLA, XYZA, xyYA, LabA, LCHabA, LuvA, LCHuvA, DIN99A, DIN99dA, DIN99oA, LMSA)

for CV in CVparametric
    @eval begin
        colortype(::$CV) = $CV
        colortype(::Type{$CV}) = $CV
        colortype{T}(::Type{$CV{T}}) = $CV
    end
end

colortype(::Type{RGBA}) = RGB
colortype(::Type{HSVA}) = HSV
colortype(::Type{HSLA}) = HSL
colortype(::Type{XYZA}) = XYZ
colortype(::Type{xyYA}) = xyY
colortype(::Type{LabA}) = Lab
colortype(::Type{LCHabA}) = LCHab
colortype(::Type{LuvA}) = Luv
colortype(::Type{LCHuvA}) = LCHuv
colortype(::Type{DIN99A}) = DIN99
colortype(::Type{DIN99dA}) = DIN99d
colortype(::Type{DIN99oA}) = DIN99o
colortype(::Type{LMSA}) = LMS

# Vector space operations
import Base: +, *

#XYZ is a linear vector space
+{T<:Number}(a::XYZ{T}, b::XYZ{T}) = XYZ(a.x+b.x, a.y+b.y, a.z+b.z)
*(c::Number, a::XYZ) = XYZ(c*a.x, c*a.y, c*a.z)

#Most color spaces are nonlinear, so do the arithmetic in XYZ and convert back
+{T<:ColorValue}(a::T, b::T) = convert(T, convert(XYZ, a) + convert(XYZ, b))
*{T<:ColorValue}(c::Number, a::T) = convert(T, c * convert(XYZ, a))


text:
# This file generates the plot in the logo
using Ripserer
using Plots
gr()
using Random
Random.seed!(1337)

function circle(n; r=0.7, center, s=0)
    points = NTuple{2,Float64}[]
    x, y = center
    for θ in range(0, 2π; length=n + 1)[circshift(1:n, s)]
        push!(points, (x + r * sin(θ), y + r * cos(θ)))
    end
    return points
end

data = [
    circle(20; center=(-1, 0), s=0)
    circle(20; center=(1, 0), s=0)
    circle(20; center=(0, √3), s=0)
]

res = ripserer(data; reps=true, dim_max=2)

ripserer_logo = plot(;
    aspect_ratio=1,
    xlim=(-2.1, 2.1),
    ylim=(-1.1, 3.1),
    xlabel=:x,
    ylabel=:y,
    ticks=nothing,
    legend=false,
    border=:none,
)
for (i, int) in enumerate(sort(res[2]; by=death, rev=true))
    if i == 1
        col = 2 # red
        α = 1
    elseif i == 2
        col = 4 # purple
        α = 1
    elseif i == 3
        col = 3 # green
        α = 1
    else
        col = 1
        α = 0.2
    end

    plot!(
        ripserer_logo,
        int,
        data;
        alpha=α,
        color=col,
        threshold_strict=death(int),
        linewidth=5,
    )
end
scatter!(ripserer_logo, data; markersize=4, color=1)
display(ripserer_logo)

diagrams_logo = plot(
    res;
    title="",
    ticks=nothing,
    legend=false,
    xlabel="",
    ylabel="",
    markersize=3,
    size=(100, 100),
)


text:
"""

# `syms_in_paths`

Collect vertices in all paths.

$(SIGNATURES)

Part of the API, Exported
"""
function syms_in_paths(paths, f, l)
  thepaths = deepcopy(paths)
  syms = Symbol[]
  for p in thepaths
    setdiff!(p, [f, l])
    append!(syms, p)
    unique!(syms)
  end
  syms
end

"""

# `syms_in_all_paths`

Check if a vertex is part of all paths.

$(SIGNATURES)

Part of the API, Exported
"""
function sym_in_all_paths(paths, sym)
  all([sym in p for p in paths])
end

"""

# `adjustment_sets`

$(SIGNATURES)

Computes the covariance adjustment vertex set.

### Required arguments
```julia
* `dag::DAG`                           : DAG
* `f::Symbol`                          : Start variable
* `l::Symbol`                          : End variable
```

### Optional arguments
```julia
* `debug::Bool`                        : Show debug trace
```

### Returns
```julia
* `adjustmentsets=Vector{Symbol}[]`    : Array of adjustment sets
```

# Extended help

### Acknowledgements

Original author                        : Rob J Goedman

### Licence

Licenced under: MIT.

Part of the api, exported.
"""
function adjustment_sets(dag::DAG, f::Symbol, l::Symbol; debug=false)

  ap  = all_paths(dag, f, l)
  paths = backdoor_paths(dag, ap, f)
  debug && display(paths)
  lsyms = syms_in_paths(paths, f, l)
  adjustmentsets = Vector{Symbol}[]
  for s in lsyms
    debug && println("checking $s in $lsyms")
    if sym_in_all_paths(paths, s)
      opensets = length(open_paths(dag, paths, [s]))
      if  opensets == 0
        debug && println("$s closes all paths.")
        append!(adjustmentsets, [[s]])
      end
    elseif length(open_paths(dag, paths, [s])) == 0
        debug && println("$s not in all paths but no open paths.")
        append!(adjustmentsets, [[s]])
    end
  end
  for s in adjustmentsets
    setdiff!(lsyms, s)
  end
  debug && lsyms |> display
  len = 2
  local csyms
  while length(lsyms) >= len
    csyms = collect(combinations(lsyms, len))
    debug && csyms |> display
    for s in csyms
      if length(open_paths(dag, paths, s)) == 0
        debug && println("$s closes all paths")
        append!(adjustmentsets, [s])
        setdiff!(lsyms, s)
      end
    end
    len += 1
  end
  debug && lsyms |> display

  adjustmentsets
end

export
  adjustment_sets,
  syms_in_paths,
  sym_in_all_paths


text:
using RDatasets, DataArrays, DataFrames, Gadfly

barley = dataset("lattice", "barley")
if isdefined(:setlevels!)
    setlevels!(barley[:Year], ["1931", "1932"])
else
    set_levels!(barley[:Year], ["1931", "1932"])
end

plot(barley,
     xgroup="Variety", ygroup="Site", x="Year", y="Yield",
     Geom.subplot_grid(Geom.line, Geom.point,
                       free_x_axis=true, free_y_axis=true))


text:
using CPUTime
using CSV
using DataFrames
using Dates

function do_example_1(path, threshold::T) where T

    now = Dates.now(Dates.UTC)

    cd(path)  # pathに移動
    file_name_all = readdir()  # カレントディレクトリにある全ファイル名取得
    csv_name = file_name_all[
        map(x->occursin(r"\.csv$",x),file_name_all)
    ]  # csvの名前だけ取得

    result = []

    # 条件を満たすcsvファイルを探す
    for name in csv_name
        df = CSV.read(name, DataFrame)
        if size(df[df.x .> threshold, :], 1) > 0
            push!(result, name)
        end
    end

    cd("..")
    if isdir("./example_rsult")
        ##
    else
        mkdir("./example_rsult")
    end
    cd("./example_rsult")

    # 結果をtxtで出力
    dt = Dates.format(now, "yyyy-mm-dd--HH-MM-SS")
    result_name = string(dt) * "_result.txt"

    open(result_name, "w+") do f
        print(f, "result is")
        for name in result
            println(f, name)
        end
    end

    cd("..")
    cd("..")
end

path = "./example_jl/example_1_data"  # csvフォルダのパス
threshold = 40.0  # xの閾値

@time do_example_1(path, threshold)


text:
using HullGMP
using Test
using LazySets

dnn1 = read_nnet("nnet/dnn1.nnet")
dnn2 = read_nnet("nnet/dnn2.nnet")
dnn3 = read_nnet("nnet/dnn3.nnet")
dnn4 = read_nnet("nnet/dnn4.nnet")
dnn5 = read_nnet("nnet/dnn5.nnet")

delta = 0.1

solver1 = MaxSens(delta, true)

solver2 = HullReach(delta, true)

inputSet1 = Hyperrectangle(fill(0.5, 2), fill(0.5, 2))
inputSet2 = Hyperrectangle(fill(0.5, 3), fill(0.5, 3))
inputSet3 = Hyperrectangle(fill(0.5, 4), fill(0.5, 4))
inputSet4 = Hyperrectangle(fill(0.5, 5), fill(0.5, 5))
inputSet5 = Hyperrectangle(fill(0.5, 6), fill(0.5, 6))

output_center = fill(0.0, 2)
output_radius = fill(3.0, 2)
outputSet = Hyperrectangle(output_center, output_radius) #output constraint

problem1 = Problem(dnn1, inputSet1, outputSet)
problem2 = Problem(dnn2, inputSet2, outputSet)
problem3 = Problem(dnn3, inputSet3, outputSet)
problem4 = Problem(dnn4, inputSet4, outputSet)
problem5 = Problem(dnn5, inputSet5, outputSet)

file = open("group1_5_dnns.txt", "a")
print(file, "Test Result of group1_5_dnns:\n\n")

#solver11
time11 = 0
solve(solver1, problem1)
for i = 1:10
    timed_result =@timed solve(solver1, problem1)
    print(file, "MaxSens - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time11 += timed_result.time
end
print(file, "dnn1 Average time: " * string(time11/10) * " s\n\n")

#solver12
time12 = 0
solve(solver1, problem2)
for i = 1:10
    timed_result =@timed solve(solver1, problem2)
    print(file, "MaxSens - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time12 += timed_result.time
end
print(file, "dnn2 Average time: " * string(time12/10) * " s\n\n")

#solver13
time13 = 0
solve(solver1, problem3)
for i = 1:10
    timed_result =@timed solve(solver1, problem3)
    print(file, "MaxSens - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time13 += timed_result.time
end
print(file, "dnn3 Average time: " * string(time13/10) * " s\n\n")

#solver14
time14 = 0
solve(solver1, problem4)
for i = 1:10
    timed_result =@timed solve(solver1, problem4)
    print(file, "MaxSens - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time14 += timed_result.time
end
print(file, "dnn4 Average time: " * string(time14/10) * " s\n\n")

#solver15
time15 = 0
solve(solver1, problem5)
for i = 1:10
    timed_result =@timed solve(solver1, problem5)
    print(file, "MaxSens - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time15 += timed_result.time
end
print(file, "dnn5 Average time: " * string(time15/10) * " s\n\n")

#solver21
time21 = 0
solve(solver2, problem1)
for i = 1:10
    timed_result =@timed solve(solver2, problem1)
    print(file, "HullReach - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time21 += timed_result.time
end
print(file, "dnn1 Average time: " * string(time21/10) * " s\n\n")

#solver22
time22 = 0
solve(solver2, problem2)
for i = 1:10
    timed_result =@timed solve(solver2, problem2)
    print(file, "HullReach - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time22 += timed_result.time
end
print(file, "dnn2 Average time: " * string(time22/10) * " s\n\n")

#solver23
time23 = 0
solve(solver2, problem3)
for i = 1:10
    timed_result =@timed solve(solver2, problem3)
    print(file, "HullReach - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time23 += timed_result.time
end
print(file, "dnn3 Average time: " * string(time23/10) * " s\n\n")

#solver24
time24 = 0
solve(solver2, problem4)
for i = 1:10
    timed_result =@timed solve(solver2, problem4)
    print(file, "HullReach - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time24 += timed_result.time
end
print(file, "dnn4 Average time: " * string(time24/10) * " s\n\n")

#solver25
time25 = 0
solve(solver2, problem5)
for i = 1:10
    timed_result =@timed solve(solver2, problem5)
    print(file, "HullReach - test " * string(i) * " - Time: " * string(timed_result.time) * " s")
    print(file, " - Output: " * string(timed_result.value) * "\n")
    global time25 += timed_result.time
end
print(file, "dnn5 Average time: " * string(time25/10) * " s\n\n")

close(file)


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

#Various sizes of hexagonal clusters that are useful in lenslet design. To see what they look like use the functions in HexTilings to make lattices with color properties and then draw them with Vis.draw
#example: Vis.draw(hex3RGB())

function hex3(scale::T =  1.0) where{T<:Real}
    clusterelements = SVector((0, 0), (-1, 0), (-1, 1))
    eltlattice = HexBasis1(scale)
    clusterbasis = LatticeBasis((-1, 2), (2, -1))
    return LatticeCluster(clusterbasis, eltlattice, clusterelements)

end
export hex3

function hex4(scale::T =  1.0) where{T<:Real}
    clusterelements = SVector((-1, 0), (-1, 1), (0, 0), (0, 1))
    eltlattice = HexBasis1(scale)
    clusterbasis = LatticeBasis((2, 0), (0, 2))
    return LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex4

hex7(scale::T = 1.0) where{T<:Real} = hexn(1,scale)
export hex7

function hex9(scale::T =  1.0) where{T<:Real}
    clusterelements = SVector((-1, -1), (-1, 0), (-2, 1), (0, -1), (0, 0), (-1, 1), (1, -1), (1, 0), (0, 1))
    eltlattice = HexBasis1(scale)
    clusterbasis = LatticeBasis((3, 0), (0, 3))
    return LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex9

function hex12(scale::T =  1.0) where{T<:Real}
    clusterelements = SVector(
        (-1, 1),(0, 1),
        (-1, 0),(0, 0),(1, 0),
        (-1, -1),(0, -1),(1, -1),(2, -1),
        (0, -2),(1, -2),(2, -2)
    )

    eltlattice = HexBasis3(scale)
    clusterbasis = LatticeBasis((2, 2), (-3, 2))
    return LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex12

#TODO fix cluster basis vectors.
"""WARNING: the cluster basis vectors are incorrect. Still need to figure out what these are. This will not tile properly"""
function hex18(scale::T =  1.0) where{T<:Real}
    throw(ErrorException("not yet working"))
    clusterelements = SVector(
        (-2,1),(-1, 1),(0, 1),(1,1),
        (-2,0),(-1, 0),(0, 0),(1, 0),(2,0),
        (-1, -1),(0, -1),(1, -1),(2, -1),
        (0, -2),(1, -2),(2, -2),
        (1,-3),(2,-3)
    )

    eltlattice = HexBasis3(scale)
    clusterbasis = LatticeBasis((2, 2), (-3, 2))
    return LatticeCluster(clusterbasis, eltlattice, clusterelements)
end
export hex18

hex19(scale::T = 1.0) where{T<:Real} = hexn(2,scale)
export hex19

hex37(scale::T = 1.0) where{T<:Real} = hexn(3,scale)
export hex37

""" function for creating symmetric clusters consisting of all lattice cells within regionsize of the origin. This includes hex1 (use regionsize = 0),hex7 (regionsize = 1),hex19 (region size = 2),hex37 (regionsize = 3), etc."""
function hexn(regionsize::Int64,scale::T = 1.0) where{T<:Real}
    eltlattice = HexBasis1(scale)
    clusterbasis = LatticeBasis((2*regionsize+1, -regionsize), (regionsize, regionsize+1))
    return LatticeCluster(clusterbasis, eltlattice, region(HexBasis1,(0,0),regionsize))
end
export hexn

function testhexclusters()
    for clusterfunc in (hex3,hex4,hex9,hex12,hex19,hex37)
        cluster = clusterfunc()
        println("$clusterfunc euclidean diameter $(euclideandiameter(cluster))")
    end
end
export testhexclusters


text:
using Roots
using Distributions
using NLsolve
using QuadGK
#using ForwardDiff
#using SpecialFunctions
#include("cr_sod_shock_main.jl")
#include("beta_inc_functions.jl")

mutable struct SodCRParameters_withCRs

    rhol::Float64
    rhor::Float64
    Pl::Float64
    Pr::Float64
    Ul::Float64
    Ur::Float64
    P_cr_l::Float64
    P_cr_r::Float64
    E_cr_l::Float64
    E_cr_r::Float64
    cl::Float64
    cr::Float64
    M::Float64
    t::Float64
    x_contact::Float64
    Pe_ratio::Float64
    γ_th::Float64
    γ_cr::Float64
    Δγ::Float64
    γ_exp::Float64
    α::Float64
    β::Float64
    η2::Float64
    # acc_function#::F
    # reacc_function

    function SodCRParameters_withCRs(;rhol::Float64=1.0,  rhor::Float64=0.125,
                                Pl::Float64=0.0,    Pr::Float64=0.0,
                                Ul::Float64=0.0,    Ur::Float64=0.0,
                                P_cr_l::Float64=0.0,      P_cr_r::Float64=0.0,
                                E_cr_l::Float64=0.0,      E_cr_r::Float64=0.0,
                                Mach::Float64=0.0,  t::Float64,
                                x_contact::Float64=70.0,
                                Pe_ratio::Float64=0.01,
                                γ_th::Float64=5.0/3.0,
                                γ_cr::Float64=4.0/3.0,
                                eff_model::Int64=-1)

        γ_exp    = ( γ_th - 1.0 )/( 2.0 * γ_th )
        η2       = (γ_th-1.0)/(γ_th+1.0)
        Δγ       = γ_th - γ_cr

        # calculate Ul and Pl depending on input
        if (Pl == 0.0) & (Ul != 0.0)
            Pl = ( γ_th - 1.0 ) * rhol * Ul
        elseif (Ul == 0.0) & (Pl != 0.0)
            Ul = Pl / ( (γ_th - 1.0) * rhol )
        else
            error("Both Ul and Pl are zero!")
        end

        # calculate B angle dependent efficiency following Pais+ 2018, MNRAS, 478, 5278
        # delta_theta = π/18.0
        # thetaB *= (π/180)
        # etaB = 0.5*( tanh( (theta_crit - thetaB)/delta_theta ) + 1.0 )

        if eff_model == -1
            acc_function = null_acc
        elseif eff_model == 0
            acc_function = KR07
        elseif eff_model == 1
            acc_function = KR13
        elseif eff_model == 2
            acc_function = R_19
        elseif eff_model == 3
            acc_function = CS15
        else
            error("Invalid DSA model selection!\n
                   Pick one of the available models, or solve a pure Hydro shock with:\n
                   SodParameters")

        end

        # calculate Ur and Pr depending on input
        if (Pr == 0.0) & (Ur != 0.0)
            Pr = ( γ_th - 1.0 ) * rhor * Ur
        elseif (Ur == 0.0) & (Pr != 0.0)
            Ur = Pr / ( (γ_th - 1.0) * rhor )
        elseif (Ur == 0.0) & (Pr == 0.0) & (Mach == 0.0)
            println("Error! Ur, Pr and Mach are zero! Can't find solution!")
        else
            println("Both Ur and Pr are zero! Will calculate them depending on Machnumber.")
            Pr = solvePrfromMach(rhol, rhor, Pl, Mach, γ_th, γ_cr, eff_function)
            Ur = Pr / ( (γ_th - 1.0) * rhor )
        end
        #
        # if Mach == 0.0
        #     Mach = solveMach(Pl, Pr, rhol, rhor, γ_th)
        # end

        cl = sqrt.(γ_eff(Pl, P_cr_l, γ_th, γ_cr) * (Pl + P_cr_l) / rhol)
        cr = sqrt.(γ_eff(Pr, P_cr_r, γ_th, γ_cr) * (Pr + P_cr_r) / rhor)

        α = (γ_cr - 1.0)/(2.0 * Δγ)
        β = (1.0 - γ_th)/(2.0 * Δγ)

        E_cr_r = P_cr_r/( ( γ_cr - 1.0 )*rhor )

        new(rhol, rhor,
            Pl, Pr,
            P_cr_l, P_cr_r,
            Ul, Ur,
            E_cr_r,
            cl, cr,
            Mach, t,
            x_contact,
            Pe_ratio,
            γ_th, γ_cr,
            Δγ,
            γ_exp,
            α, β,
            η2)#,
            #acc_function,
            #acc_function)
    end
end

"""
        Integral function
"""

@inline A(P::Float64, ρ::Float64, γ::Float64) = γ * P * ρ^(-γ)
@inline a(P::Float64, γ::Float64) = γ * P

@inline function x(ρ::Float64, P_th::Float64, P_cr::Float64, par::SodCRParameters_withCRs)
    return  a(P_th, par.γ_th)/(a(P_th, par.γ_th) + a(P_cr, par.γ_cr))
end

@inline function reduced_Beta(x, α, β)
    return 1.0/beta(α, β) * x^(α - 1.0) * ( 1.0 - x )^(β - 1.0)
end

@inline function γ_eff(P_th::Float64, P_cr::Float64, γ_th::Float64, γ_cr::Float64)
    return (γ_cr*P_cr + γ_th*P_th)/(P_th + P_cr)
end

@inline function incomplete_beta(a::Float64, b::Float64, t::Float64)
    return t^(a-1.0) * (1.0-t)^(b-1.0)
end

@inline function I(rho::Float64, P_th::Float64, P_cr::Float64, par::SodCRParameters_withCRs)

    x_ = x(rho, P_th, P_cr,  par)
    A_cr = A(P_cr, rho, par.γ_cr)
    A_th = A(P_th, rho, par.γ_th)
    B(x) = incomplete_beta(par.α, par.β, x)

    result_B = quadgk(B, 0, x_, rtol=1.e-4)

    return √(A_cr)/par.Δγ * (A_cr/A_th)^par.α * result_B
end

"""
    Pressure
"""
function solveP3(par::SodCRParameters_withCRs, sol::SodCRSolution)
    # Solves the pressure of the middle state in a shocktube

    Pcr3 = par.P_cr_l * ( sol.rho3/par.rhol )^par.γ_cr
    P3 = par.Pl * ( sol.rho3/par.rhol )^par.γ_th

    return P3, Pcr3
end

function solveP2(x::Float64; par::SodCRParameters_withCRs)
    # solves the pressure along the rarefaction wave
    return par.Pl * ( -par.η2 * x / ( par.cl * par.t ) + ( 1.0 - par.η2 ) )^( 1.0/par.γ_exp )
end

function solveP(x::Float64, par::SodCRParameters_withCRs, sol::SodCRSolution)
    # returns P values according to position in shocktube

    if x <= -par.cl * par.t
        return (par.Pl + par.P_cr_l),
                par.Pl, par.P_cr_l,
                (1.0 - par.Pe_ratio)*par.P_cr_l,
                par.Pe_ratio*par.P_cr_l

    elseif -par.cl * par.t < x <= -sol.vt * par.t
        return 0.0, 0.0, 0.0, 0.0, 0.0

    elseif -sol.vt * par.t < x <= sol.v34 * par.t
        return (sol.P3_th + sol.P3_cr),
                sol.P3_th, sol.P3_cr,
                (1.0 - par.Pe_ratio)*sol.P3_cr,
                par.Pe_ratio*sol.P3_cr

    elseif sol.v34 * par.t < x <= sol.vs * par.t
        return (sol.P4_th + sol.P4_cr),
                sol.P4_th, sol.P4_cr,
                (1.0 - par.Pe_ratio)*par.P4_cr,
                par.Pe_ratio*par.P4_cr

    elseif sol.vs * par.t < x
        return (par.Pr + par.P_cr_r),
                par.Pr, par.P_cr_r,
                (1.0 - par.Pe_ratio)*par.P_cr_r,
                par.Pe_ratio*par.P_cr_r
    else
        println("Error! x = $x  t = $(par.t)")
        return 0.0, 0.0, 0.0, 0.0, 0.0
    end
end

"""
    Density
"""
function solveRho2(x::Float64, par::SodCRParameters_withCRs)
    # solves density along the rarefaction wave

    #f(rho) = I(rho) - I(par.rhol) + x/par.t + sqrt.( A())
    return 0.0

end

function rho34_solver!(F, x, par::SodCRParameters_withCRs)

    xs(rho4)        = rho4/par.rhor
    xr(rho3)        = rho3/par.rhol

    P34(rho3)       = par.P_cr_l*xr(rho3)^par.γ_cr + par.Pl*xr(rho3)^par.γ_th

    Pcr4(rho4)      = par.P_cr_r*xs(rho4)^par.γ_cr

    ε2(rho4, rho3)  = par.E_cr_r*xs(rho4)^par.γ_cr + 1.0/(par.γ_th - 1.0) *
                     ( P34(rho3) - Pcr4(rho4) )

    Pth3(rho3)      = par.Pl * xr(rho3)^par.γ_th
    Pcr3(rho3)      = par.P_cr_l * xs(rho3)^par.γ_cr

    F[1] = (P34(x[1]) - (par.Pr + par.P_cr_r)) * ( xs(x[2]) - 1.0 ) -
            par.rhor * xs(x[2]) * ( I(par.rhol, par.Pl, par.P_cr_l, par) - I(x[1], Pth3(x[1]), Pcr3(x[1]), par) )^2

    F[2] = (P34(x[1]) + (par.Pr + par.P_cr_r)) * ( xs(x[2]) - 1.0 ) +
            2.0 * ( xs(x[2]) * (par.Ur + par.E_cr_r) - ε2(x[2], x[1]) )

end

function solveRho34(par::SodCRParameters_withCRs)

    rho34_helper!(F, x) = rho34_solver!(F, x, par)

    initial_x = [par.rhol, par.rhor]
    nlsolve(rho34_helper!, initial_x)

        #   rho3           rho4
    return initial_x[1], initial_x[2]

end

function solveRho(x::Float64, par::SodCRParameters_withCRs, sol::SodCRSolution)
    # returns P values according to position in shocktube

    if x <= -par.cl * par.t
        return par.rhol
    elseif -par.cl * par.t < x <= -sol.vt * par.t
        return solveRho2(x, par)
    elseif -sol.vt * par.t < x <= sol.v34 * par.t
        return sol.rho3
    elseif sol.v34 * par.t < x <= sol.vs * par.t
        return sol.rho4
    elseif sol.vs * par.t < x
        return par.rhor
    else
        println("Error!")
        return 0.0
    end
end

"""
    Velocity
"""
function solveV34(par::SodCRParameters_withCRs, sol::SodCRSolution)
    # solves velocity along the isopressure region 2-3
    return sqrt.( (sol.P34_tot - (par.Pr + par.P_cr_r ) ) * (sol.rho4 - par.rhor)/sol.rho4*par.rhor )
end

function solveV2(x::Float64, par::SodCRParameters_withCRs)
    # solves the velocity along the rarefaction wave
    return 0.0
end

function solveVs(par::SodCRParameters_withCRs, sol::SodCRSolution)
    # solves the shock velocity
    return sol.v34 * sol.rho4 / ( sol.rho4 - par.rhor )
end

function solveVt(par::SodCRParameters_withCRs, sol::SodCRSolution)
    # solves the velocity of the tail of the rarefaction wave
    return I(sol.rho3, sol.P3_th, sol.P3_cr, par) - I(par.rhol, par.Pl, par.P_cr_l, par) +
            sqrt.(
                A(sol.P3_cr, sol.rho3, par.γ_cr)*sol.rho3^(par.γ_cr - 1.0) +
                A(sol.P3_th, sol.rho3, par.γ_th)*sol.rho3^(par.γ_th - 1.0) )
end

function solveV(x::Float64, par::SodCRParameters_withCRs, sol::SodCRSolution)
    # solves the velocity along the shocktube
    if x <= -par.cl * par.t
        return 0.0
    elseif -par.cl * par.t < x <= -sol.vt * par.t
        return solveV2(x, par)
    elseif -sol.vt * par.t < x <= sol.vs * par.t
        return sol.v34
    elseif sol.vs * par.t < x
        return 0.0
    else
        println("Error!")
        return 0.0
    end
end

"""
    Main function for solver
"""

function solveSodShockCR_withPrepopulation(x::Array{Float64,1}; par::SodCRParameters_withCRs)

    # set up datatype to store riemann solution
    sol = SodCRSolution(x)

    # transform into rest-frame of contact discontiuity
    x_in = sol.x .- par.x_contact

    # solve Pressure
    sol.rho3, sol.rho4 = solveRho34(par)

    sol.xs  = sol.rho4/par.rhor
    sol.xr  = sol.rho3/par.rhol

    sol.P34_tot = par.P_cr_l*sol.xr^par.γ_cr + par.Pl*sol.xr^par.γ_th

    sol.P3_cr = par.P_cr_l*(sol.rho3/par.rhol)^par.γ_cr
    sol.P3_th = sol.P34_tot - sol.P3_cr

    sol.P4_cr = par.P_cr_r*sol.xs^par.γ_cr
    sol.P4_th = sol.P34_tot - sol.P4_cr

    # solve velocity
    sol.v34 = solveV34(par, sol)

    sol.vs = solveVs(par, sol)
    sol.vt = solveVt(par, sol)

    for i=1:length(sol.x)
        sol.v[i] = solveV(x_in[i], par, sol)
    end

    for i = 1:length(sol.x)
        sol.P_tot[i], sol.P_th[i], sol.P_cr[i], sol.P_cr_p[i], sol.P_cr_e[i] = solveP(x_in[i], par, sol)
    end

    for i = 1:length(sol.x)
        sol.rho[i] = solveRho(x_in[i], par, sol)
    end

    sol.Mach = sol.vs/par.cr

    sol.U_tot = sol.P_tot ./ ( (par.γ_th - 1.0) .* sol.rho )
    sol.U_th = sol.P_th ./ ( (par.γ_th - 1.0) .* sol.rho )
    sol.E_cr_p = sol.P_cr_p ./ ( (par.γ_cr - 1.0) .* sol.rho )
    sol.E_cr_e = sol.P_cr_e ./ ( (par.γ_cr - 1.0) .* sol.rho )

    return sol
end

# x_in = collect(50.0:0.01:100.0)
#
# Pl = 63.499
# Pr = 0.1
# Pcrl = 0.3*Pl
# Pcrr = 0.3*Pr
# par = SodCRParameters_withCRs(Pl=Pl, Pr=Pr, P_cr_l=Pcrl, P_cr_r=Pcrr, t=1.5, eff_model=-1)
#
# sol = solveSodShockCR_withPrepopulation(x_in, par=par)
#
# B(x) = incomplete_beta(par.α, par.β, x)
#
# result_B = quadgk(B, 0, 4.74, rtol=1.e-4)

max_stars_repo_name: LudwigBoess/GadJet

max_issues_repo_name: LudwigBoess/GadJet

max_forks_repo_name: LudwigBoess/GadJet


text:
# Statistics
function Statistics.mean(f::Function, sim::EPOCHSimulation; cond=x->true)
    ThreadsX.map(sim.files) do file
        println("Loading $(file.name)")
        qunatity = f(file)
        z = zero(eltype(qunatity))
        (qunatity ./ length(qunatity)) |>
            Filter(cond) |>
            foldxt(+, simd=true, init=z)
    end
end

# FileTrees
# FileTrees._maketree(node::SDFFile) = File(nothing, basename(node.name), node)

pretty_summarysize(x) = Base.format_bytes(Base.summarysize(x))

get_parameter(sim::EPOCHSimulation, p::Symbol) = getindex(sim.param, p)
get_parameter(sim::EPOCHSimulation, p::Symbol, c::Symbol) = getindex(get_parameter(sim, p), c)

function domain_length(sdf, dir::Symbol)
    d_min = Symbol(string(dir)*"_min")
    d_max = Symbol(string(dir)*"_max")

    get_parameter(sdf, d_max) - get_parameter(sdf, d_min)
end

function domain_volume(sdf)
    nd = ndims(sdf)

    Δx = domain_length(sdf, :x)

    if nd == 1
        Δx
    elseif nd == 2
        Δy = domain_length(sdf, :y)
        Δx * Δy
    else
        Δy = domain_length(sdf, :y)
        Δz = domain_length(sdf, :z)
        Δx * Δy * Δz
    end
end

function cell_volume(sdf)
    nd = ndims(sdf)
    vol = domain_volume(sdf)
    nx = get_parameter(sdf, :nx)

    if nd == 1
        vol / nx
    elseif nd == 2
        ny = get_parameter(sdf, :ny)
        vol / (nx * ny)
    else
        ny = get_parameter(sdf, :ny)
        nz = get_parameter(sdf, :nz)
        vol / (nx * ny * nz)
    end
end

function cell_length(sdf, direction::Symbol)
    domain_length(sdf, direction)/get_parameter(sdf, Symbol("n"*"$direction"))
end

function timestep(sim)
    if haskey(sim, :dt_multiplier)
        C = get_parameter(sim, :dt_multiplier)
    else
        C = 0.95
    end

    n = ndims(sim)
    c = 2.99792458e8u"m/s" # https://github.com/Warwick-Plasma/epoch/blob/066d3dafd5582360c2e9136d2020cef404bb84f5/epoch1d/src/constants.F90#L147

    Δx = cell_length(sim, :x)

    if n == 1
        Δt = Δx / c
    elseif n == 2
        Δy = cell_length(sim, :y)
        Δt = inv(c * √(1/Δx^2 + 1/Δy^2))
    else
        Δy = cell_length(sim, :y)
        Δz = cell_length(sim, :z)

        Δt = inv(c * √(1/Δx^2 + 1/Δy^2 + 1/Δz^2))
    end

    return C * Δt
end


text:
using Test, YaoBlocks, YaoArrayRegister

reg = rand_state(10)
block = kron(4, 2=>X)
c = concentrate(10, block, [1, 3, 9, 2])
@test nqubits(c) == 10
@test nactive(c) == 4
@test isunitary(c) == true
@test isreflexive(c) == true
@test ishermitian(c) == true
blk = kron(4, 2=>Rx(0.3))
@test chsubblocks(c, [blk]) |> subblocks |> first == blk

@test apply!(copy(reg), c) == apply!(copy(reg), kron(10, 3=>X))
@test apply!(rand_state(12, nbatch=10) |> focus!(Tuple(1:10)...), c) |> nactive == 10

@testset "test repeat" begin
    c = concentrate(8, repeat(5, H), 1:5)
    r = rand_state(8)
    r1 = copy(r) |> c
    r2 = copy(r) |> repeat(8, H, 1:5)
    @test r1 ≈ r2
end

@testset "mat" begin
    cc = concentrate(5, kron(X,X), (1,3))
    @test applymatrix(cc) ≈ mat(cc)
end


text:
using DelimitedFiles, PyPlot, Statistics
#pyplot()

labels = ["Planet b","c","d","e","f","g","h"]
color = ["C0","C1","C2","C3","C4","C5","C6"]
nbg = readdlm("nbg_tts.txt", ',')[:,1:3000]
reb = readdlm("reb_tts.txt", ',')[:,1:3000]
fast = readdlm("ttvfast_tts.txt", ',')[:,1:3000]

fig,axes = subplots(2,2,figsize=(10,5))
reb_diff = nbg .- reb
#plot()
ax = axes[1]
for i in 1:3
    ind = findfirst(isequal(0.0),reb[i,:])
    ind = ind === nothing ? ind = 3000 : ind = ind - 5
    tt = nbg[i,1:ind] ; pmean = mean(tt[2:ind] .- tt[1:ind-1]);
    iref = collect(0:1:ind-1)
    println(labels[i]," ",ind," ",maximum(tt))
    tref = mean(tt .- pmean .* iref) .+ pmean .* iref ; ttv = (tt .- tref) .* 24*60
    ax.plot(nbg[i,1:ind], ttv, label=labels[i],color=color[i])
    ttv = (reb[i,1:ind] .- tref) .* 24*60
    ax.plot(nbg[i,1:ind], ttv, color="k",linestyle=":")
end
ax.legend()
ax.set_title("TTV of inner 3 planets")
ax.set_xlabel("Time [Days]")
ax.set_ylabel("TTV [minutes]")
ax.axis([0,4000,-800,800])
ax = axes[2]
for i in 4:7
    ind = findfirst(isequal(0.0),reb[i,:])
    ind = ind === nothing ? ind = 3000 : ind = ind - 5
    tt = nbg[i,1:ind] ; pmean = mean(tt[2:ind] .- tt[1:ind-1]);
    iref = collect(0:1:ind-1)
    tref = mean(tt .- pmean .* iref) .+ pmean .* iref ; ttv = (tt .- tref) .* 24*60
    ax.plot(nbg[i,1:ind], ttv, label=labels[i],color=color[i])
    ttv = (fast[i,1:ind] .- tref) .* 24*60
    ax.plot(nbg[i,1:ind], ttv, color="k",linestyle=":")
end
ax.set_title("TTV of outer 4 planets")
ax.set_xlabel("Time [Days]")
ax.set_ylabel("TTV [minutes]")
ax.legend(ncol=2)
ax.axis([0,4000,-800,800])

ax = axes[4]
for i in 1:7
    ind = findfirst(isequal(0.0),reb[i,:])
    ind = ind === nothing ? ind = 3000 : ind = ind - 5
    ax.plot(nbg[i,1:ind], reb_diff[i,1:ind] .* 24*60*60, label=labels[i],color=color[i],".")
end

ax.set_title("Time difference with REBOUND")
ax.set_ylabel("Transit time diff [sec]")
ax.set_xlabel("Time [Days]")
ax.axis([0,4000,-0.000004,0.000004])
#ax.legend(ncol=2,loc="upper left")
#read(stdin,Char)
#savefig("nbg_vs_reb_transit_times.pdf")

fast_diff = nbg .- fast
#plot()
ax = axes[3]
for i in 1:7
    ind = findfirst(isequal(0.0),fast[i,:])
    ind = ind === nothing ? ind = 3000 : ind = ind - 5
    ax.plot(nbg[i,1:ind], fast_diff[i,1:ind] .* 24*60*60, label=labels[i],color=color[i])
end
ax.set_title("Time difference with TTVFast")
ax.set_ylabel("Transit time diff [sec]")
ax.set_xlabel("Time [Days]")
ax.axis([0,4000,-0.008,0.008])
#ax.legend(ncol=2,loc="upper left")
#savefig("nbg_vs_ttvfast_transit_times.pdf")

fig.tight_layout()
savefig("nbg_vs_ttvfast_vs_rebound_transit_times.pdf")

max_forks_repo_name: langfzac/nbg-papers


text:
using Documenter, Diana

makedocs(format = Documenter.HTML(),
         modules = [Diana],
         sitename = "Diana.jl"
         )
deploydocs(
	repo = "github.com/codeneomatrix/Diana.jl.git"
	#devbranch = "master",
    #devurl = "dev",
    #versions = ["stable" => "v^", "v#.#", "devurl" => devurl]
    )


text:
module Testing
using PhaseSpaceIO: IAEAParticle, ParticleType, GeneralizedEGSParticle,
EGSHeader, EGSParticle, photon, electron, positron, compute_u_v_w
using PhaseSpaceIO: Latch
export arbitrary
export testdatapath

function testdatapath(args...)
    joinpath(@__DIR__, "..", "test", "assets", args...)
end

function arbitrary(::Type{IAEAParticle{Nf,Ni}}) where {Nf, Ni}
    typ = rand([instances(ParticleType)...])
    E = 100*rand()
    weight = rand()
    x = randn(Float32)
    y = randn(Float32)
    z = randn(Float32)
    u = randn(Float32)
    v = randn(Float32)
    w = randn(Float32)
    scale = 1f0/sqrt(u^2 + v^2 + w^2)
    u *= scale
    v *= scale
    u,v,w = compute_u_v_w(u, v, sign(w))
    new_history  = rand(Bool)
    extra_floats = tuple(randn(Float32, Nf)...)
    extra_ints   = tuple(rand(Int32, Ni)...)
    IAEAParticle{Nf,Ni}(typ,
                    E,
                    weight,
                    x,y,z,
                    u,v,w,
                    new_history,
                    extra_floats,
                    extra_ints)
end

function arbitrary(H::Type{EGSHeader{P}};
    particlecount = Int32(rand(1:typemax(Int32))),
    photoncount = Int32(rand(1:particlecount)),
    max_E_kin = Float32(25*rand()),
    min_E_kin_electrons= Float32(rand()),
    originalcount = Float32(10^(10*rand()))
    ) where {P}
    H(particlecount,
    photoncount,
    max_E_kin,
    min_E_kin_electrons,
    originalcount)
end

function arbitrary(::Type{Latch})
    Latch(charge=rand(-1:1))
end

arbitrary_z(::Type{Nothing}) = nothing
arbitrary_z(::Type{Float32}) = randn(Float32)

function arbitrary(P::Type{GeneralizedEGSParticle{ZLAST,Z}};
        latch=arbitrary(Latch),
        E = rand(Float32),
        x = randn(Float32),
        y = randn(Float32),
        z = arbitrary_z(Z),
        ) where {ZLAST,Z}

    u = randn()
    v = randn()
    w = randn()
    scale = 1/sqrt(u^2 + v^2 + w^2)
    u = Float32(scale * u)
    v = Float32(scale * v)
    w = Float32(scale * w)
    new_history = rand(Bool)
    weight = 10*rand(Float32)
    if ZLAST == Nothing
        zlast = nothing
    else
        zlast = randn(Float32)
    end
    GeneralizedEGSParticle(
        latch=latch,
        E=E,
        weight=weight,
        x=x,
        y=y,
        z=z,
        u=u,
        v=v,
        w=w,
        new_history=new_history,
        zlast=zlast,
   )
end

end


text:
using Distributed
using Test

(@isdefined verbose) || (verbose = false)

srv_proc = addprocs(1)[]

@everywhere begin
    using Pkg
    Pkg.activate("LibModbus")
    Pkg.instantiate()
    using LibModbus
    include("unit_test_constants.jl")
    include("unit_test_functions.jl")
end

include("unit_test_client.jl")

fetch(@spawnat srv_proc @eval Main mb_mapping=init_server_mapping())
fetch(@spawnat srv_proc @eval Main ctx_server = TcpContext("127.0.0.1", 1502))

remote_do(()->start_unit_test_server(ctx_server, mb_mapping), srv_proc)

unit_test_client(verbose)

rmprocs(srv_proc)


text:
## Sorting

function issorted(x::RLEVector, order::Ordering)
  issorted(x.runvalues)
end

function reverse(x::RLEVector, start=1, stop=length(x))
  rle = RLEVector(reverse(x.runvalues), cumsum(reverse(widths(x))))
  return(rle)
end

function reverse!(x::RLEVector, start=1, stop=length(x))
  reverse!(x.runvalues)
  x.runends[:] = cumsum(reverse(widths(x)))
  return(x)
end

function permute_runs(x::RLEVector, indices)
  RLEVector(rvalue(x)[indices], cumsum(widths(x)[indices]))
end

function sort(x::RLEVector{T1,T2}) where {T1,T2}
  ord = sortperm(x.runvalues)
  rle = RLEVector{T1,T2}( x.runvalues[ord], cumsum(widths(x)[ord]) )
  return(rle)
end

function sort!(x::RLEVector)
  ord = sortperm(x.runvalues)
  x.runvalues[:] = x.runvalues[ord]
  x.runends[:] = cumsum(widths(x)[ord])
  return(x)
end

function sortperm(x::RLEVector)
  ord = sortperm(x.runvalues)
  RLEVector(ord, cumsum(widths(x)[ord]))
end


text:
# Betweenness centrality measures
# TODO - weighted, separate unweighted, edge betweenness

doc"""
    betweenness_centrality(g, k=0; normalize=true, endpoints=false)

Calculates the [betweenness centrality](https://en.wikipedia.org/wiki/Centrality#Betweenness_centrality) of
the graph `g`, or, optionally, of a random subset of `k` vertices. Can
optionally include endpoints in the calculations. Normalization is enabled by
default.

Betweenness centrality is defined as:

$bc(v) = \frac{1}{\mathcal{N}} \sum_{s \neq t \neq v}
        \frac{\sigma_{st}(v)}{\sigma_{st}}$.

 **Parameters**

g: SimpleGraph
    A Graph, directed or undirected.

k: Integer, optional
    Use `k` nodes sample to estimate the betweenness centrality. If none,
    betweenness centrality is computed using the `n` nodes in the graph.

normalize: bool, optional
    If true, the betweenness values are normalized by the total number
    of possible distinct paths between all pairs in the graphs. For an undirected graph,
    this number if `((n-1)*(n-2))/2` and for a directed graph, `(n-1)*(n-2)`
    where `n` is the number of nodes in the graph.

endpoints: bool, optional
    If true, endpoints are included in the shortest path count.

**Returns**

betweenness: Array{Float64}
    Betweenness centrality value per node id.

**References**

[1] Brandes 2001 & Brandes 2008
"""
function betweenness_centrality(
    g::SimpleGraph,
    k::Integer=0;
    normalize=true,
    endpoints=false)

    n_v = nv(g)
    isdir = is_directed(g)

    betweenness = zeros(n_v)
    if k == 0
        nodes = 1:n_v
    else
        nodes = sample!([1:n_v;], k)   #112
    end
    for s in nodes
        if degree(g,s) > 0  # this might be 1?
            state = dijkstra_shortest_paths(g, s; allpaths=true)
            if endpoints
                _accumulate_endpoints!(betweenness, state, g, s)
            else
                _accumulate_basic!(betweenness, state, g, s)
            end
        end
    end

    _rescale!(betweenness,
              n_v,
              normalize,
              isdir,
              k)

    return betweenness
end

function _accumulate_basic!(
    betweenness::Vector{Float64},
    state::DijkstraState,
    g::SimpleGraph,
    si::Integer
    )

    n_v = length(state.parents) # this is the ttl number of vertices
    δ = zeros(n_v)
    σ = state.pathcounts
    P = state.predecessors

    # make sure the source index has no parents.
    P[si] = []
    # we need to order the source nodes by decreasing distance for this to work.
    S = sortperm(state.dists, rev=true)
    for w in S
        coeff = (1.0 + δ[w]) / σ[w]
        for v in P[w]
            if v > 0
                δ[v] += (σ[v] * coeff)
            end
        end
        if w != si
            betweenness[w] += δ[w]
        end
    end
end

function _accumulate_endpoints!(
    betweenness::Vector{Float64},
    state::DijkstraState,
    g::SimpleGraph,
    si::Integer
    )

    n_v = nv(g) # this is the ttl number of vertices
    δ = zeros(n_v)
    σ = state.pathcounts
    P = state.predecessors
    v1 = [1:n_v;]
    v2 = state.dists
    S = sortperm(state.dists, rev=true)
    s = g.vertices[si]
    betweenness[s] += length(S) - 1    # 289

    for w in S
        coeff = (1.0 + δ[w]) / σ[w]
        for v in P[w]
            δ[v] += σ[v] * coeff
        end
        if w != si
            betweenness[w] += (δ[w] + 1)
        end
    end
end

function _rescale!(betweenness::Vector{Float64}, n::Int, normalize::Bool, directed::Bool, k::Int)
    if normalize
        if n <= 2
            do_scale = false
        else
            do_scale = true
            scale = 1.0 / ((n - 1) * (n - 2))
        end
    else
        if !directed
            do_scale = true
            scale = 1.0 / 2.0
        else
            do_scale = false
        end
    end
    if do_scale
        if k > 0
            scale = scale * n / k
        end
        for v = 1:length(betweenness)
            betweenness[v] *= scale
        end
    end
end