text:
# Parts of this code were taken / derived from NetworkX. See LICENSE for
# licensing details.

"""
    complete_graph(n)

Create an undirected [complete graph](https://en.wikipedia.org/wiki/Complete_graph)
with `n` vertices.

# Examples
```jldoctest
julia> complete_graph(5)
{5, 10} undirected simple Int64 graph

julia> complete_graph(Int8(6))
{6, 15} undirected simple Int8 graph
```
"""
function complete_graph(n::T) where {T <: Integer}
    n <= 0 && return SimpleGraph{T}(0)
    ne = Int(n * (n - 1) ÷ 2)
    fadjlist = Vector{Vector{T}}(undef, n)
    @inbounds @simd for u = 1:n
        listu = Vector{T}(undef, n-1)
        listu[1:(u-1)] = 1:(u - 1)
        listu[u:(n-1)] = (u + 1):n
        fadjlist[u] = listu
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    complete_bipartite_graph(n1, n2)

Create an undirected [complete bipartite graph](https://en.wikipedia.org/wiki/Complete_bipartite_graph)
with `n1 + n2` vertices.

# Examples
```jldoctest
julia> complete_bipartite_graph(3, 4)
{7, 12} undirected simple Int64 graph

julia> complete_bipartite_graph(Int8(3), Int8(4))
{7, 12} undirected simple Int8 graph
```
"""
function complete_bipartite_graph(n1::T, n2::T) where {T <: Integer}
    (n1 < 0 || n2 < 0) && return SimpleGraph{T}(0)
    Tw = widen(T)
    nw = Tw(n1) + Tw(n2)
    n = T(nw)  # checks if T is large enough for n1 + n2

    ne = Int(n1) * Int(n2)

    fadjlist = Vector{Vector{T}}(undef, n)
    range1 = 1:n1
    range2 = (n1 + 1):n
    @inbounds @simd for u in range1
        fadjlist[u] = Vector{T}(range2)
    end
    @inbounds @simd for u in range2
        fadjlist[u] = Vector{T}(range1)
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    complete_multipartite_graph(partitions)

Create an undirected [complete bipartite graph](https://en.wikipedia.org/wiki/Complete_bipartite_graph)
with `sum(partitions)` vertices. A partition with `0` vertices is skipped.

### Implementation Notes
Preserves the eltype of the partitions vector. Will error if the required number of vertices
exceeds the eltype.

# Examples
```jldoctest
julia> complete_multipartite_graph([1,2,3])
{6, 11} undirected simple Int64 graph

julia> complete_multipartite_graph(Int8[5,5,5])
{15, 75} undirected simple Int8 graph
```
"""
function complete_multipartite_graph(partitions::AbstractVector{T}) where {T <: Integer}
    any(x -> x < 0, partitions) && return SimpleGraph{T}(0)
    length(partitions) == 1 && return SimpleGraph{T}(partitions[1])
    length(partitions) == 2 && return complete_bipartite_graph(partitions[1], partitions[2])

    n = sum(partitions)

    ne = 0
    for p in partitions # type stability fails if we use sum and a generator here
        ne += p*(Int(n)-p) # overflow if we don't convert to Int
    end
    ne = div(ne, 2)

    fadjlist = Vector{Vector{T}}(undef, n)
    cur = 1
    for p in partitions
        currange = cur:(cur+p-1) # all vertices in the current partition
        lowerrange = 1:(cur-1)   # all vertices lower than the current partition
        upperrange = (cur+p):n   # all vertices higher than the current partition
        @inbounds @simd for u in currange
            fadjlist[u] = Vector{T}(undef, length(lowerrange) + length(upperrange))
            fadjlist[u][1:length(lowerrange)] = lowerrange
            fadjlist[u][(length(lowerrange)+1):end] = upperrange
        end
        cur += p
    end

    return SimpleGraph{T}(ne, fadjlist)
end

"""
    turan_graph(n, r)

Creates a [Turán Graph](https://en.wikipedia.org/wiki/Tur%C3%A1n_graph), a complete
multipartite graph with `n` vertices and `r` partitions.

# Examples
```jldoctest
julia> turan_graph(6, 2)
{6, 9} undirected simple Int64 graph

julia> turan_graph(Int8(7), 2)
{7, 12} undirected simple Int8 graph
```
"""
function turan_graph(n::Integer, r::Integer)
    !(1 <= r <= n) && throw(DomainError("n=$n and r=$r are invalid, must satisfy 1 <= r <= n"))
    T = typeof(n)
    partitions = Vector{T}(undef, r)
    c = cld(n,r)
    f = fld(n,r)
    @inbounds @simd for i in 1:(n%r)
        partitions[i] = c
    end
    @inbounds @simd for i in ((n%r)+1):r
        partitions[i] = f
    end
    return complete_multipartite_graph(partitions)
end

"""
    complete_digraph(n)

Create a directed [complete graph](https://en.wikipedia.org/wiki/Complete_graph)
with `n` vertices.

# Examples
```jldoctest
julia> complete_digraph(5)
{5, 20} directed simple Int64 graph

julia> complete_digraph(Int8(6))
{6, 30} directed simple Int8 graph
```
"""
function complete_digraph(n::T) where {T <: Integer}
    n <= 0 && return SimpleDiGraph{T}(0)

    ne = Int(n * (n - 1))
    fadjlist = Vector{Vector{T}}(undef, n)
    badjlist = Vector{Vector{T}}(undef, n)
    @inbounds @simd for u = 1:n
        listu = Vector{T}(undef, n-1)
        listu[1:(u-1)] = 1:(u - 1)
        listu[u:(n-1)] = (u + 1):n
        fadjlist[u] = listu
        badjlist[u] = deepcopy(listu)
    end
    return SimpleDiGraph(ne, fadjlist, badjlist)
end

"""
    star_graph(n)

Create an undirected [star graph](https://en.wikipedia.org/wiki/Star_(graph_theory))
with `n` vertices.

# Examples
```jldoctest
julia> star_graph(3)
{3, 2} undirected simple Int64 graph

julia> star_graph(Int8(10))
{10, 9} undirected simple Int8 graph
```
"""
function star_graph(n::T) where {T <: Integer}
    n <= 0 && return SimpleGraph{T}(0)

    ne = Int(n - 1)
    fadjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = Vector{T}(2:n)
    @inbounds @simd for u = 2:n
        fadjlist[u] = T[1]
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    star_digraph(n)

Create a directed [star graph](https://en.wikipedia.org/wiki/Star_(graph_theory))
with `n` vertices.

# Examples
```jldoctest
julia> star_digraph(3)
{3, 2} directed simple Int64 graph

julia> star_digraph(Int8(10))
{10, 9} directed simple Int8 graph
```
"""
function star_digraph(n::T) where {T <: Integer}
    n <= 0 && return SimpleDiGraph{T}(0)

    ne = Int(n - 1)
    fadjlist = Vector{Vector{T}}(undef, n)
    badjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = Vector{T}(2:n)
    @inbounds badjlist[1] = T[]
    @inbounds @simd for u = 2:n
        fadjlist[u] = T[]
        badjlist[u] = T[1]
    end
    return SimpleDiGraph(ne, fadjlist, badjlist)
end

"""
    path_graph(n)

Create an undirected [path graph](https://en.wikipedia.org/wiki/Path_graph)
with `n` vertices.

# Examples
```jldoctest
julia> path_graph(5)
{5, 4} undirected simple Int64 graph

julia> path_graph(Int8(10))
{10, 9} undirected simple Int8 graph
```
"""
function path_graph(n::T) where {T <: Integer}
    n <= 1 && return SimpleGraph(n)

    ne = Int(n - 1)
    fadjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = T[2]
    @inbounds fadjlist[n] = T[n - 1]

    @inbounds @simd for u = 2:(n-1)
        fadjlist[u] = T[u - 1, u + 1]
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    path_digraph(n)

Creates a directed [path graph](https://en.wikipedia.org/wiki/Path_graph)
with `n` vertices.

# Examples
```jldoctest
julia> path_digraph(5)
{5, 4} directed simple Int64 graph

julia> path_digraph(Int8(10))
{10, 9} directed simple Int8 graph
```
"""
function path_digraph(n::T) where {T <: Integer}
    n <= 1 && return SimpleDiGraph(n)

    ne = Int(n - 1)
    fadjlist = Vector{Vector{T}}(undef, n)
    badjlist = Vector{Vector{T}}(undef, n)

    @inbounds fadjlist[1] = T[2]
    @inbounds badjlist[1] = T[]
    @inbounds fadjlist[n] = T[]
    @inbounds badjlist[n] = T[n - 1]

    @inbounds @simd for u = 2:(n-1)
        fadjlist[u] = T[u + 1]
        badjlist[u] = T[u - 1]
    end
    return SimpleDiGraph(ne, fadjlist, badjlist)
end

"""
    cycle_graph(n)

Create an undirected [cycle graph](https://en.wikipedia.org/wiki/Cycle_graph)
with `n` vertices.

# Examples
```jldoctest
julia> cycle_graph(3)
{3, 3} undirected simple Int64 graph

julia> cycle_graph(Int8(5))
{5, 5} undirected simple Int8 graph
```
"""
function cycle_graph(n::T) where {T <: Integer}
    n <= 1 && return SimpleGraph(n)
    n == 2 && return SimpleGraph(Edge{T}.([(1, 2)]))

    ne = Int(n)
    fadjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = T[2, n]
    @inbounds fadjlist[n] = T[1, n-1]

    @inbounds @simd for u = 2:(n-1)
        fadjlist[u] = T[u-1, u+1]
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    cycle_digraph(n)

Create a directed [cycle graph](https://en.wikipedia.org/wiki/Cycle_graph)
with `n` vertices.

# Examples
```jldoctest
julia> cycle_digraph(3)
{3, 3} directed simple Int64 graph

julia> cycle_digraph(Int8(5))
{5, 5} directed simple Int8 graph
```
"""
function cycle_digraph(n::T) where {T <: Integer}
    n <= 1 && return SimpleDiGraph(n)
    n == 2 && return SimpleDiGraph(Edge{T}.([(1, 2), (2, 1)]))

    ne = Int(n)
    fadjlist = Vector{Vector{T}}(undef, n)
    badjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = T[2]
    @inbounds badjlist[1] = T[n]
    @inbounds fadjlist[n] = T[1]
    @inbounds badjlist[n] = T[n-1]

    @inbounds @simd for u = 2:(n-1)
        fadjlist[u] = T[u + 1]
        badjlist[u] = T[u + -1]
    end
    return SimpleDiGraph(ne, fadjlist, badjlist)
end

"""
    wheel_graph(n)

Create an undirected [wheel graph](https://en.wikipedia.org/wiki/Wheel_graph)
with `n` vertices.

# Examples
```jldoctest
julia> wheel_graph(5)
{5, 8} undirected simple Int64 graph

julia> wheel_graph(Int8(6))
{6, 10} undirected simple Int8 graph
```
"""
function wheel_graph(n::T) where {T <: Integer}
    n <= 1 && return SimpleGraph(n)
    n <= 3 && return cycle_graph(n)

    ne = Int(2 * (n - 1))
    fadjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = Vector{T}(2:n)
    @inbounds fadjlist[2] = T[1, 3, n]
    @inbounds fadjlist[n] = T[1, 2, n - 1]

    @inbounds @simd for u = 3:(n-1)
        fadjlist[u] = T[1, u - 1, u + 1]
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    wheel_digraph(n)

Create a directed [wheel graph](https://en.wikipedia.org/wiki/Wheel_graph)
with `n` vertices.

# Examples
```jldoctest
julia> wheel_digraph(5)
{5, 8} directed simple Int64 graph

julia> wheel_digraph(Int8(6))
{6, 10} directed simple Int8 graph
```
"""
function wheel_digraph(n::T) where {T <: Integer}
    n <= 2 && return path_digraph(n)
    n == 3 && return SimpleDiGraph(Edge{T}.([(1,2),(1,3),(2,3),(3,2)]))

    ne = Int(2 * (n - 1))
    fadjlist = Vector{Vector{T}}(undef, n)
    badjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = Vector{T}(2:n)
    @inbounds badjlist[1] = T[]
    @inbounds fadjlist[2] = T[3]
    @inbounds badjlist[2] = T[1, n]
    @inbounds fadjlist[n] = T[2]
    @inbounds badjlist[n] = T[1, n - 1]

    @inbounds @simd for u = 3:(n-1)
        fadjlist[u] = T[u + 1]
        badjlist[u] = T[1, u - 1]
    end
    return SimpleDiGraph(ne, fadjlist, badjlist)
end

"""
    grid(dims; periodic=false)

Create a ``|dims|``-dimensional cubic lattice, with length `dims[i]`
in dimension `i`.

### Optional Arguments
- `periodic=false`: If true, the resulting lattice will have periodic boundary
condition in each dimension.

# Examples
```jldoctest
julia> grid([2,3])
{6, 7} undirected simple Int64 graph

julia> grid(Int8[2, 2, 2], periodic=true)
{8, 12} undirected simple Int8 graph

julia> grid((2,3))
{6, 7} undirected simple Int64 graph
```
"""
function grid(dims::AbstractVector{T}; periodic=false) where {T <: Integer}
    # checks if T is large enough for product(dims)
    Tw = widen(T)
    n = one(T)
    for d in dims
        d <= 0 && return SimpleGraph{T}(0)
        nw = Tw(n) * Tw(d)
        n = T(nw)
    end

    if periodic
        g = cycle_graph(dims[1])
        for d in dims[2:end]
            g = cartesian_product(cycle_graph(d), g)
        end
    else
        g = path_graph(dims[1])
        for d in dims[2:end]
            g = cartesian_product(path_graph(d), g)
        end
    end
    return g
end
grid(dims::Tuple; periodic=false) = grid(collect(dims); periodic=periodic)

"""
    binary_tree(k::Integer)

Create a [binary tree](https://en.wikipedia.org/wiki/Binary_tree)
of depth `k`.

# Examples
```jldoctest
julia> binary_tree(4)
{15, 14} undirected simple Int64 graph

julia> binary_tree(Int8(5))
{31, 30} undirected simple Int8 graph
```
"""
function binary_tree(k::T) where {T <: Integer}
    k <= 0 && return SimpleGraph(0)
    k == 1 && return SimpleGraph(1)
    if LightGraphs.isbounded(k) && BigInt(2) ^ k - 1 > typemax(k)
        throw(DomainError(k, "2^k - 1 not representable by type $T"))
    end
    n = T(2^k - 1)

    ne = Int(n - 1)
    fadjlist = Vector{Vector{T}}(undef, n)
    @inbounds fadjlist[1] = T[2, 3]
    @inbounds for i in 1:(k - 2)
        @simd for j in (2^i):(2^(i + 1) - 1)
            fadjlist[j] = T[j ÷ 2, 2j, 2j + 1]
        end
    end
    i = k - 1
    @inbounds @simd for j in (2^i):(2^(i + 1) - 1)
        fadjlist[j] = T[j ÷ 2]
    end
    return SimpleGraph(ne, fadjlist)
end

"""
    double_binary_tree(k::Integer)

Create a double complete binary tree with `k` levels.

### References
- Used as an example for spectral clustering by Guattery and Miller 1998.

# Examples
```jldoctest
julia> double_binary_tree(4)
{30, 29} undirected simple Int64 graph

julia> double_binary_tree(Int8(5))
{62, 61} undirected simple Int8 graph
```
"""
function double_binary_tree(k::Integer)
    gl = binary_tree(k)
    gr = binary_tree(k)
    g = blockdiag(gl, gr)
    add_edge!(g, 1, nv(gl) + 1)
    return g
end

"""
    roach_graph(k)

Create a Roach graph of size `k`.

### References
- Guattery and Miller 1998

# Examples
```jldoctest
julia> roach_graph(10)
{40, 48} undirected simple Int64 graph
```
"""
function roach_graph(k::Integer)
    dipole = complete_graph(2)
    nopole = SimpleGraph(2)
    antannae = crosspath(k, nopole)
    body = crosspath(k, dipole)
    roach = blockdiag(antannae, body)
    add_edge!(roach, nv(antannae) - 1, nv(antannae) + 1)
    add_edge!(roach, nv(antannae), nv(antannae) + 2)
    return roach
end

"""
    clique_graph(k, n)

Create a graph consisting of `n` connected `k`-cliques.

# Examples
```jldoctest
julia> clique_graph(4, 10)
{40, 70} undirected simple Int64 graph

julia> clique_graph(Int8(10), Int8(4))
{40, 184} undirected simple Int8 graph
```
"""
function clique_graph(k::T, n::T) where {T <: Integer}
    Tw = widen(T)
    knw = Tw(k) * Tw(n)
    kn = T(knw)  # checks if T is large enough for k * n

    g = SimpleGraph(kn)
    for c = 1:n
        for i = ((c - 1) * k + 1):(c * k - 1), j = (i + 1):(c * k)
            add_edge!(g, i, j)
        end
    end
    for i = 1:(n - 1)
        add_edge!(g, (i - 1) * k + 1, i * k + 1)
    end
    add_edge!(g, 1, (n - 1) * k + 1)
    return g
end

"""
    ladder_graph(n)

Create a [ladder graph](https://en.wikipedia.org/wiki/Ladder_graph) consisting of `2n` nodes and `3n-2` edges.

### Implementation Notes
Preserves the eltype of `n`. Will error if the required number of vertices
exceeds the eltype.

# Examples
```jldoctest
julia> ladder_graph(3)
{6, 7} undirected simple Int64 graph

julia> ladder_graph(Int8(4))
{8, 10} undirected simple Int8 graph
```
"""
function ladder_graph(n::T) where {T <: Integer}
    n <= 0 && return SimpleGraph{T}(0)
    n == 1 && return path_graph(T(2))
    Tw = widen(T)
    temp = T(Tw(n)+Tw(n)) # test to check if T is large enough

    fadjlist = Vector{Vector{T}}(undef, 2*n)
    @inbounds @simd for i in 2:(n-1)
        fadjlist[i]   = T[i-1, i+1, i+n]
        fadjlist[n+i] = T[i, n+i-1, n+i+1]
    end
    fadjlist[1]   = T[2, n+1]
    fadjlist[n+1] = T[1, n+2]
    fadjlist[n]   = T[n-1, 2*n]
    fadjlist[2*n] = T[n, 2*n-1]

    return SimpleGraph(3*n-2, fadjlist)
end

"""
    circular_ladder_graph(n)

Create a [circular ladder graph](https://en.wikipedia.org/wiki/Ladder_graph#Circular_ladder_graph) consisting of `2n` nodes and `3n` edges.
This is also known as the [prism graph](https://en.wikipedia.org/wiki/Prism_graph).

### Implementation Notes
Preserves the eltype of the partitions vector. Will error if the required number of vertices
exceeds the eltype.
`n` must be at least 3 to avoid self-loops and multi-edges.

# Examples
```jldoctest
julia> circular_ladder_graph(3)
{6, 9} undirected simple Int64 graph

julia> circular_ladder_graph(Int8(4))
{8, 12} undirected simple Int8 graph
```
"""
function circular_ladder_graph(n::Integer)
    n < 3 && throw(DomainError("n=$n must be at least 3"))
    g = ladder_graph(n)
    add_edge!(g, 1, n)
    add_edge!(g, n+1, 2*n)
    return g
end

"""
    barbell_graph(n1, n2)

Create a [barbell graph](https://en.wikipedia.org/wiki/Barbell_graph) consisting of a clique of size `n1` connected by an edge to a clique of size `n2`.

### Implementation Notes
Preserves the eltype of `n1` and `n2`. Will error if the required number of vertices
exceeds the eltype.
`n1` and `n2` must be at least 1 so that both cliques are non-empty.
The cliques are organized with nodes `1:n1` being the left clique and `n1+1:n1+n2` being the right clique. The cliques are connected by and edge `(n1, n1+1)`.

# Examples
```jldoctest
julia> barbell_graph(3, 4)
{7, 10} undirected simple Int64 graph

julia> barbell_graph(Int8(5), Int8(5))
{10, 21} undirected simple Int8 graph
```
"""
function barbell_graph(n1::T, n2::T) where {T <: Integer}
    (n1 < 1 || n2 < 1) && throw(DomainError("n1=$n1 and n2=$n2 must be at least 1"))

    n = Base.Checked.checked_add(n1, n2) # check for overflow
    fadjlist = Vector{Vector{T}}(undef, n)

    ne = Int(n1)*(n1-1)÷2 + Int(n2)*(n2-1)÷2

    @inbounds @simd for u = 1:n1
        listu = Vector{T}(undef, n1-1)
        listu[1:(u-1)] = 1:(u-1)
        listu[u:(n1-1)] = (u+1):n1
        fadjlist[u] = listu
    end

    @inbounds for u = 1:n2
        listu = Vector{T}(undef, n2-1)
        listu[1:(u-1)] = (n1+1):(n1+(u-1))
        listu[u:(n2-1)] = (n1+u+1):(n1+n2)
        fadjlist[n1+u] = listu
    end

    g = SimpleGraph(ne, fadjlist)
    add_edge!(g, n1, n1+1)
    return g
end

"""
    lollipop_graph(n1, n2)

Create a [lollipop graph](https://en.wikipedia.org/wiki/Lollipop_graph) consisting of a clique of size `n1` connected by an edge to a path of size `n2`.

### Implementation Notes
Preserves the eltype of `n1` and `n2`. Will error if the required number of vertices
exceeds the eltype.
`n1` and `n2` must be at least 1 so that both the clique and the path have at least one vertex.
The graph is organized with nodes `1:n1` being the clique and `n1+1:n1+n2` being the path. The clique is connected to the path by an edge `(n1, n1+1)`.

# Examples
```jldoctest
julia> lollipop_graph(2, 5)
{7, 6} undirected simple Int64 graph

julia> lollipop_graph(Int8(3), Int8(4))
{7, 7} undirected simple Int8 graph
```
"""
function lollipop_graph(n1::T, n2::T) where {T <: Integer}
    (n1 < 1 || n2 < 1) && throw(DomainError("n1=$n1 and n2=$n2 must be at least 1"))

    if n1 == 1
        return path_graph(T(n2+1))
    elseif n1 > 1 && n2 == 1
        g = complete_graph(n1)
        add_vertex!(g)
        add_edge!(g, n1, n1+1)
        return g
    end

    n = Base.Checked.checked_add(n1, n2) # check for overflow
    fadjlist = Vector{Vector{T}}(undef, n)

    ne = Int(Int(n1)*(n1-1)÷2 + n2-1)

    @inbounds @simd for u = 1:n1
        listu = Vector{T}(undef, n1-1)
        listu[1:(u-1)] = 1:(u-1)
        listu[u:(n1-1)] = (u+1):n1
        fadjlist[u] = listu
    end

    @inbounds fadjlist[n1+1] = T[n1+2]
    @inbounds fadjlist[n1+n2] = T[n1+n2-1]

    @inbounds @simd for u = (n1+2):(n1+n2-1)
        fadjlist[u] = T[u-1, u+1]
    end

    g = SimpleGraph(ne, fadjlist)
    add_edge!(g, n1, n1+1)
    return g
end


text:
@testset "Capillary pressure" begin
    # Struct
    sys = Swalbe.SysConst(Lx=5, Ly=5, n=3, m=2, γ=1.0, hmin=0.1, hcrit=0.1)
    state = Swalbe.Sys(sys, "CPU")
    state.height .= reshape(collect(1.0:25),5,5)
    # With arguments struct
    f = reshape(collect(1.0:25),5,5)
    f_float = reshape(collect(1.0f0:25.0f0),5,5)
    res = zeros(5,5)
    @testset "No contact angle" begin
        Swalbe.filmpressure!(res, f, 1.0, 0.0, 3, 2, 0.1, 0.1)
        Swalbe.filmpressure!(state, sys, 0.0)
        sol = [-30.0 -5.0 -5.0 -5.0 20;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -20.0 5.0 5.0 5.0 30.0;
               ]

        @test all(isapprox.(res, sol; atol=1e-10))
        @test all(isapprox.(state.pressure, sol; atol=1e-10))

        Swalbe.filmpressure!(res, f, 0.0)
        @test all(isapprox.(res, 0.01 .* sol; atol=1e-10))
    end
    @testset "No height gradient" begin
        nograd = ones(5,5)
        state.height .= 1.0
        sys2 = Swalbe.SysConst(Lx=5, Ly=5, n=3, m=2, γ=1.0, hmin=0.1, hcrit=0.0, θ=1/2)
        Swalbe.filmpressure!(res, nograd, 1.0, 1/2, 3, 2, 0.1, 0.0)
        Swalbe.filmpressure!(state, sys2, 1/2)
        @test all(isapprox.(res, -2(0.1^2-0.1); atol=1e-10))
        @test all(isapprox.(state.pressure, -2(0.1^2-0.1); atol=1e-10))
        Swalbe.filmpressure!(state, sys2)
        @test all(isapprox.(state.pressure, -2(0.1^2-0.1); atol=1e-10))
    end
    @testset "Gradient and contact angle" begin
        state.height .= reshape(collect(1.0:25),5,5)
        sys2 = Swalbe.SysConst(Lx=5, Ly=5, n=3, m=2, γ=1.0, hmin=0.1, hcrit=0.0)
        Swalbe.filmpressure!(state, sys2, 1/2)
        Swalbe.filmpressure!(res, f, 1.0, 1/2, 3, 2, 0.1, 0.0)
        sol = [30.0 5.0 5.0 5.0 -20;
               25.0 0.0 0.0 0.0 -25.0;
               25.0 0.0 0.0 0.0 -25.0;
               25.0 0.0 0.0 0.0 -25.0;
               20.0 -5.0 -5.0 -5.0 -30.0]

        @test all(isapprox.(res, -1 .* (sol .+ 20 .* ((0.1 ./ f).^3 .- (0.1 ./ f).^2)); atol=1e-10))
        @test all(isapprox.(state.pressure, -1 .* (sol .+ 20 .* ((0.1 ./ f).^3 .- (0.1 ./ f).^2)); atol=1e-10))

    end
    @testset "Gradient and contact angle Float32" begin
        Swalbe.filmpressure!(res, f_float, 1.0f0, 0.5f0, 3, 2, 0.1f0, 0.0f0)
        sol = [30.0f0 5.0f0 5.0f0 5.0f0 -20f0;
               25.0f0 0.0f0 0.0f0 0.0f0 -25.0f0;
               25.0f0 0.0f0 0.0f0 0.0f0 -25.0f0;
               25.0f0 0.0f0 0.0f0 0.0f0 -25.0f0;
               20.0f0 -5.0f0 -5.0f0 -5.0f0 -30.0f0]

        @test all(isapprox.(res, -1 .* (sol .+ 20 .* ((0.1f0 ./ f_float).^3 .- (0.1f0 ./ f_float).^2)); atol=1e-6))
    end
    dgrad = zeros(5,5,8)
    @testset "No contact angle circshift!" begin
        state.height .= reshape(collect(1.0:25),5,5)
        Swalbe.filmpressure!(state, sys, 0.0)
        Swalbe.filmpressure!(res, f, dgrad, 1.0, 0.0, 3, 2, 0.1, 0.1)
        sol = [-30.0 -5.0 -5.0 -5.0 20;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -25.0 0.0 0.0 0.0 25.0;
               -20.0 5.0 5.0 5.0 30.0;
               ]

        @test all(isapprox.(res, sol; atol=1e-10))
        @test all(isapprox.(state.pressure, sol; atol=1e-10))
    end
end

@testset "Capillary pressure 1D" begin
    # Struct
    sys = Swalbe.SysConst_1D(L=30, n=3, m=2, γ=1.0, hmin=0.1, hcrit=0.1)
    state = Swalbe.Sys(sys)
    state.height .= collect(1.0:30)
    sys2 = Swalbe.SysConst_1D(L=30, n=3, m=2, γ=1.0, hmin=0.1, hcrit=0.0, θ=1/2)
    state2 = Swalbe.Sys(sys2)
    # Without the struct
    f = collect(1.0:30)
    f_float = collect(1.0f0:30.0f0)
    res = zeros(30)
    dummy = zeros(30,3)
    @testset "No contact angle" begin
        Swalbe.filmpressure!(res, f, dummy, 1.0, 0.0, 3, 2, 0.1, 0.1)
        Swalbe.filmpressure!(state, sys, 0.0)
        # println("My result: $res")
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        @test all(res .== -sol)
        @test all(state.pressure .== -sol)
    end
    @testset "No height gradient" begin
        nograd = ones(30)
        state2.height .= 1.0
        Swalbe.filmpressure!(res, nograd, dummy, 1.0, 1/2, 3, 2, 0.1, 0.0)
        Swalbe.filmpressure!(state2, sys2, 1/2)
        @test all(isapprox.(res, -2(0.1^2-0.1); atol=1e-10))
        @test all(isapprox.(state2.pressure, -2(0.1^2-0.1); atol=1e-10))

    end
    @testset "Gradient and contact angle" begin
        Swalbe.filmpressure!(res, f, dummy, 1.0, 1/2, 3, 2, 0.1, 0.0)
        Swalbe.filmpressure!(state, sys2, 1/2)
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        @test all(isapprox.(res, -1 .* (sol .+ 20 .* ((0.1 ./ f).^3 .- (0.1 ./ f).^2)); atol=1e-10))
        @test all(isapprox.(state.pressure, -1 .* (sol .+ 20 .* ((0.1 ./ f).^3 .- (0.1 ./ f).^2)); atol=1e-10))
        # TODO: Fix me
        # state2.height .= sol
        # Swalbe.filmpressure!(state2, sys2)
        # @test all(isapprox.(state2.pressure, -1 .* (sol .+ 20 .* ((0.1 ./ f).^3 .- (0.1 ./ f).^2)); atol=1e-10))

    end
end

@testset "Capillary pressure 1D with Gamma" begin
    f = collect(1.0:30)
    res = zeros(30)
    rho = fill(0.1, 30)
    dummy = zeros(30,2)
    @testset "No Gamma no theta" begin
        Swalbe.filmpressure!(res, f, dummy, rho, 1.0, 0.0, 3, 2, 0.1, 0.1)
        # println("My result: $res")
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        @test all(res .== -sol)
    end
    @testset "No Gamma no laplacian" begin
        nograd = ones(30)
        Swalbe.filmpressure!(res, nograd, dummy, rho, 1.0, 1/2, 3, 2, 0.1, 0.0)
        @test all(isapprox.(res, -2(0.1^2-0.1); atol=1e-10))
    end
    @testset "Gamma but rho=0 no laplacian" begin
        nograd = ones(30)
        Swalbe.filmpressure!(res, nograd, dummy, zeros(30), 1.0, 1/2, 3, 2, 0.1, 0.0)
        @test all(isapprox.(res, -2(0.1^2-0.1); atol=1e-10))
    end
    @testset "Gamma, rho, gradient and contact angle" begin
        rho = fill(0.1, 30)
        Swalbe.filmpressure!(res, f, dummy, rho, 1.0, 1/2, 3, 2, 0.1, 0.0, Gamma=0.1)
        sol = zeros(30)
        sol[1] = 30
        sol[end] = -30
        @test all(isapprox.(res, -(1 + 0.01) .* (sol .+ 20 .* ((0.1 ./ f).^3 .- (0.1 ./ f).^2)); atol=1e-10))
    end
end

@testset "Power function" begin
    arg = 2
    argsimple = 5.0
    argsimplefp32 = 5.0f0
    argvec = [2.0 3.0 4.0]
    n = [3 6]
    r = [8 64]
    r1 = [125.0 15625.0]
    r2 = [125.0f0 15625.0f0]
    r3 = [8.0 27.0 64.0; 64.0 729.0 4096.0]
    for i in 1:2
        res0 = Swalbe.power_broad(arg, n[i])
        @test res0 == r[i]
        res0 = Swalbe.power_broad(argsimple, n[i])
        @test res0 == r1[i]
        res1 = Swalbe.power_broad(argsimplefp32, n[i])
        @test res1 == r2[i]
        res2 = Swalbe.power_broad.(argvec, n[i])
        for j in 1:3
            @test res2[j] .== r3[i,j]
        end
    end
end

ext: jl


text:
function fn_test()
    @test isa(filesep(), Char);

    @test paths_equal("abc", "abc/")
    @test !paths_equal("abc", "abcd")

    fPath = "/abd/def.txt";
    @test change_extension(fPath, "doc") == "/abd/def.doc"
    @test change_extension(fPath, ".doc") == "/abd/def.doc"
    @test change_extension("/abc/def", "doc") == "/abc/def.doc"

    @test add_extension(fPath, ".new") == fPath
    @test add_extension("abc/def", "txt") == "abc/def.txt"

    @test extensions_equal("txt", ".txt")
    @test !extensions_equal(".txt", ".txt2")

    @test has_extension("abc/def.txt", ".txt")
    @test has_extension("abc/def.txt", "txt")
    @test !has_extension("abc/def.txt", "txt2")
end

@testset "File names" begin
    fn_test()
end

# --------------


text:
using Plots,LaTeXStrings
pgfplotsx()
default(legend=nothing,aspect_ratio=1)

plt = plot(showaxis=false,grid=false,xaxis=:none,yaxis=false)

plot!([0,.4, .4, 0],[.9, .9, .1, .1],l=nothing,seriestype=:shape,fillcolor=RGB(1,1,.7))
plot!([1, .6, .6, 1],[.9, .9, .1, .1],l=nothing,seriestype=:shape,fillcolor=RGB(1,1,.7))

plot!([.86, .7],[.75, .2],l=(:dash,2,:black))
plot!([.2, .3],[.75, .5],l=(:dash,2,:black))

plot!([.2, .86],[.75,.75],l=(2,:darkblue),m=5,markercolor=:black,arrow=true)
plot!([.2, .7],[.75,.2],l=(2,:darkred),arrow=true)
plot!([.3, .7],[.5,.2],l=(2,:darkblue),m=5,markercolor=:black,arrow=true)

annotate!(.875,.75,text(L"y",:left))
annotate!(.72,.18,text(L"\tilde{y}",:left))
annotate!(.18,.77,text(L"x",:right))
annotate!(.29,.47,text(L"\tilde{x}",:right))
annotate!(.5,.79,text(L"f(x)",:top))
annotate!(.51,.52,text(L"\tilde{f}(x)",:bottom))
annotate!(.48,.25,text(L"f(\tilde{x})",:top))
annotate!(.79,.475,text("error",:left))
annotate!(.26,.6,text("backward error",:right))
annotate!(.2,.95,text("data",:center))
annotate!(.8,.95,text("results",:center))

xlims!(0,1)
ylims!(0,1)

savefig(plt,"backwarderror.svg")

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
using Lasso
DATADIR = joinpath(dirname(@__FILE__), "data")

# Test against fused lasso results from glmgen
lakehuron = readcsvmat(joinpath(DATADIR, "LakeHuron.csv");header=true)[:, 3]
@testset "FusedLasso" begin
	@testset "λ = $(lambda)" for lambda in (10, 1, 0.1)
		@test coef(fit(FusedLasso, lakehuron, lambda)) ≈ vec(readcsvmat(joinpath(DATADIR, "LakeHuron_lambda_$lambda.csv")))
	end
end


text:
using Random: rand, RandomDevice
using Primes: nextprime
"""
    symmetric_key(bytes::Int64)

Generates a symmetric key with the number of bytes given by `bytes`.
"""
symmetric_key(bytes::Int64) = rand(RandomDevice(), UInt8, bytes)

function asymmetric_key(bits=20, e=2^16+1)
    lower = big(1) << (bits - 1)
    upper = (big(1) << bits) - 1

    p = nextprime(rand(RandomDevice(), lower:upper))
    q = nextprime(rand(RandomDevice(), lower:upper))

    n = p * q
    # lambda = lcm(p-1, q-1)
    d = invmod(e, (p - 1) * (q - 1))

    public_key = (n, e)
    private_key = (n, d)

    return (public_key, private_key)
end

sbox(bytes) = map(byte -> [
    0x63, 0x7c, 0x77, 0x7b, 0xf2, 0x6b, 0x6f, 0xc5, 0x30, 0x01, 0x67, 0x2b, 0xfe, 0xd7, 0xab, 0x76,
    0xca, 0x82, 0xc9, 0x7d, 0xfa, 0x59, 0x47, 0xf0, 0xad, 0xd4, 0xa2, 0xaf, 0x9c, 0xa4, 0x72, 0xc0,
    0xb7, 0xfd, 0x93, 0x26, 0x36, 0x3f, 0xf7, 0xcc, 0x34, 0xa5, 0xe5, 0xf1, 0x71, 0xd8, 0x31, 0x15,
    0x04, 0xc7, 0x23, 0xc3, 0x18, 0x96, 0x05, 0x9a, 0x07, 0x12, 0x80, 0xe2, 0xeb, 0x27, 0xb2, 0x75,
    0x09, 0x83, 0x2c, 0x1a, 0x1b, 0x6e, 0x5a, 0xa0, 0x52, 0x3b, 0xd6, 0xb3, 0x29, 0xe3, 0x2f, 0x84,
    0x53, 0xd1, 0x00, 0xed, 0x20, 0xfc, 0xb1, 0x5b, 0x6a, 0xcb, 0xbe, 0x39, 0x4a, 0x4c, 0x58, 0xcf,
    0xd0, 0xef, 0xaa, 0xfb, 0x43, 0x4d, 0x33, 0x85, 0x45, 0xf9, 0x02, 0x7f, 0x50, 0x3c, 0x9f, 0xa8,
    0x51, 0xa3, 0x40, 0x8f, 0x92, 0x9d, 0x38, 0xf5, 0xbc, 0xb6, 0xda, 0x21, 0x10, 0xff, 0xf3, 0xd2,
    0xcd, 0x0c, 0x13, 0xec, 0x5f, 0x97, 0x44, 0x17, 0xc4, 0xa7, 0x7e, 0x3d, 0x64, 0x5d, 0x19, 0x73,
    0x60, 0x81, 0x4f, 0xdc, 0x22, 0x2a, 0x90, 0x88, 0x46, 0xee, 0xb8, 0x14, 0xde, 0x5e, 0x0b, 0xdb,
    0xe0, 0x32, 0x3a, 0x0a, 0x49, 0x06, 0x24, 0x5c, 0xc2, 0xd3, 0xac, 0x62, 0x91, 0x95, 0xe4, 0x79,
    0xe7, 0xc8, 0x37, 0x6d, 0x8d, 0xd5, 0x4e, 0xa9, 0x6c, 0x56, 0xf4, 0xea, 0x65, 0x7a, 0xae, 0x08,
    0xba, 0x78, 0x25, 0x2e, 0x1c, 0xa6, 0xb4, 0xc6, 0xe8, 0xdd, 0x74, 0x1f, 0x4b, 0xbd, 0x8b, 0x8a,
    0x70, 0x3e, 0xb5, 0x66, 0x48, 0x03, 0xf6, 0x0e, 0x61, 0x35, 0x57, 0xb9, 0x86, 0xc1, 0x1d, 0x9e,
    0xe1, 0xf8, 0x98, 0x11, 0x69, 0xd9, 0x8e, 0x94, 0x9b, 0x1e, 0x87, 0xe9, 0xce, 0x55, 0x28, 0xdf,
    0x8c, 0xa1, 0x89, 0x0d, 0xbf, 0xe6, 0x42, 0x68, 0x41, 0x99, 0x2d, 0x0f, 0xb0, 0x54, 0xbb, 0x16,
][byte+1], bytes)  # Julia is 1-indexed

invsbox(bytes) = map(byte -> [
    0x52, 0x09, 0x6a, 0xd5, 0x30, 0x36, 0xa5, 0x38, 0xbf, 0x40, 0xa3, 0x9e, 0x81, 0xf3, 0xd7, 0xfb,
    0x7c, 0xe3, 0x39, 0x82, 0x9b, 0x2f, 0xff, 0x87, 0x34, 0x8e, 0x43, 0x44, 0xc4, 0xde, 0xe9, 0xcb,
    0x54, 0x7b, 0x94, 0x32, 0xa6, 0xc2, 0x23, 0x3d, 0xee, 0x4c, 0x95, 0x0b, 0x42, 0xfa, 0xc3, 0x4e,
    0x08, 0x2e, 0xa1, 0x66, 0x28, 0xd9, 0x24, 0xb2, 0x76, 0x5b, 0xa2, 0x49, 0x6d, 0x8b, 0xd1, 0x25,
    0x72, 0xf8, 0xf6, 0x64, 0x86, 0x68, 0x98, 0x16, 0xd4, 0xa4, 0x5c, 0xcc, 0x5d, 0x65, 0xb6, 0x92,
    0x6c, 0x70, 0x48, 0x50, 0xfd, 0xed, 0xb9, 0xda, 0x5e, 0x15, 0x46, 0x57, 0xa7, 0x8d, 0x9d, 0x84,
    0x90, 0xd8, 0xab, 0x00, 0x8c, 0xbc, 0xd3, 0x0a, 0xf7, 0xe4, 0x58, 0x05, 0xb8, 0xb3, 0x45, 0x06,
    0xd0, 0x2c, 0x1e, 0x8f, 0xca, 0x3f, 0x0f, 0x02, 0xc1, 0xaf, 0xbd, 0x03, 0x01, 0x13, 0x8a, 0x6b,
    0x3a, 0x91, 0x11, 0x41, 0x4f, 0x67, 0xdc, 0xea, 0x97, 0xf2, 0xcf, 0xce, 0xf0, 0xb4, 0xe6, 0x73,
    0x96, 0xac, 0x74, 0x22, 0xe7, 0xad, 0x35, 0x85, 0xe2, 0xf9, 0x37, 0xe8, 0x1c, 0x75, 0xdf, 0x6e,
    0x47, 0xf1, 0x1a, 0x71, 0x1d, 0x29, 0xc5, 0x89, 0x6f, 0xb7, 0x62, 0x0e, 0xaa, 0x18, 0xbe, 0x1b,
    0xfc, 0x56, 0x3e, 0x4b, 0xc6, 0xd2, 0x79, 0x20, 0x9a, 0xdb, 0xc0, 0xfe, 0x78, 0xcd, 0x5a, 0xf4,
    0x1f, 0xdd, 0xa8, 0x33, 0x88, 0x07, 0xc7, 0x31, 0xb1, 0x12, 0x10, 0x59, 0x27, 0x80, 0xec, 0x5f,
    0x60, 0x51, 0x7f, 0xa9, 0x19, 0xb5, 0x4a, 0x0d, 0x2d, 0xe5, 0x7a, 0x9f, 0x93, 0xc9, 0x9c, 0xef,
    0xa0, 0xe0, 0x3b, 0x4d, 0xae, 0x2a, 0xf5, 0xb0, 0xc8, 0xeb, 0xbb, 0x3c, 0x83, 0x53, 0x99, 0x61,
    0x17, 0x2b, 0x04, 0x7e, 0xba, 0x77, 0xd6, 0x26, 0xe1, 0x69, 0x14, 0x63, 0x55, 0x21, 0x0c, 0x7d,
][byte+1], bytes)

"""
Generate the subkeys for the AES encryption standard given a key.
    https://embeddedsw.net/Cipher_Reference_Home.html
    https://en.wikipedia.org/wiki/AES_key_schedule
    https://cryptography.fandom.com/wiki/Rijndael_key_schedule
"""
function rijndael_key_schedule(key::Vector{UInt8}, r)
    rot_word(bytes::Vector{UInt8}) = UInt8[bytes[2], bytes[3], bytes[4], bytes[1]]

    function genrcon(n::Integer)  # no longer needed
        result = UInt8[0x01]
        for i in 2:n
            rc_i_minus_one = result[end]
            rc_i = 0x02*rc_i_minus_one
            if rc_i_minus_one < 0x80
                push!(result, rc_i)
            else
                push!(result, xor(rc_i, 0x1b))
            end
        end
        return result
    end

    function nextrcon(rc_i_minus_one)
        rc_i = 0x02*rc_i_minus_one
        if rc_i_minus_one < 0x80
            return rc_i
        else
            return xor(rc_i, 0x1b)
        end
    end

    w = [key[i:(i+3)] for i in 1:4:length(key)]
    n = length(w)
    #rcon = genrcon(4*(r+1) % 256)  # rcon_i = x^(254+i) mod x^8 + x^4 + x^3 + x + 1
    rcon = 0x01

    for i in (n+1):(4*(r+1))
        if (i-1) % n == 0
            push!(w, xor.(w[i-n], sbox(rot_word(w[i-1])), [rcon, 0x00, 0x00, 0x00]))
            rcon = nextrcon(rcon)
        elseif (i-1)% n == 4 && n > 6
            push!(w, xor.(w[i-n], sbox(w[i-1])))
        else
            push!(w, xor.(w[i-n], w[i-1]))
        end
    end

    return w
end


text:
using DataStructures

Zygote.@adjoint! function heappush!(xs, x, ord)
  heappush!(xs, x, ord), Δ -> begin
    (nothing, nothing, nothing)
  end
end

Zygote.@adjoint! function heappop!(xs, ord)
  heappop!(xs, ord), Δ -> begin
    (nothing, nothing)
  end
end


text:
"""
$(SIGNATURES)

Maxwell's diffusive boundary flux

- @args: distribution functions and their slopes at left/right sides of interface
- @args: particle velocity quadrature points and weights
- @args: time step

1D continuum
"""
function flux_boundary_maxwell!(
    fw::AV{T1},
    bc::AV{T2},
    w::AV{T3},
    inK,
    γ,
    dt,
    rot, # 1 / -1
) where {T1<:FN,T2<:Real,T3<:Real}

    @assert length(bc) == 3

    primL, primR = ifelse(rot == 1, (bc, conserve_prim(w, γ)), (conserve_prim(w, γ), bc))

    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    Muv_L = moments_conserve(MuL1, Mxi1, 1, 0)
    Muv_R = moments_conserve(MuR2, Mxi2, 1, 0)

    ρ = ifelse(rot == 1, -primR[1] * Muv_R[1] / Muv_L[1], -primL[1] * Muv_L[1] / Muv_R[1])

    @. fw = ρ * (Muv_L + Muv_R) * dt

    return nothing

end

"""
$(SIGNATURES)

1D2F1V
"""
function flux_boundary_maxwell!(
    fw::AV{T1},
    fh::T2,
    fb::T2,
    bc::AV{T3},
    h::T4,
    b::T4,
    u::T5,
    ω::T5,
    inK,
    dt,
    rot = 1,
) where {T1<:FN,T2<:AV{<:FN},T3<:Real,T4<:AV{<:FN},T5<:AV{<:FN}} # 1D2F1V

    @assert length(bc) == 3

    δ = heaviside.(u .* rot)
    SF = sum(ω .* u .* h .* (1.0 .- δ))
    SG = (bc[end] / π)^0.5 * sum(ω .* u .* exp.(-bc[end] .* (u .- bc[2]) .^ 2) .* δ)
    prim = [-SF / SG; bc[2:end]]

    H = maxwellian(u, prim)
    B = H .* inK ./ (2.0 * prim[end])

    hWall = H .* δ .+ h .* (1.0 .- δ)
    bWall = B .* δ .+ b .* (1.0 .- δ)

    fw[1] = discrete_moments(hWall, u, ω, 1) * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * dt
    fw[3] =
        (
            0.5 * discrete_moments(hWall .* u .^ 2, u, ω, 1) +
            0.5 * discrete_moments(bWall, u, ω, 1)
        ) * dt

    @. fh = u * hWall * dt
    @. fb = u * bWall * dt

    return nothing

end

"""
$(SIGNATURES)

2D Continuum
"""
function flux_boundary_maxwell!(
    fw::AV{T1},
    bc::AV{T2},
    w::AV{T3},
    inK,
    γ,
    dt,
    len,
    rot,
) where {T1<:Real,T2<:Real,T3<:Real}

    @assert length(bc) == 4

    primL, primR = ifelse(rot == 1, (bc, conserve_prim(w, γ)), (conserve_prim(w, γ), bc))

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    Muv_L = moments_conserve(MuL1, Mv1, Mxi1, 1, 0, 0)
    Muv_R = moments_conserve(MuR2, Mv2, Mxi2, 1, 0, 0)

    ρ = ifelse(rot == 1, -primR[1] * Muv_R[1] / Muv_L[1], -primL[1] * Muv_L[1] / Muv_R[1])

    @. fw = ρ * (Muv_L + Muv_R) * dt * len

    return nothing

end

"""
$(SIGNATURES)

2D1F2V
"""
function flux_boundary_maxwell!(
    fw::AV{T1},
    fh::AM{T2},
    bc::AV{T3},
    h::AM{T4},
    u::T5,
    v::T5,
    ω::T5,
    dt,
    len,
    rot = 1,
) where {T1<:FN,T2<:FN,T3<:Real,T4<:FN,T5<:AM{<:FN}}

    @assert length(bc) == 4

    δ = heaviside.(u .* rot)

    SF = sum(ω .* u .* h .* (1.0 .- δ))
    SG =
        (bc[end] / π) *
        sum(ω .* u .* exp.(-bc[end] .* ((u .- bc[2]) .^ 2 .+ (v .- bc[3]) .^ 2)) .* δ)
    prim = [-SF / SG; bc[2:end]]

    H = maxwellian(u, v, prim)

    hWall = H .* δ .+ h .* (1.0 .- δ)

    fw[1] = discrete_moments(hWall, u, ω, 1) * len * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * len * dt
    fw[3] = discrete_moments(hWall .* u, v, ω, 1) * len * dt
    fw[4] = 0.5 * discrete_moments(hWall .* (u .^ 2 .+ v .^ 2), u, ω, 1) * len * dt

    @. fh = u * hWall * len * dt

    return nothing

end

"""
$(SIGNATURES)

2D2F2V
"""
function flux_boundary_maxwell!(
    fw::AV{T1},
    fh::T2,
    fb::T2,
    bc::AV{T3},
    h::T4,
    b::T4,
    u::T5,
    v::T5,
    ω::T5,
    inK,
    dt,
    len,
    rot = 1,
) where {T1<:FN,T2<:AM{<:FN},T3<:Real,T4<:AM{<:FN},T5<:AM{<:FN}}

    @assert length(bc) == 4

    δ = heaviside.(u .* rot)

    SF = sum(ω .* u .* h .* (1.0 .- δ))
    SG =
        (bc[end] / π) *
        sum(ω .* u .* exp.(-bc[end] .* ((u .- bc[2]) .^ 2 .+ (v .- bc[3]) .^ 2)) .* δ)
    prim = [-SF / SG; bc[2:end]]

    H = maxwellian(u, v, prim)
    B = H .* inK ./ (2.0 * prim[end])

    hWall = H .* δ .+ h .* (1.0 .- δ)
    bWall = B .* δ .+ b .* (1.0 .- δ)

    fw[1] = discrete_moments(hWall, u, ω, 1) * len * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * len * dt
    fw[3] = discrete_moments(hWall .* u, v, ω, 1) * len * dt
    fw[4] =
        (
            0.5 * discrete_moments(hWall .* (u .^ 2 .+ v .^ 2), u, ω, 1) +
            0.5 * discrete_moments(bWall, u, ω, 1)
        ) *
        len *
        dt

    @. fh = u * hWall * len * dt
    @. fb = u * bWall * len * dt

    return nothing

end

"""
$(SIGNATURES)

1F3V
"""
function flux_boundary_maxwell!(
    fw::AV{T1},
    ff::AA{T2,3},
    bc::AV{T3},
    f::AA{T4,3},
    u::T5,
    v::T5,
    w::T5,
    ω::T5,
    dt,
    area,
    rot = 1,
) where {T1<:FN,T2<:FN,T3<:Real,T4<:FN,T5<:AA{<:FN,3}}

    @assert length(bc) == 5

    δ = heaviside.(u .* rot)
    SF = sum(ω .* u .* f .* (1.0 .- δ))
    SG =
        (bc[end] / π)^1.5 * sum(
            ω .* u .*
            exp.(
                -bc[end] .* ((u .- bc[2]) .^ 2 .+ (v .- bc[3]) .^ 2 .+ (w .- bc[4]) .^ 2),
            ) .* δ,
        )
    prim = [-SF / SG; bc[2:end]]

    M = maxwellian(u, v, w, prim)
    fWall = M .* δ .+ f .* (1.0 .- δ)

    fw[1] = discrete_moments(fWall, u, ω, 1) * area * dt
    fw[2] = discrete_moments(fWall, u, ω, 2) * area * dt
    fw[3] = discrete_moments(fWall .* u, v, ω, 1) * area * dt
    fw[4] = discrete_moments(fWall .* u, w, ω, 1) * area * dt
    fw[5] =
        (0.5 * discrete_moments(fWall .* (u .^ 2 .+ v .^ 2 + w .^ 2), u, ω, 1)) * area * dt

    @. ff = u * fWall * area * dt

    return nothing

end

"""
$(SIGNATURES)

Specular reflection boundary flux

- @args: fluxes of conservative variables and distribution functions
- @args: distribution functions
- @args: velocity quadrature points and weights
- @args: time step
"""
function flux_boundary_specular!(
    fw::AV{T1},
    ff::AV{T2},
    f::AV{T3},
    u::T4,
    ω::T4,
    dt,
) where {T1<:Real,T2<:FN,T3<:FN,T4<:AV{<:FN}} # 1D1F1V

    fWall = similar(f)
    for i in eachindex(f)
        fWall[i] = f[end-i+1]
    end

    fw[1] = discrete_moments(fWall, u, ω, 1) * dt
    fw[2] = discrete_moments(fWall, u, ω, 2) * dt
    fw[3] = 0.5 * discrete_moments(fWall .* u .^ 2, u, ω, 1) * dt

    @. ff = u * fWall * dt

    return nothing

end

"""
$(SIGNATURES)

1D2F1V
"""
function flux_boundary_specular!(
    fw::AV{T1},
    fh::T2,
    fb::T2,
    h::T3,
    b::T3,
    u::T4,
    ω::T4,
    dt,
) where {T1<:Real,T2<:AV{<:FN},T3<:AV{<:FN},T4<:AV{<:FN}}

    hWall = similar(h)
    bWall = similar(b)
    for i in eachindex(h)
        hWall[i] = h[end-i+1]
        bWall[i] = b[end-i+1]
    end

    fw[1] = discrete_moments(hWall, u, ω, 1) * dt
    fw[2] = discrete_moments(hWall, u, ω, 2) * dt
    fw[3] =
        (
            0.5 * discrete_moments(hWall .* u .^ 2, u, ω, 1) +
            0.5 * discrete_moments(bWall, u, ω, 1)
        ) * dt

    @. fh = u * hWall * dt
    @. fb = u * bWall * dt

    return nothing

end


text:
# The mglplot macro which can be used similar to the mgl scripting language
macro mglplot(graph, commands)
    # Check if the format provided is right
    if typeof(graph)    != Symbol ||
       typeof(commands) != Expr
        Base.error("Arguments have wrong type. Need (Symbol, Expr),"  *
              "given was ($(typeof(graph)), $(typeof(commands)))")
    end
    if commands.head != :vcat
        Base.error("Can't understand the sequence of mgl commands. " *
                   "Expected symbol ':vcat'.")
    end
    # Initialize empty quote block that is to be filled with mgl calls
    block = Expr(:block)
    # Go through every line (= row) in the given block
    # Change the function name from 'name' to 'MathGL.name'
    # And append the row to the block to be executed
    for row in commands.args
        if typeof(row) == Symbol
            fname = Expr(:., :MathGL, Expr(:quote, :($(row))))
            push!(block.args, Expr(:call, fname, graph))
        elseif typeof(row) != Expr # Do nothing now, maybe TODO
        elseif row.head != :row # Just do what is given as statement, if it is no row
            push!(block.args, row)
        else
            # Replace ":(=)" symbols with ":kw" symbols so that they can be put in a function
            args = [ typeof(arg)    != Expr ?
                     arg : arg.head != :(=) ?
                     arg : Expr(:kw, arg.args...)
                     for arg in row.args          ]

            fname = Expr(:., :MathGL, Expr( :quote, :($(args[1]))))
            push!(block.args, Expr(:call, fname, graph, args[2:end]...))
        end
    end
    return esc(block)
end

export @mglplot


text:
using Combinatorics
include("intcode.jl")

"""Compute the thruster signal given prog to run on amplifiers and phase setting list"""
function thruster_signal(prog, phase_settings)
    previous_output = 0
    for phase_setting in phase_settings
        previous_output = first(run_intcode!(copy(prog), [phase_setting, previous_output]))
    end
    return previous_output
end

"""Compute the thruster signal in feedback mode after all programs halt"""
function feedback_signal(prog, phase_settings)
    io_channels = [Channel(4) for i in 1:5]
    for (io_channel, phase_setting) in zip(io_channels, phase_settings) put!(io_channel, phase_setting) end
    put!(first(io_channels), 0)
    @sync for i in 1:5
        @async run_intcode!(copy(prog), io_channels[i], io_channels[i%5 + 1])
    end
    return take!(first(io_channels))
end

# Part 1: run the amplifier program with all permutations of 01234 as phase settings
prog = open(f -> parse.(Int32, split(readline(f), ",")), "day07_input.txt")
signals = [thruster_signal(prog, phase_settings) for phase_settings in permutations(0:4)]
println("Highest thruster signal: $(maximum(signals))")

# Part 2: this time we keep running in feedback mode
signals = [feedback_signal(prog, phase_settings) for phase_settings in permutations(5:9)]
println("Highest feedback thruster signal: $(maximum(signals))")

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
using Test
using SumOfSquares
using DynamicPolynomials

function rearrangement_test(optimizer, config::MOIT.Config)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @polyvar x y z
    # If x, y is in the same order as y, z, then by the rearrangement inequality
    # xy + yz ≤ xz + y^2
    # Then, if x and z are nonnegative, by the Geometric Mean-Quadratic Mean inequality,
    # xz ≤ (x^2 + z^2)/2
    # Therefore,
    # 2xy + 2yz ≤ x^2 + z^2 + 2y^2 over the set:
    S = @set y ≤ x && y ≤ z && x ≥ 0 && z ≥ 0
    # In fact, it is nonnegative everywhere so let's add terms `x + z`:
    cref = @constraint(model, x^2 + z^2 + 2y^2 - 2x * y - 2y * z + x + z in SOSCone(),
                       domain = S, sparsity = Sparsity.Variable(), maxdegree = 3)

    optimize!(model)

    @test termination_status(model) == MOI.OPTIMAL

    @test primal_status(model) == MOI.FEASIBLE_POINT

    p = gram_matrix(cref)
    @test p isa SumOfSquares.SparseGramMatrix
    @test length(p.sub_gram_matrices) == 2
    @test getmat(p.sub_gram_matrices[1]) ≈ [1 -1 0; -1 1 0; 0 0 0] atol=9atol rtol=9rtol
    @test p.sub_gram_matrices[1].basis.monomials == [y, z, 1]
    @test getmat(p.sub_gram_matrices[2]) ≈ [1 -1 0; -1 1 0; 0 0 0] atol=9atol rtol=9rtol
    @test p.sub_gram_matrices[2].basis.monomials == [x, y, 1]

    λ = lagrangian_multipliers(cref)
    @test length(λ) == 4
    @test λ[1] isa SumOfSquares.SparseGramMatrix
    @test length(λ[1].sub_gram_matrices) == 1
    @test λ[1].sub_gram_matrices[1].basis.monomials == [x, y, 1]
    @test λ[2] isa SumOfSquares.SparseGramMatrix
    @test length(λ[2].sub_gram_matrices) == 1
    @test λ[2].sub_gram_matrices[1].basis.monomials == [y, z, 1]
    @test λ[3] isa SumOfSquares.SparseGramMatrix
    @test length(λ[3].sub_gram_matrices) == 1
    @test λ[3].sub_gram_matrices[1].basis.monomials == [x, y, 1]
    @test λ[4] isa SumOfSquares.SparseGramMatrix
    @test length(λ[4].sub_gram_matrices) == 1
    @test λ[4].sub_gram_matrices[1].basis.monomials == [y, z, 1]

    ν = moment_matrix(cref)
    @test ν isa SparseMomentMatrix
    @test length(ν.sub_moment_matrices) == 2
    @test ν.sub_moment_matrices[1].basis.monomials == [y, z, 1]
    @test ν.sub_moment_matrices[2].basis.monomials == [x, y, 1]
end
sd_tests["rearrangement"] = rearrangement_test

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
#=
do:
- Julia version:
- Author: anderscs
- Date: 2018-11-07
=#
using PyCall

pushfirst!(PyVector(pyimport("sys")["path"]), "")

models = pyimport("odin.models")
def = pyimport("odin.utils.default")
actions = pyimport("odin.actions")

model_wrapper = models.load_model("pix2pix")
#args, model_wrapper = def.default_chainer()

actions.actions["train_model"]()

max_stars_repo_name: Anderssorby/odin

max_issues_repo_name: Anderssorby/odin

max_forks_repo_name: Anderssorby/odin


text:
using Distributed

# Thea: set this path to point to the correct module file
include("SeisNoiseAWS.jl")
@everywhere include("SeisNoiseAWS.jl")

# Thea: if you want to run tests, you can hard code a path here to a data folder that contains params.csv and stations.csv
default_path = "/Users/dmikesell/Documents/GitHub/esip_julia_aws/dev/20200615_102819/41808/"

# Example of the CLI call of this Julia file.
#
# $ julia run_SeisNoiseAWS.jl /Users/dmikesell/Documents/GitHub/esip_julia_aws/dev/20200615_102819/41808/
#
# or using multiple threads (4 in this example)
#
# $ julia -p 4 run_SeisNoiseAWS.jl /Users/dmikesell/Documents/GitHub/esip_julia_aws/dev/20200615_102819/41808/

# ==============================================================================
# For testing in Atom without calling from Linux CLI.
# clearconsole()
if isempty(ARGS)
    basepath = default_path
else
    basepath = ARGS[1]
end
# ==============================================================================
using Logging
# Setup the log file
log_file = joinpath(basepath, "xcor_log.txt")
io = open(log_file, "w")
logger = SimpleLogger(io, Logging.Info)
global_logger(logger) # Set the global logger to logger
# ==============================================================================
# ncores = length(Sys.cpu_info())
# ntheads = Threads.nthreads()
# addprocs(ntheads) # start DEFAULT parallel threads
# addprocs(4)
num_procs = nprocs()
if num_procs == 1
    @info("Number of workers: $num_procs")
else
    @info("Number of workers: $(num_procs-1)", )
end
# ==============================================================================
using Dates # even though loaded in SeisNoiseAWS module, needs to be loaded here too.
using SeisIO
using Combinatorics
# ==============================================================================
# Load the input files from the web application
#
# Build the COR output directory
COR_DIR = SeisNoiseAWS.build_cor_directory(basepath)
# Get Station information
df_st = SeisNoiseAWS.load_station_file( joinpath(basepath, "stations.csv") )
# Get XCOR parameters
df_params = SeisNoiseAWS.load_parameter_file( joinpath(basepath, "params.csv") )
# ==============================================================================
# Build the data availability table
df_date = SeisNoiseAWS.build_availability_df(df_st)
# Make boolean vectors from each column
columns = collect( eachcol(df_date) )
# Set the file length for downloading data
t_duration = Dates.Second(Dates.Day(1)).value # [s] one day in seconds
# ==============================================================================
# Loop over days
t0a = now() # start timer
# for col = 7:7#20 # end-1 # skip the first column, which is channel tag
col = 10
# col=2 set up some if statements about num_pairs
# col=8 has a trace that is shorter than the cc_len (so need try/catch)
# col=7 has only a single station (the try/catch from above works)

# Get the date information
the_date    = String(names(df_date)[col  ])  # start date
the_date_p1 = String(names(df_date)[col+1])  # end date
@info("Processing: $the_date to $the_date_p1")

# Create the output folder for this date
folder_date = Dates.format(DateTime(String(names(df_date)[col])),"yyyy-mm-dd")
output_dir = joinpath(COR_DIR, folder_date)
@debug("Writing daily correlations to $output_dir")

# Download the station data for this date
stations = columns[1][columns[col]] # array of channel_tags on this day
@debug("Downloading data for $stations")
S = SeisData() # create empty seisdata object
# Download data if more than 1 channel
if length(stations) > 1
    # This does the instrument deconvolution on the fly
    get_data!(S, "FDSN", stations, src="IRIS", s=the_date, t=t_duration, rr=true)
end

# Build the cross-correlation pairs
xcor_pairs = collect(combinations(S, 2))

# We need a flag to build the autocorrelation pairs
# Then we need to create those and append to xcor_pairs

# Do the pair-wise cross correlation
if isempty(xcor_pairs)
    @info("Number of channels downloaded: $(length(S))")
    @info("No channel pairs to correlate.")
elseif any(x->x=="XX.FAIL..001", S.id)
    @info("Bad station request: nothing to correlate.")
else
    numpairs = length(xcor_pairs)
    @debug("Number of channel pair correlations: $numpairs")
    t1 = now()  # start timer for the parallel computation for this day
    # Do the channel pair correlations in parallel for a single day
    pmap(SeisNoiseAWS.xcor,
        xcor_pairs,
        fill(df_params, numpairs),
        fill(output_dir, numpairs))
    # SeisNoiseAWS.xcor(xcor_pairs, df_params, output_dir)

    cor_time = now() - t1 # compute time for this day
    @debug("Compute time for this correlation day: $cor_time")
end # if statement over xcor_pairs

# end # the FOR loop over days

run_time = now() - t0a # total run time timer
@info("Compute time for all correlations: $run_time")

close(io) # stop logger

exit()


text:
# This file is a part of ForwardDiffPullbacks.jl, licensed under the MIT License (MIT).

using ForwardDiffPullbacks
using Test

using Distributions
using ForwardDiff, ChainRulesCore, Zygote

include("testfuncs.jl")

@testset "Zygote" begin
    function cr_fwd_and_back(f, xs, ΔΩ)
        y, back = ChainRulesCore.rrule(f, xs...)
        back_thunks = back(ΔΩ)
        Δx = map(unthunk, back_thunks)
        @assert first(Δx) == NoTangent()
        y, Base.tail(Δx)
    end

    function zg_fwd_and_back(f, xs, ΔΩ)
        y, back = Zygote.pullback(f, xs...)
        Δx = back(ΔΩ)
        y, Δx
    end

    function cr_bc_fwd_and_back(f, Xs, ΔΩA)
        y, back = ChainRulesCore.rrule(Base.broadcasted, f, Xs...)
        back_thunks = back(ΔΩA)
        Δx = map(unthunk, back_thunks)
        @assert first(Δx) == NoTangent()
        y, Base.tail(Δx)
    end

    function zg_bc_fwd_and_back(f, Xs, ΔΩA)
        y, back = Zygote.pullback((Xs...) -> f.(Xs...), Xs...)
        Δx = back(ΔΩA)
        y, Δx
    end
    function zg_bcpre_fwd_and_back(f, Xs, ΔΩA)
        y, back = Zygote.pullback(Base.broadcasted, f, Xs...)
        Δx = back(ΔΩA)
        y, Δx
    end

    @test @tinferred(cr_fwd_and_back(fwddiff(f), xs, ΔΩ)) isa Tuple{Tuple{Float32, Float32, Int}, Tuple{Float32, Tuple{Float32, Float32}, SVector{3, Float32}}}
    @test @tinferred(zg_fwd_and_back(fwddiff(f), xs, ΔΩ)) isa Tuple{Tuple{Float32, Float32, Int}, Tuple{Float32, Tuple{Float32, Float32}, SVector{3, Float32}}}

    @test cr_fwd_and_back(fwddiff(f), xs, ΔΩ) == ((139, 783, 42), (280, (600, 1040), SVector(1600, 2280, 3080)))
    @test zg_fwd_and_back(fwddiff(f), xs, ΔΩ) == ((139, 783, 42), (280, (600, 1040), SVector(1600, 2280, 3080))) # == zg_fwd_and_back(f, xs, ΔΩ)

    function f_loss_1(xs...)
        r = fwddiff(f)(xs...)
        @assert sum(r) < 10000
        sum(r[1])
    end

    function f_loss_3(xs...)
        r = fwddiff(f)(xs...)
        @assert sum(r) < 10000
        sum(r[3])
    end

    function f_loss_3z(xs...)
        r = f(xs...)
        @assert sum(r) < 10000
        sum(r[3])
    end

    Xs = map(x -> fill(x, 5), xs)
    ΔΩA = fill(ΔΩ, 5)

    @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, (Xs[1], Ref(xs[2]), Xs[1]), Val(3), ΔΩA)) == fill(280, 5)
    @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, (Xs[1], Ref(xs[2]), Xs[2]), Val(3), ΔΩA)) == fill((600, 1040), 5)
    @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, (Xs[1], Ref(xs[2]), Xs[3]), Val(3), ΔΩA)) == fill(SVector(1600, 2280, 3080), 5)

    for args in (Xs, (Xs[1], Ref(xs[2]), Xs[3]), map(Ref, xs))
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(1), ΔΩA)) == fill(280, 5)
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(2), ΔΩA)) == fill((600, 1040), 5)
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(3), ΔΩA)) == fill(SVector(1600, 2280, 3080), 5)
    end

    for args in (Xs, (Xs[1], Ref(xs[2]), Xs[3]))
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(1), Ref(ΔΩ))) == fill(280, 5)
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(2), Ref(ΔΩ))) == fill((600, 1040), 5)
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(3), Ref(ΔΩ))) == fill(SVector(1600, 2280, 3080), 5)
    end

    let args = map(Ref, xs), ΔY = Ref(ΔΩ)
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(1), Ref(ΔΩ))) == 280
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(2), Ref(ΔΩ))) == (600, 1040)
        @test @tinferred(ForwardDiffPullbacks.forwarddiff_bc_fwd_back(f, args, Val(3), Ref(ΔΩ))) == SVector(1600, 2280, 3080)
    end

    # @tinferred cr_bc_fwd_and_back(fwddiff(f), Xs, ΔΩA)
    # zg_bc_fwd_and_back(fwddiff(f), Xs, ΔΩA)
    # zg_bc_fwd_and_back(f, Xs, ΔΩA)

    @testset "disttrafo" begin
        function disttrafo(trg_dist::Distribution, src_dist::Distribution, x::Real, prev_ladj::Union{Real,Missing})
            # @info "disttrafo for real x and optional prev_ladj"
            u = cdf(src_dist, x)
            y = quantile(trg_dist, u)
            ladj = if isnothing(prev_ladj)
                nothing
            else
                logladj_1 = + logpdf(src_dist, x)
                logladj_2 = - logpdf(trg_dist, y)
                logladj_1 + logladj_2 + prev_ladj
            end
            (y = y, ladj = ladj)
        end

        function disttrafo(trg_dist::Distribution, src_dist::Distribution, x_dual::ForwardDiff.Dual{TAG}, prev_ladj::Missing) where TAG
            # @info "disttrafo for dual x and missing prev_ladj"
            x = ForwardDiff.value(x_dual)
            u = cdf(src_dist, x)
            y = quantile(trg_dist, u)
            dudx = pdf(src_dist, x)
            dydu = inv(pdf(trg_dist, y))
            y_dual = ForwardDiff.Dual{TAG}(x, dydu * dudx * ForwardDiff.partials(x_dual))
            (y = y_dual, ladj = nothing)
        end

        # Must work:
        ForwardDiff.derivative(x -> disttrafo(Beta(), Normal(), x, missing).y, -0.5)

        # Requires Zygote to utilize thunks (https://github.com/FluxML/Zygote.jl/pull/966):
        if isdefined(Zygote, Symbol("@_adjoint_keepthunks"))
            Zygote.gradient(x -> fwddiff(disttrafo)(Beta(), Normal(), x, missing).y, -0.5)
        end

        # Expected to fail for Beta():
        # Zygote.gradient(x -> disttrafo(Beta(), Normal(), x, missing).y, x)

        # disttrafo(trg_dist::Distribution, src_dist::Distribution, x::Dual, prev_ladj::Missing) = (y = quantile(trg_dist, cdf(src_dist, x)), ladj = missing)

        trg_d = Weibull()
        src_d = Normal()
        x = 0.5

        ΔΩ = (y = 1, ladj = nothing)
        #=@tinferred=# cr_fwd_and_back(fwddiff(disttrafo), (trg_d, src_d, x, missing), ΔΩ)
        #=@tinferred=# zg_fwd_and_back(fwddiff(disttrafo), (trg_d, src_d, x, missing), ΔΩ)
        #compare with (ignore missings):
        zg_fwd_and_back(disttrafo, (trg_d, src_d, x, missing), ΔΩ)

        ΔΩ = (y = 1, ladj = 42)
        #=@tinferred=# cr_fwd_and_back(fwddiff(disttrafo), (trg_d, src_d, x, 7), ΔΩ)
        #=@tinferred=# zg_fwd_and_back(fwddiff(disttrafo), (trg_d, src_d, x, 7), ΔΩ)
        #compare with (use deep approx):
        zg_fwd_and_back(disttrafo, (trg_d, src_d, x, 7), ΔΩ)

        n = 100
        trg_D = fill(trg_d, n)
        src_D = fill(src_d, n)
        X = randn(n)

        ΔΩs = fill((y = 1, ladj = nothing), n)
        #=@tinferred=#(cr_bc_fwd_and_back(fwddiff(disttrafo), (trg_D, src_D, X, missing), ΔΩs))
        #=@tinferred=#(zg_bc_fwd_and_back(fwddiff(disttrafo), (trg_D, src_D, X, missing), ΔΩs))
        #compare with (ignore missings):
        zg_bc_fwd_and_back(disttrafo, (trg_D, src_D, X, missing), ΔΩs)

        ΔΩs = fill((y = 1, ladj = 42), n)
        #=@tinferred=#(cr_bc_fwd_and_back(fwddiff(disttrafo), (trg_D, src_D, X, 7), ΔΩs))
        #=@tinferred=#(zg_bc_fwd_and_back(fwddiff(disttrafo), (trg_D, src_D, X, 7), ΔΩs))
        #compare with (use deep approx):
        zg_bc_fwd_and_back(disttrafo, (trg_D, src_D, X, 7), ΔΩs)

        #=
        using BenchmarkTools
        @benchmark disttrafo.($trg_D, $src_D, $X, $7)
        @benchmark cr_bc_fwd_and_back(fwddiff(disttrafo), ($trg_D, $src_D, $X, 7), $ΔΩs)
        @benchmark zg_bc_fwd_and_back(fwddiff(disttrafo), ($trg_D, $src_D, $X, 7), $ΔΩs)
        @benchmark zg_bc_fwd_and_back(disttrafo, ($trg_D, $src_D, $X, 7), $ΔΩs)
        =#

        function dummy_loss(trg_d::Distribution, src_d::Distribution, x::Real, prev_ladj::Union{Real,Missing})
            y, ladj = fwddiff(disttrafo)(trg_d, src_d, x, prev_ladj)
            ismissing(ladj) ? 2 * y : typeof(y)(y^2 + ladj^2)
        end

        #=@tinferred=# Zygote.gradient(dummy_loss, trg_d, src_d, x, 42)
        #=@tinferred=# Zygote.gradient(dummy_loss, trg_d, src_d, x, missing)

        function bc_dummy_loss(trg_D::AbstractVector{<:Distribution}, src_D::AbstractVector{<:Distribution}, X::AbstractVector{<:Real}, prev_ladj::Union{Real,Missing})
            Y_ladj = fwddiff(disttrafo).(trg_D, src_D, X, prev_ladj)
            Y = (x -> x.y).(Y_ladj)
            ladj = (x -> x.ladj).(Y_ladj)
            T = eltype(Y)
            any(ismissing, ladj) ? T(sum(Y)) : T(sum(Y)^2 + sum(ladj)^2)
        end

        #=@tinferred=# bc_dummy_loss(trg_D, src_D, X, 42)
        #=@tinferred=# bc_dummy_loss(trg_D, src_D, X, missing)

        #=@tinferred=# Zygote.gradient(bc_dummy_loss, trg_D, src_D, X, 42)
        #=@tinferred=# Zygote.gradient(bc_dummy_loss, trg_D, src_D, X, missing)

        #=@tinferred=# Zygote.gradient(X -> bc_dummy_loss(trg_D, src_D, X, 42), X)

        # Requires Zygote to utilize thunks (https://github.com/FluxML/Zygote.jl/pull/966):
        if isdefined(Zygote, Symbol("@_adjoint_keepthunks"))
            trg_D = fill(Beta(), n)
            Zygote.gradient(X -> bc_dummy_loss(trg_D, src_D, X, missing), X)
        end
    end
end


text:
"""
   check(cond, msg)

If `cond` evaluates false throw `ArgumentError` and print evaluation of `cond`.
"""
macro check(cond::Expr, msg)
    head = lstrip(repr(cond), ':')
    head = head * " == false"
    args = ()
    for (i,a) in enumerate(cond.args[2:end])
        lhs = lstrip(repr(a), ':')
        symbol = (i == length(cond.args)-1) ? "└ " : "├ "
        args  = (args..., :("\n   " * $symbol * $lhs * " = " * repr($(esc(a)))))
    end
    return :($(esc(cond)) ||
             throw(ArgumentError($msg * "\n  " * $head * $(args...))))
end

"""
    remove_namespace(namespace, x)

If first namespace of `x` is `namespace` then remove it.
```
remove_namespace(A, A₊B₊x) -> B₊x
remove_namespace(B, A₊B₊x) -> A₊B₊x
````
"""
remove_namespace(namespace, name::T) where T = T(replace(String(name), Regex("^$(namespace)₊") => ""))
remove_namespace(namespace, x::Sym) = rename(x, remove_namespace(namespace, x.name))
remove_namespace(namespace, x::Term) = rename(x, remove_namespace(namespace, operation(x).name))

"""
    remove_namespace(x)

Strips `x` of its first namespace. `A₊B₊x -> B₊x`
"""
remove_namespace(name::T) where T = T(replace(String(name), Regex("^(.+?)₊") => ""))
remove_namespace(x::Sym) = rename(x, remove_namespace(x.name))
remove_namespace(x::Term) = rename(x, remove_namespace(operation(x).name))

eqsubstitute(eq::Equation, rules) = substitute(eq.lhs, rules) ~ substitute(eq.rhs, rules)

uniquenames(syms) = allunique(getname.(syms))

"""
    eq_type(eq::Equation)

Checks the type of the equation. Returns:
- `(:diffeq, lhs_variable)` for differential equations
- `(:explicit_algebraic, lhs_variable)` for explicit algebraic equations
- `(:implicit_algebraic, lhs_variable)` for implicit algebraic equations
"""
function eq_type(eq::Equation)
    if eq.lhs isa Term && operation(eq.lhs) isa Differential
        vars = get_variables(eq.lhs)
        @check length(vars) == 1 "Diff. eq $eq has more than on variable in lhs!"
        return (:diffeq, vars[1])
    elseif eq.lhs isa Sym || eq.lhs isa Term
        vars = get_variables(eq.lhs)
        @check length(vars) == 1 "Algebraic eq $eq has more than on variable in lhs!"
        if vars[1] ∈ Set(get_variables(eq.rhs))
            return (:implicit_algebraic, vars[1])
        else
            return (:explicit_algebraic, vars[1])
        end
    elseif isequal(eq.lhs, 0)
        return (:implicit_algebraic, nothing)
    else
        throw(ArgumentError("Uknnown equation type $eq"))
    end
end

"""
    lhs_var(eq::Equation)

Returns the variable on the lhs of the equation for equations.
"""
lhs_var(eq::Equation) = eq_type(eq)[2]

"""
    recusive_substitute(term, rules::Dict)

Apply substitutions until there is no more change in `term`.
"""
function recursive_substitute(term, rules::Dict)
    new_term = substitute(term, rules)
    while !isequal(new_term, term)
        term = new_term
        new_term = substitute(term, rules)
    end
    return new_term
end


text:
function gimbal_plant(;name, J=1.0, k_s=1.0, k_v=1.0, T_c=1.0, ω_brk=1.0, T_brk=1.0)
    @parameters t _J _k_s _k_v _T_c _ω_brk _T_brk
    @variables T(t) ω_rel(t) θ_rel(t) ω_body(t) ω(t)
    D = Differential(t)
    eqs = [
        D(θ_rel) ~ ω_rel
        ω ~ ω_rel + ω_body
        D(ω) ~ (T-_k_s*θ_rel-_k_v*ω_rel
                -_T_c*tanh(10*ω_rel/_ω_brk)
                -sqrt(2ℯ)*(_T_brk-_T_c)*exp(-(ω_rel/(_ω_brk*sqrt(2)))^2)*(ω_rel/(_ω_brk*sqrt(2)))
                )/_J
    ]
    ODESystem(eqs, t; name, defaults =
        [ω=>0.0, ω_body=>0.0, ω_rel=>0.0,
         θ_rel=>0.0, T=>0.0,
         _T_c=>T_c, _T_brk=>T_brk, _ω_brk=>ω_brk,
         _J=>J, _k_s=>k_s, _k_v=>k_v])
end

export gimbal_plant


text:
maxdims["dmrg_2d_qns"] = 1_000:1_000:5_000
descriptions["dmrg_2d_qns"] = "DMRG, 2D Hubbard model (U/t = 4)\n(Nx, Ny) = (8, 4), 10 sweeps\nReal space, periodic in the y-direction"

function runbenchmark(::Val{:dmrg_2d_qns};
                      maxdim::Int, nsweeps::Int = 10,
                      outputlevel::Int = 0, cutoff::Float64 = 0.0,
                      Nx::Int = 8, Ny::Int = 4,
                      U::Float64 = 4.0, t::Float64 = 1.0,
                      splitblocks::Bool = false)
  N = Nx * Ny
  sweeps = Sweeps(nsweeps)
  maxdims = min.(maxdim, [100, 200, 400, 800, 2000, 3000, maxdim])
  maxdim!(sweeps, maxdims...)
  cutoff!(sweeps, cutoff)
  noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  sites = siteinds("Electron", N; conserve_qns = true)
  lattice = square_lattice(Nx, Ny; yperiodic = true)
  ampo = AutoMPO()
  for b in lattice
    ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
    ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
    ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
    ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
  end
  for n in 1:N
    ampo .+= U, "Nupdn", n
  end
  H = MPO(ampo,sites)
  if splitblocks
    H = ITensors.splitblocks(linkinds, H)
  end
  state = [isodd(n) ? "↑" : "↓" for n in 1:N]
  ψ0 = productMPS(sites, state)
  energy, ψ = dmrg(H, ψ0, sweeps;
                   outputlevel = outputlevel)
  if outputlevel > 0
    @show nsweeps
    @show Nx, Ny
    @show U, t
    @show cutoff
    @show energy
    @show flux(ψ)
  end
  return maxlinkdim(ψ)
end


text:
abstract type AbstractGeometricalAxisWeights{T} <: AbstractArray{T, 2} end

function sizeof(gw::AbstractGeometricalAxisWeights{T})::Int where {T}
    return sizeof(gw.weights)
end
@inline function size(gw::AbstractGeometricalAxisWeights{T}) where {T}
    return size(gw.weights)
end

@inline getindex(gw::AbstractGeometricalAxisWeights{T}, I::Vararg{Int, N}) where {T, N} = getindex(gw.weights, I...)

struct GeometricalCartesianAxisWeights{T} <: AbstractGeometricalAxisWeights{T}
    weights::Array{T, 2}
end

struct GeometricalAzimutalAxisWeights{T} <: AbstractGeometricalAxisWeights{T}
    weights::Array{T, 2}
end

struct GeometricalRadialAxisWeights{T} <: AbstractGeometricalAxisWeights{T}
    weights::Array{T, 2}
end

function GeometricalCartesianAxisWeights( ax::DiscreteAxis{T, BL, BR} )::GeometricalCartesianAxisWeights{T} where {T, BL, BR}
    axticks::Vector{T} = collect(ax.ticks)
    ax_ext::Vector{T} = get_extended_ticks(ax)
    Δax_ext::Vector{T} = diff(ax_ext)
    Δax_ext_inv::Vector{T} = inv.(Δax_ext)
    mpax::Vector{T} = midpoints(ax_ext)
    Δmpax::Vector{T} = diff(mpax)
    Δmpax_inv::Vector{T} = inv.(Δmpax)
    Δhmpaxr::Vector{T} = mpax[2:end] - axticks
    Δhmpaxl::Vector{T} = axticks - mpax[1:end - 1]
    waxr::Vector{T} = Δmpax_inv .* Δhmpaxr
    waxl::Vector{T} = Δmpax_inv .* Δhmpaxl

    w::Array{T, 2} = zeros(T, 4, length(waxr) + 1)
    w[1, 1:length(waxr)] = waxr
    w[2, 1:length(waxr)] = waxl
    w[3, 1:length(waxr)] = Δmpax
    w[4, :] = Δax_ext_inv

    return GeometricalCartesianAxisWeights{T}( w )
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module APITests
import ..Pkg # ensure we are using the correct Pkg

using Pkg, Test, REPL
import Pkg.Types.PkgError, Pkg.Resolve.ResolverError
using UUIDs

using ..Utils
@testset "Pkg.activate" begin
    isolate() do; cd_tempdir() do tmp
        path = pwd()
        Pkg.activate(".")
        mkdir("Foo")
        cd(mkdir("modules")) do
            Pkg.generate("Foo")
        end
        Pkg.develop(Pkg.PackageSpec(path="modules/Foo")) # to avoid issue #542
        Pkg.activate("Foo") # activate path Foo over deps Foo
        @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
        Pkg.activate(".")
        rm("Foo"; force=true, recursive=true)
        Pkg.activate("Foo") # activate path from developed Foo
        @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
        Pkg.activate(".")
        Pkg.activate("./Foo") # activate empty directory Foo (sidestep the developed Foo)
        @test Base.active_project() == joinpath(path, "Foo", "Project.toml")
        Pkg.activate(".")
        Pkg.activate("Bar") # activate empty directory Bar
        @test Base.active_project() == joinpath(path, "Bar", "Project.toml")
        Pkg.activate(".")
        Pkg.add("Example") # non-deved deps should not be activated
        Pkg.activate("Example")
        @test Base.active_project() == joinpath(path, "Example", "Project.toml")
        Pkg.activate(".")
        cd(mkdir("tests"))
        Pkg.activate("Foo") # activate developed Foo from another directory
        @test Base.active_project() == joinpath(path, "modules", "Foo", "Project.toml")
        Pkg.activate() # activate home project
        @test Base.ACTIVE_PROJECT[] === nothing
    end end
end

include("FakeTerminals.jl")
import .FakeTerminals.FakeTerminal

@testset "Pkg.precompile_script" begin
    function fake_repl(@nospecialize(f); options::REPL.Options=REPL.Options(confirm_exit=false))
        # Use pipes so we can easily do blocking reads
        # In the future if we want we can add a test that the right object
        # gets displayed by intercepting the display
        input = Pipe()
        output = Pipe()
        err = Pipe()
        Base.link_pipe!(input, reader_supports_async=true, writer_supports_async=true)
        Base.link_pipe!(output, reader_supports_async=true, writer_supports_async=true)
        Base.link_pipe!(err, reader_supports_async=true, writer_supports_async=true)

        repl = REPL.LineEditREPL(FakeTerminal(input.out, output.in, err.in), true)
        repl.options = options

        f(input.in, output.out, repl)
        t = @async begin
            close(input.in)
            close(output.in)
            close(err.in)
        end
        @test read(err.out, String) == ""
        #display(read(output.out, String))
        Base.wait(t)
        nothing
    end
    pwd_before = pwd()
    fake_repl() do stdin_write, stdout_read, repl
        repltask = @async REPL.run_repl(repl)

        for line in split(Pkg.precompile_script, "\n"; keepempty=false)
            sleep(0.1)
            # Consume any extra output
            if bytesavailable(stdout_read) > 0
                copyback = readavailable(stdout_read)
                #@info(copyback)
            end

            # Write the line
            write(stdin_write, line, "\n")

            # Read until some kind of prompt
            readuntil(stdout_read, "\n")
            readuntil(stdout_read, ">")
            #@info(line)
        end

        write(stdin_write, "\x04")
        wait(repltask)
    end
    cd(pwd_before) # something in the precompile_script changes the working directory
end

@testset "Pkg.precompile" begin
    # sequential precompile, depth-first
    isolate() do; cd_tempdir() do tmp
        Pkg.activate(".")
        cd(mkdir("packages")) do
            Pkg.generate("Dep1")
            Pkg.generate("Dep2")
            Pkg.generate("Dep3")
            Pkg.generate("Dep4")
            Pkg.generate("Dep5")
            Pkg.generate("NoVersion")
            open(joinpath("NoVersion","Project.toml"), "w") do io
                write(io, "name = \"NoVersion\"\nuuid = \"$(UUIDs.uuid4())\"")
            end
            Pkg.generate("BrokenDep")
            open(joinpath("BrokenDep","src","BrokenDep.jl"), "w") do io
                write(io, "module BrokenDep\nerror()\nend")
            end
        end
        Pkg.develop(Pkg.PackageSpec(path="packages/Dep1"))

        Pkg.activate("Dep1")
        Pkg.develop(Pkg.PackageSpec(path="packages/Dep2"))
        Pkg.activate("Dep2")
        Pkg.develop(Pkg.PackageSpec(path="packages/Dep3"))

        Pkg.activate(".")
        Pkg.resolve()
        Pkg.precompile()

        iob = IOBuffer()
        ENV["JULIA_PKG_PRECOMPILE_AUTO"]=1
        println("Auto precompilation enabled")
        Pkg.develop(Pkg.PackageSpec(path="packages/Dep4"))
        Pkg.develop(Pkg.PackageSpec(path="packages/NoVersion")) # a package with no version number
        Pkg.build(io=iob) # should trigger auto-precomp
        @test occursin("Precompiling", String(take!(iob)))
        Pkg.precompile(io=iob)
        @test !occursin("Precompiling", String(take!(iob))) # test that the previous precompile was a no-op

        Pkg.precompile("Dep4", io=iob)
        @test !occursin("Precompiling", String(take!(iob))) # should be a no-op
        Pkg.precompile(["Dep4", "NoVersion"], io=iob)
        @test !occursin("Precompiling", String(take!(iob))) # should be a no-op

        ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0
        println("Auto precompilation disabled")
        Pkg.develop(Pkg.PackageSpec(path="packages/Dep5"))
        Pkg.precompile(io=iob)
        @test occursin("Precompiling", String(take!(iob)))
        @test isempty(Pkg.API.pkgs_precompile_suspended)

        ENV["JULIA_PKG_PRECOMPILE_AUTO"]=1
        Pkg.develop(Pkg.PackageSpec(path="packages/BrokenDep"))
        Pkg.build(io=iob) # should trigger auto-precomp and soft-error
        @test occursin("Precompiling", String(take!(iob)))
        broken_packages = Pkg.API.pkgs_precompile_suspended
        @test length(broken_packages) == 1
        Pkg.activate("newpath")
        Pkg.precompile(io=iob)
        @test !occursin("Precompiling", String(take!(iob))) # test that the previous precompile was a no-op
        @test isempty(Pkg.API.pkgs_precompile_suspended)

        Pkg.activate(".") # test that going back to the project restores suspension list
        Pkg.update("BrokenDep", io=iob) # should trigger auto-precomp but do nothing due to error suspension
        @test !occursin("Precompiling", String(take!(iob)))
        @test length(Pkg.API.pkgs_precompile_suspended) == 1

        @test_throws PkgError Pkg.precompile() # calling precompile should retry any suspended, and throw on errors
        @test Pkg.API.pkgs_precompile_suspended == broken_packages

        ptoml = joinpath("packages","BrokenDep","Project.toml")
        lines = readlines(ptoml)
        open(joinpath("packages","BrokenDep","src","BrokenDep.jl"), "w") do io
            write(io, "module BrokenDep\n\nend") # remove error
        end
        open(ptoml, "w") do io
            for line in lines
                if startswith(line, "version = \"0.1.0\"")
                    println(io, replace(line, "version = \"0.1.0\"" => "version = \"0.1.1\"", count=1)) # up version
                else
                    println(io, line)
                end
            end
        end
        Pkg.update("BrokenDep") # should trigger auto-precomp including the fixed BrokenDep
        Pkg.precompile(io=iob)
        @test !occursin("Precompiling", String(take!(iob))) # test that the previous precompile was a no-op

        # https://github.com/JuliaLang/Pkg.jl/pull/2142
        Pkg.build(; verbose=true)

        ENV["JULIA_PKG_PRECOMPILE_AUTO"]=0
    end end
    # ignoring circular deps, to avoid deadlock
    isolate() do; cd_tempdir() do tmp
        Pkg.activate(".")
        cd(mkdir("packages")) do
            Pkg.generate("CircularDep1")
            Pkg.generate("CircularDep2")
            Pkg.generate("CircularDep3")
        end
        Pkg.develop(Pkg.PackageSpec(path="packages/CircularDep1"))
        Pkg.develop(Pkg.PackageSpec(path="packages/CircularDep2"))
        Pkg.develop(Pkg.PackageSpec(path="packages/CircularDep3"))

        Pkg.activate("CircularDep1")
        Pkg.develop(Pkg.PackageSpec(path="packages/CircularDep2"))
        Pkg.activate("CircularDep2")
        Pkg.develop(Pkg.PackageSpec(path="packages/CircularDep3"))
        Pkg.activate("CircularDep3")
        Pkg.develop(Pkg.PackageSpec(path="packages/CircularDep1"))

        Pkg.activate(".")
        Pkg.resolve()

        ## Tests when circularity is in dependencies
        @test_logs (:warn, r"Circular dependency detected") Pkg.precompile()

        ## Tests when circularity goes through the active project
        Pkg.activate("CircularDep1")
        @test_logs (:warn, r"Circular dependency detected") Pkg.precompile()
        Pkg.activate(".")
        Pkg.activate("CircularDep2")
        @test_logs (:warn, r"Circular dependency detected") Pkg.precompile()
        Pkg.activate(".")
        Pkg.activate("CircularDep3")
        @test_logs (:warn, r"Circular dependency detected") Pkg.precompile()
    end end
end

@testset "Pkg.API.check_package_name: Error message if package name ends in .jl" begin
    @test_throws Pkg.Types.PkgError("`Example.jl` is not a valid package name. Perhaps you meant `Example`") Pkg.API.check_package_name("Example.jl")
end

@testset "issue #2587, PR #2589: `Pkg.PackageSpec` accepts `Union{UUID, AbstractString, Nothing}` for `uuid`" begin
    @testset begin
        xs = [
            Pkg.PackageSpec(uuid = Base.UUID(0)),
            Pkg.PackageSpec(uuid = Base.UUID("00000000-0000-0000-0000-000000000000")),
            Pkg.PackageSpec(uuid = "00000000-0000-0000-0000-000000000000"),
            Pkg.PackageSpec(uuid = strip("00000000-0000-0000-0000-000000000000")), # `strip` returns a `SubString{String}`, which is an `AbstractString` but is not a `String`
        ]
        for x in xs
            @test x isa Pkg.PackageSpec
            @test x.uuid isa Base.UUID
            @test x.uuid == Base.UUID(0)
        end
    end

    @testset begin
        xs = [
            Pkg.PackageSpec(),
            Pkg.PackageSpec(uuid = nothing),
        ]
        for x in xs
            @test x isa Pkg.PackageSpec
            @test x.uuid === nothing
        end
    end
end

end # module APITests


text:
"""
Functionality for dealing with HTTP cookies.
"""
module Cookies

import Revise
import HTTP
import Genie, Genie.Encryption, Genie.HTTPUtils

"""
    get(payload::Union{HTTP.Response,HTTP.Request}, key::Union{String,Symbol}, default::T; encrypted::Bool = true)::T where T

Attempts to get the Cookie value stored at `key` within `payload`.
If the `key` is not set, the `default` value is returned.

# Arguments
- `payload::Union{HTTP.Response,HTTP.Request}`: the request or response object containing the Cookie headers
- `key::Union{String,Symbol}`: the name of the cookie value
- `default::T`: default value to be returned if no cookie value is set at `key`
- `encrypted::Bool`: if `true` the value stored on the cookie is automatically decrypted
"""
function get(payload::Union{HTTP.Response,HTTP.Request}, key::Union{String,Symbol}, default::T; encrypted::Bool = true)::T where T
  val = get(payload, key, encrypted = encrypted)
  val === nothing ? default : parse(T, val)
end

"""
    get(res::HTTP.Response, key::Union{String,Symbol}) :: Union{Nothing,String}

Retrieves a value stored on the cookie as `key` from the `Respose` object.

# Arguments
- `payload::Union{HTTP.Response,HTTP.Request}`: the request or response object containing the Cookie headers
- `key::Union{String,Symbol}`: the name of the cookie value
- `encrypted::Bool`: if `true` the value stored on the cookie is automatically decrypted
"""
function get(res::HTTP.Response, key::Union{String,Symbol}; encrypted::Bool = true) :: Union{Nothing,String}
  (haskey(HTTPUtils.Dict(res), "Set-Cookie") || haskey(HTTPUtils.Dict(res), "set-cookie")) ?
    nullablevalue(res, key, encrypted = encrypted) :
      nothing
end

"""
    get(req::Request, key::Union{String,Symbol}) :: Union{Nothing,String}

Retrieves a value stored on the cookie as `key` from the `Request` object.

# Arguments
- `req::HTTP.Request`: the request or response object containing the Cookie headers
- `key::Union{String,Symbol}`: the name of the cookie value
- `encrypted::Bool`: if `true` the value stored on the cookie is automatically decrypted
"""
function get(req::HTTP.Request, key::Union{String,Symbol}; encrypted::Bool = true) :: Union{Nothing,String}
  (haskey(HTTPUtils.Dict(req), "cookie") || haskey(HTTPUtils.Dict(req), "Cookie")) ?
    nullablevalue(req, key, encrypted = encrypted) :
      nothing
end

"""
    set!(res::HTTP.Response, key::Union{String,Symbol}, value::Any, attributes::Dict; encrypted::Bool = true) :: HTTP.Response

Sets `value` under the `key` label on the `Cookie`.

# Arguments
- `res::HTTP.Response`: the HTTP.Response object
- `key::Union{String,Symbol}`: the key for storing the cookie value
- `value::Any`: the cookie value
- `attributes::Dict`: additional cookie attributes, such as `Path` or `HttpOnly`
- `encrypted::Bool`: if `true` the value is stored encoded
"""
function set!(res::HTTP.Response, key::Union{String,Symbol}, value::Any, attributes::Dict = Dict(); encrypted::Bool = true) :: HTTP.Response
  normalized_attrs = Dict{Symbol,Any}()
  for (k,v) in attributes
    normalized_attrs[Symbol(lowercase(string(k)))] = v
  end

  value = string(value)
  encrypted && (value = Genie.Encryption.encrypt(value))
  cookie = HTTP.Cookies.Cookie(string(key), value; normalized_attrs...)

  headers = Dict(res.headers)
  if haskey(headers, "Set-Cookie")
    headers["Set-Cookie"] *= "\nSet-Cookie: " * HTTP.Cookies.String(cookie, false) * "; "
  elseif haskey(headers, "set-cookie")
    headers["set-cookie"] *= "\nset-cookie: " * HTTP.Cookies.String(cookie, false) * "; "
  else
    headers["Set-Cookie"] = HTTP.Cookies.String(cookie, false) * "; "
  end
  res.headers = [(k => v) for (k,v) in headers]

  res
end

"""
    Dict(req::Request) :: Dict{String,String}

Extracts the `Cookie` and `Set-Cookie` data from the `Request` and `Response` objects and converts it into a Dict.
"""
function Base.Dict(r::Union{HTTP.Request,HTTP.Response}) :: Dict{String,String}
  d = Dict{String,String}()
  headers = Dict(r.headers)

  h = if haskey(headers, "Cookie")
    split(headers["Cookie"], ";")
  elseif haskey(headers, "cookie")
    split(headers["cookie"], ";")
  elseif haskey(headers, "Set-Cookie")
    split(headers["Set-Cookie"], ";")
  elseif haskey(headers, "set-cookie")
    split(headers["set-cookie"], ";")
  else
    []
  end

  for cookie in h
    cookie_parts = split(cookie, "=")
    if length(cookie_parts) == 2
      d[strip(cookie_parts[1])] = cookie_parts[2]
    else
      d[strip(cookie_parts[1])] = ""
    end
  end

  d
end

### PRIVATE ###

"""
    nullablevalue(payload::Union{HTTP.Response,HTTP.Request}, key::Union{String,Symbol}; encrypted::Bool = true)

Attempts to retrieve a cookie value stored at `key` in the `payload object` and returns a `Union{Nothing,String}`

# Arguments
- `payload::Union{HTTP.Response,HTTP.Request}`: the request or response object containing the Cookie headers
- `key::Union{String,Symbol}`: the name of the cookie value
- `encrypted::Bool`: if `true` the value stored on the cookie is automatically decrypted
"""
function nullablevalue(payload::Union{HTTP.Response,HTTP.Request}, key::Union{String,Symbol}; encrypted::Bool = true) :: Union{Nothing,String}
  for cookie in split(Dict(payload)["cookie"], ';')
    cookie = strip(cookie)
    if startswith(lowercase(cookie), lowercase(string(key)))
      value = split(cookie, '=')[2] |> String
      encrypted && (value = Genie.Encryption.decrypt(value))

      return string(value)
    end
  end

  nothing
end

"""
    getcookies(req::HTTP.Request) :: Vector{HTTP.Cookies.Cookie}

Extracts cookies from within `req`
"""
function getcookies(req::HTTP.Request) :: Vector{HTTP.Cookies.Cookie}
  HTTP.Cookies.cookies(req)
end

"""
    getcookies(req::HTTP.Request) :: Vector{HTTP.Cookies.Cookie}

Extracts cookies from within `req`, filtering them by `matching` name.
"""
function getcookies(req::HTTP.Request, matching::String) :: Vector{HTTP.Cookies.Cookie}
  HTTP.Cookies.readcookies(req.headers, matching)
end

end


text:
###
### LongSequence specific specializations of src/biosequence/transformations.jl
###

"""
    resize!(seq, size)

Resize a biological sequence `seq`, to a given `size`.
"""
function Base.resize!(seq::LongSequence{A}, size::Integer) where {A}
    if size < 0
        throw(ArgumentError("size must be non-negative"))
    end
    orphan!(seq, size)
    resize!(seq.data, seq_data_len(A, size + seq.part.start - 1))
    seq.part = seq.part.start:seq.part.start+size-1
    return seq
end

function Base.filter!(f::Function, seq::LongSequence{A}) where {A}
    orphan!(seq)

    len = 0
    next = bitindex(seq, 1)
    j = index(next)
    datum::UInt64 = 0
    for i in 1:lastindex(seq)
        x = inbounds_getindex(seq, i)
        if f(x)
            datum |= enc64(seq, x) << offset(next)
            len += 1
            #TODO: Resolve use of bits_per_symbol.
            next += bits_per_symbol(A())
            if index(next) != j
                seq.data[j] = datum
                datum = 0
                j = index(next)
            end
        end
    end
    if offset(next) > 0
        seq.data[j] = datum
    end
    resize!(seq, len)

    return seq
end

function Base.map!(f::Function, seq::LongSequence)
    orphan!(seq)
    for i in 1:lastindex(seq)
        unsafe_setindex!(seq, f(inbounds_getindex(seq, i)), i)
    end
    return seq
end

"""
    reverse!(seq::LongSequence)

Reverse a biological sequence `seq` in place.
"""
function Base.reverse!(seq::LongSequence)
    orphan!(seq) # TODO: Is the orphan call really nessecery given the indexing calls will also call orphan?
    @inbounds for i in 1:div(lastindex(seq), 2)
	    x = seq[i]
        i′ = lastindex(seq) - i + 1
	    seq[i] = seq[i′]
	    seq[i′] = x
    end
    return seq
end

"""
    Base.reverse(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet}

Create a reversed copy of a LongSequence representing a DNA or RNA sequence.

This version of the reverse method is optimized for speed, taking advantage of
bit-parallel operations.
"""
function Base.reverse(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet}
    data = Vector{UInt64}(undef, seq_data_len(A, length(seq)))
    i = 1
    next = lastbitindex(seq)
    stop = bitindex(seq, 0)
    r = rem(offset(next) + bits_per_symbol(seq), 64)
    if r == 0
        @inbounds while next - stop > 0
            x = seq.data[index(next)]
            data[i] = reversebits(x, BitsPerSymbol(seq))
            i += 1
            next -= 64
        end
    else
        @inbounds while next - stop > 64
            j = index(next)
            x = (seq.data[j] << (64 - r)) | (seq.data[j - 1] >> r)
            data[i] = reversebits(x, BitsPerSymbol(seq))
            i += 1
            next -= 64
        end
        if next - stop > 0
            j = index(next)
            x = seq.data[j] << (64 - r)
            if r < next - stop
                x |= seq.data[j - 1] >> r
            end
            data[i] = reversebits(x, BitsPerSymbol(seq))
        end
    end
    return LongSequence{A}(data, 1:length(seq), false)
end

"""
    complement!(seq)

Make a complement sequence of `seq` in place.
"""
function complement!(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet}
    orphan!(seq)
    next = firstbitindex(seq)
    stop = bitindex(seq, lastindex(seq) + 1)
    seqdata = seq.data
    @inbounds while next < stop
        x = seqdata[index(next)]
        seqdata[index(next)] = complement_bitpar(x, Alphabet(seq))
        next += 64
    end
    return seq
end

###
### Shuffle
###

function Random.shuffle!(seq::LongSequence)
    orphan!(seq) # TODO: Is this call to orphan nessecery, given setindex should call `orphan!` for us?
    # Fisher-Yates shuffle
    for i in 1:lastindex(seq) - 1
        j = rand(i:lastindex(seq))
        seq[i], seq[j] = seq[j], seq[i]
    end
    return seq
end


text:
using Base.Test
using IntrinsicPersistence
using Distances
using LightGraphs, SimpleWeightedGraphs
using StaticArrays

macro test_inferred(ex)
    :(@test @inferred($(esc(ex))) ≠ nothing)
end

include("TriangleIterator.jl")
include("persistence.jl")


text:
"""
Test script for loopless constraints in MinFlux.
Model construction and package imports can be found in `testmodel_construction.jl`.
"""

include("testmodel_construction.jl")

minflux_bounds = MinFlux_preprocessing(m, reaction_dG)
substrates = ["m1"]
targets = ["m6"]
cofactors = ["co1", "co2", "h"] #allow cofactors

bounds = Dict("substrates" => -1, "reactants" => (-15, 15))

os = build_OptStoic_model(m, substrates, targets, cofactors, Gurobi, energy_dc = formation_dG, variable_bounds = bounds)
os = optimize_OptStoic(os)
os_solution = solution_OptStoic(os)

if !isnothing(os_solution)
    m, mf = build_MinFlux_model(m, os_solution, minflux_bounds, optimizer = Gurobi, loopless = true)
    mf = MinFlux(mf)
    MinFlux_solution(m, mf)
end #Output should be the same as for `optimization.jl`

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
"""
    SeisPlotCoordinates(headers;<keyword arguments>)

Plot shot and receiver coordinates or fold map. Optionally can print to a file.

# Arguments
* `headers` :a vector where each element is of type Header

# Keyword arguments
* `style="sxsygxgy"` or "fold"
* `cmap="Greys"` or "PuOr"
* `aspect="auto"`: color image aspect ratio.
* `pclip=98`: percentile for determining clip.
* `vmin="NULL"`: minimum value used in colormapping data.
* `vmax="NULL"`: maximum value used in colormapping data.
* `title=" "`: title of plot.
* `xlabel=" "`: label on x-axis.
* `xunits=" "`: units of y-axis.
* `ylabel=" "`: label on x-axis.
* `yunits=" "`: units of y-axis.
* `ox=0`: first point of x-axis.
* `dx=1`: increment of x-axis.
* `oy=0`: first point of y-axis.
* `dy=1`: increment of y-axis.
* `dpi=100`: dots-per-inch of figure.
* `wbox=6`: width of figure in inches.
* `hbox=6`: height of figure in inches.
* `name="NULL"`: name of the figure to save (only if `name` is given).
* `interpolation="none"`: interpolation method for colormapping data.
* `titlesize=16`: size of title.
* `labelsize=14`: size of labels on axis.
* `ticksize=11`: size of labels on ticks.
* `fignum="NULL"`: number of figure.

*Credits: Aaron Stanton, 2015*

"""
function SeisPlotCoordinates(headers;style="sxsygxgy",cmap="Greys",aspect="auto",pclip=98,vmin="NULL",vmax="NULL",title=" ",xlabel=" ",xunits=" ",ylabel=" ",yunits=" ",ox=0,dx=1,oy=0,dy=1,dpi=100,wbox=6,hbox=6,name="NULL",interpolation="none",titlesize=20,labelsize=15,ticksize=10,fignum="NULL")

	plt.ion()
	if (fignum == "NULL")
		fig = plt.figure(figsize=(wbox, hbox), dpi=dpi, facecolor="w", edgecolor="k")
	else
		fig = plt.figure(num=fignum,figsize=(wbox, hbox), dpi=dpi, facecolor="w", edgecolor="k")
	end
	if (style == "sxsygxgy")
		sx = Seismic.ExtractHeader(headers,"sx")
		sy = Seismic.ExtractHeader(headers,"sy")
		gx = Seismic.ExtractHeader(headers,"gx")
		gy = Seismic.ExtractHeader(headers,"gy")
		im = plt.scatter(gx,gy,marker="^",s=5,color="b");
		im = plt.scatter(sx,sy,marker="*",s=8,color="r");
		#plt.axis([ox,ox + (size(in,2)-1)*dx,oy + (size(in,1)-1)*dy,oy])
	elseif (style == "fold")
		imx = Seismic.ExtractHeader(headers,"imx")
		imy = Seismic.ExtractHeader(headers,"imy")
		imx_min = minimum(imx)
		imx_max = maximum(imx)
		nmx = imx_max - imx_min + 1
		imy_min = minimum(imy)
		imy_max = maximum(imy)
		nmy = imy_max - imy_min + 1
		fold = zeros(nmy,nmx)
		for itrace = 1 : length(imx)
			imx_index = imx[itrace] - imx_min + 1
			imy_index = imy[itrace] - imy_min + 1
			fold[imy_index,imx_index] += 1
		end
		if (vmin=="NULL" || vmax=="NULL")
			if (pclip<=100)
				a = -quantile(abs(fold[:]),(pclip/100))
			else
				a = -quantile(abs(fold[:]),1)*pclip/100
			end
			b = -a
		else
			a = vmin
			b = vmax
		end
		im = plt.imshow(fold,cmap=cmap,vmin=a,vmax=b,extent=[ox,ox + (size(fold,2)-1)*dx,oy + (size(fold,1)-1)*dy,oy],aspect=aspect)
		ax = plt.gca()
		ax[:invert_yaxis]();
		plt.colorbar()
	else
		println("no other plotting style defined yet")
	end

	plt.title(title)
	plt.xlabel(join([xlabel " " xunits]))
	plt.ylabel(join([ylabel " " yunits]))

	if (name == "NULL")
		plt.show()
	else
		plt.savefig(name,dpi=dpi)
		plt.close()
	end
	return im
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
using JSObjectLiteral
using Test

## @js tests

e = π
g = "gee!"

a = @js {
  a: 1,
  b: [2, 3 * 3],
  c : {
    d: "doubly-quoted string",
    e
  },
  f: g
}

b = Dict{String, Any}(
  "a" => 1,
  "b" => [2, 9],
  "c" => Dict{String, Any}(
    "d" => "doubly-quoted string",
    "e" => π
  ),
  "f" => "gee!"
)

@test a == b

@test Dict{String, Any}( "a" => "a") == @js { a: "a" }

a = @js { b: { c: { d: 4 } } }

@test a == @js { b.c.d: 4 }
@test @js(a.b.c.d) == 4
@test @js(a.b.c) == Dict("d" => 4)
@test @js(a.b) == @js { c: { d: 4 } }

a = @js { f: x -> x^2, b: { g: √ } }

@test a["f"](5) == 25
@test @js(a.b.g)(49) == 7

@test Dict("a" => "b") == @js { a: "b" }

b = @js { c: 3 }
@test Dict("a" => b) == @js { a: b }

a = 101
@test Dict("a" => a) == @js { a }
@test Dict("a" => a+1 ) == @js { a: a + 1 }

e = exp(1)
@test Dict("a" => a, "b" => Dict("c" => Dict("d" => 5), "e" => e), "f" => [1, Dict("g" => 2)]) == @js { a, b: { c: { d: 5 }, e }, f: [1, { g: 2} ] }
@test Dict("a" => Dict("b" => 3)) == @js { a.b: 3 }
@test [a, b, 3] == @js [ a, b, 3 ]

c = [1, 2, @js({d: 3})]
a = @js { c }
@test @js(a.c[3].d) == 3
@test @js(a.c[a.c[3].d].d) == 3

@js a = { b.c.d: 4 }
@test @js(a.b) == @js { c: { d: 4 } }
@test @js(a.b.c) == @js { d: 4 }

@js a = { b: "c" }
@test a == @js { b: "c" }

@js a.b = "d"
@test a == @js { b: "d" }

@js a.b = { c: "d" }
@test a == @js { b.c: "d" }

@js a.b.c = "e"
@test a == @js { b: { c: "e" } }

@js b = { c: 3 }
@js a.b = b.c
@test a == @js { b: 3 }

@js a, b = [ { a: 1 }, { b: 2} ]
@test a == @js { a: 1 }
@test b == @js { b: 2 }

@js dict = { d: π, e: [ 1, { f: 2}, 3], c: sin }
@js { c, d, e } = dict
@test c == sin
@test d == π
@test e == [ 1, Dict("f" => 2), 3]

## issue #2
@test Dict("Content-Type" => "application/json") == @js { "Content-Type": "application/json"}

## JSObject tests
a = JSObject(@js({ b: { c: { d: 4 } } }))
@test a.b.c.d == 4
@test a["b"].c["d"] == 4

a.c = JSObject(@js({ d: { e: [5, 6] } }))
@test a.c.d.e[2] == 6

@test stripobject(a) == @js { b: { c: { d: 4 } }, c: { d: { e: [5, 6] } } }

b = JSObject(@js([ { c: 3 }, { d: [ 4, { e: 5 } ] } ]))
@test b[1].c == 3
@test b[2].d[1] == 4
@test b[2].d[2].e == 5

b[1] = @js([6, { f: 7}])
@test b[1][1] == 6
@test b[1][2].f == 7

@test stripobject(b) == @js [ [6, { f: 7} ], { d: [4, { e: 5 } ] } ]


text:
using BenchmarkTools
using LoopVectorization
using Flux
using ExplainableAI
import ExplainableAI: modify_layer, lrp!

on_CI = haskey(ENV, "GITHUB_ACTIONS")

include("../test/vgg11.jl")
vgg11 = VGG11(; pretrain=false)
model = flatten_model(strip_softmax(vgg11.layers))

T = Float32
img = rand(MersenneTwister(123), T, (224, 224, 3, 1))

# Benchmark custom LRP composite
function LRPCustom(model::Chain)
    return LRP(
        model,
        [ZBoxRule(zero(T), oneunit(T)), repeat([GammaRule()], length(model.layers) - 1)...],
    )
end

# Use one representative algorithm of each type
algs = Dict(
    "Gradient" => Gradient,
    "InputTimesGradient" => InputTimesGradient,
    "LRPZero" => LRPZero,
    "LRPCustom" => LRPCustom, #modifies weights
    "SmoothGrad" => model -> SmoothGrad(model, 10),
    "IntegratedGradients" => model -> IntegratedGradients(model, 10),
)

# Define benchmark
contruct_analyzer(alg, model) = alg(model) # for use with @benchmarkable macro

SUITE = BenchmarkGroup()
SUITE["VGG"] = BenchmarkGroup([k for k in keys(algs)])
for (name, alg) in algs
    SUITE["VGG"][name] = BenchmarkGroup(["construct analyzer", "analyze"])
    SUITE["VGG"][name]["construct analyzer"] = @benchmarkable contruct_analyzer(
        $(alg), $(model)
    )

    analyzer = alg(model)
    SUITE["VGG"][name]["analyze"] = @benchmarkable analyze($(img), $(analyzer))
end

# Rules benchmarks – use wrapper to trigger AD fallback
struct TestWrapper{T}
    layer::T
end
(w::TestWrapper)(x) = w.layer(x)
modify_layer(r::AbstractLRPRule, w::TestWrapper) = modify_layer(r, w.layer)
lrp!(Rₖ, rule::ZBoxRule, w::TestWrapper, aₖ, Rₖ₊₁) = lrp!(Rₖ, rule, w.layer, aₖ, Rₖ₊₁)

# generate input for conv layers
insize = (64, 64, 3, 1)
in_dense = 500
out_dense = 100
aₖ = randn(T, insize)

layers = Dict(
    "MaxPool" => (MaxPool((3, 3); pad=0), aₖ),
    "Conv" => (Conv((3, 3), 3 => 2), aₖ),
    "Dense" => (Dense(in_dense, out_dense, relu), randn(T, in_dense, 1)),
    "WrappedDense" =>
        (TestWrapper(Dense(in_dense, out_dense, relu)), randn(T, in_dense, 1)),
)
rules = Dict(
    "ZeroRule" => ZeroRule(),
    "EpsilonRule" => EpsilonRule(),
    "GammaRule" => GammaRule(),
    "ZBoxRule" => ZBoxRule(zero(T), oneunit(T)),
)

SUITE["Layer"] = BenchmarkGroup([k for k in keys(layers)])
for (layername, (layer, aₖ)) in layers
    SUITE["Layer"][layername] = BenchmarkGroup([k for k in keys(rules)])
    Rₖ = similar(aₖ)
    Rₖ₊₁ = layer(aₖ)
    for (rulename, rule) in rules
        SUITE["Layer"][layername][rulename] = @benchmarkable lrp!(
            $(Rₖ), $(rule), $(layer), $(aₖ), $(Rₖ₊₁)
        )
    end
end