text:
export doubledescription
# Naive implementation of the double description Algorithm
# See JuliaPolyhedra/ConvexHull.jl for an efficient implementation

polytypefor(::Type{T}) where {T <: Integer} = Rational{BigInt}
polytypefor(::Type{Float32}) = Float64
polytypefor(::Type{T}) where {T} = T

polyvectortype(a) = a
# TODO sparse halfspaces does not mean sparse points
polyvectortype(::Type{<:SparseVector{T}}) where T = Vector{T}

dualtype(RepT::Type{<:Representation}) = dualtype(RepT, polyvectortype(vectortype(RepT)))
function dualfullspace(rep::Representation, d::FullDim, ::Type{T}) where T
    dualfullspace(rep, d, T, polyvectortype(similar_type(vectortype(rep), d, T)))
end
function dualfullspace(rep::Representation{T}) where T
    dualfullspace(rep, FullDim(rep), polytypefor(T))
end

"""
    doubledescription(h::HRepresentation)

Computes the V-representation of the polyhedron represented by `h` using the Double-Description algorithm [MRTT53, FP96].
It maintains a list of the points, rays and lines that are in the resulting representation
and in the hyperplane corresponding to each halfspace and then add their intersection for each halfspace
with [`Polyhedra.add_intersection!`](@ref).
The resulting V-representation has no redundancy.

    doubledescription(V::VRepresentation)

Computes the H-representation of the polyhedron represented by `v` using the Double-Description algorithm [MRTT53, FP96].
Currently, this fallbacks to lifting it to the V-representation
of a cone by homogenization, interpreting it
as the H-representation of the dual cone so that
it can rely on `doubledescription(::HRepresentation)`.

[MRTT53] Motzkin, T. S., Raiffa, H., Thompson, G. L. and Thrall, R. M.
The double description method
*Contribution to the Theory of Games*, *Princeton University Press*, **1953**

[FP96] Fukuda, K. and Prodon, A.
Double description method revisited
*Combinatorics and computer science*, *Springer*, **1996**, 91-111
"""
function doubledescription end

function print_v_summary(v::VRep)
    print("$(npoints(v)) points, $(nrays(v)) rays and $(nlines(v)) lines")
end

function intersect_and_remove_redundancy(v, hs, h; verbose=0)
    if eltype(hs) <: HalfSpace
        str = "halfspace"
    else
        str = "hyperplane"
    end
    for (i, hel) in enumerate(hs)
        if verbose >= 1
            println("Intersecting $str $i/$(length(hs))")
        end
        v_int = v ∩ hel
        if verbose >= 3
            print("Removing duplicates: ")
            print_v_summary(v_int)
            println(".")
        end
        # removeduplicates is cheaper than removevredundancy since the latter
        # needs to go through all the hrep element
        # FIXME not sure what to do here but it will be revamped by
        #       https://github.com/JuliaPolyhedra/Polyhedra.jl/pull/195 anyway
        v_dup = removeduplicates(v_int, OppositeMockOptimizer)
        if verbose >= 3
            print("Removing redundancy: ")
            print_v_summary(v_dup)
            println(".")
        end
        v = removevredundancy(v_dup, h)
        if verbose >= 2
            print("After intersection:  ")
            print_v_summary(v)
            println(".")
        end
    end
    return v
end

function slow_doubledescription(h::HRepresentation, solver = nothing; kws...)
    # The redundancy of V-elements are measured using
    # the number of hyperplane they are in. If there are
    # redundant hyperplanes, it does not matter since no point
    # will be inside them but if there are duplicates it is a problem

    # FIXME Note that removevredundancy, uses `h` which contains all hyperplanes and halfspaces
    # already taken into account but also all the other ones. We should check that this
    # is the right idea.

    # FIXME not sure what to do here but it will be revamped by
    #       https://github.com/JuliaPolyhedra/Polyhedra.jl/pull/195 anyway
    h = removeduplicates(h, solver === nothing ? OppositeMockOptimizer : solver)
    v = dualfullspace(h)
    checkvconsistency(v)
    v = intersect_and_remove_redundancy(v, hyperplanes(h), h; kws...)
    v = intersect_and_remove_redundancy(v, halfspaces(h), h; kws...)
    return v
end

struct CutoffPointIndex
    cutoff::Int
    index::Int
end
Base.show(io::IO, p::CutoffPointIndex) = print(io, "p[$(p.cutoff), $(p.index)]")
struct CutoffRayIndex
    cutoff::Int
    index::Int
end
Base.show(io::IO, r::CutoffRayIndex) = print(io, "r[$(r.cutoff), $(r.index)]")
struct DoubleDescriptionData{PointT, RayT, LineT, HST}
    fulldim::Int
    halfspaces::Vector{HST}
    # Elements ordered by first halfspace cutting it off
    points::Vector{PointT}
    pz::Vector{BitSet}
    cutpoints::Vector{Vector{PointT}}
    cutpz::Vector{Vector{BitSet}}
    pin::Vector{Vector{CutoffPointIndex}}
    rays::Vector{RayT}
    rz::Vector{BitSet}
    cutrays::Vector{Vector{RayT}}
    cutrz::Vector{Vector{BitSet}}
    rin::Vector{Vector{CutoffRayIndex}}
    lines::Vector{LineT}
    cutline::Vector{Union{Nothing, LineT}}
    lineray::Vector{Union{Nothing, CutoffRayIndex}}
    nlines::Vector{Int}
end
function Base.show(io::IO, data::DoubleDescriptionData)
    println(io, "DoubleDescriptionData in $(data.fulldim) dimension:")
    println(io, data.points)
    println(io, data.rays)
    println(io, data.lines)
    for i in reverse(eachindex(data.cutpoints))
        println(io, " Halfspace $i: $(data.halfspaces[i]):")
        if !isempty(data.cutpoints[i])
            println(io, "  Cut points:")
            for j in eachindex(data.cutpoints[i])
                println(io, "  $j: ", data.cutpoints[i][j], " zero at: ", data.cutpz[i][j])
            end
        end
        if !isempty(data.pin[i])
            println(io, "  In: ", data.pin[i])
        end
        if !isempty(data.cutrays[i])
            println(io, "  Cut rays:")
            for j in eachindex(data.cutrays[i])
                println(io, "  $j: ", data.cutrays[i][j], " zero at: ", data.cutrz[i][j])
            end
        end
        if !isempty(data.rin[i])
            println(io, "  In: ", data.rin[i])
        end
        if data.cutline[i] !== nothing
            println(io, "  Cut line: ", data.cutline[i])
            if data.lineray[i] !== nothing
                println(io, "  Line ray: ", data.lineray[i])
            end
        end
        if !iszero(data.nlines[i])
            println(io, "  $(data.nlines[i]) uncut lines left")
        end
    end
end

function DoubleDescriptionData{PointT, RayT, LineT}(fulldim::Integer, hyperplanes, halfspaces) where {PointT, RayT, LineT}
    n = length(halfspaces)
    m = length(hyperplanes)
    return DoubleDescriptionData{PointT, RayT, LineT, eltype(halfspaces)}(
        fulldim,
        halfspaces,
        PointT[],
        BitSet[],
        [PointT[] for i in 1:n],
        [BitSet[] for i in 1:n],
        [CutoffPointIndex[] for i in 1:n],
        RayT[],
        BitSet[],
        [RayT[] for i in 1:n],
        [BitSet[] for i in 1:n],
        [CutoffRayIndex[] for i in 1:n],
        LineT[],
        Union{Nothing, LineT}[nothing for i in 1:(m + n)],
        Union{Nothing, CutoffRayIndex}[nothing for i in 1:n],
        zeros(Int, m + n)
    )
end
function tight_halfspace_indices(data::DoubleDescriptionData, p::CutoffPointIndex)
    if iszero(p.cutoff)
        return data.pz[p.index]
    else
        return data.cutpz[p.cutoff][p.index]
    end
end
function tight_halfspace_indices(data::DoubleDescriptionData, r::CutoffRayIndex)
    if iszero(r.cutoff)
        return data.rz[r.index]
    else
        return data.cutrz[r.cutoff][r.index]
    end
end
function Base.getindex(data::DoubleDescriptionData, p::CutoffPointIndex)
    if iszero(p.cutoff)
        return data.points[p.index]
    else
        return data.cutpoints[p.cutoff][p.index]
    end
end
function Base.getindex(data::DoubleDescriptionData, r::CutoffRayIndex)
    if iszero(r.cutoff)
        return data.rays[r.index]
    else
        return data.cutrays[r.cutoff][r.index]
    end
end

function _bitdot_range(b1::BitSet, b2::BitSet, i, n)
    # They share the hyperplance `i` so we could
    # start `count` at `1` and the iteration at `i + 1`
    # but in the case where `i` is redundant, it was
    # removed from the `BitSet` and it should not be counted
    # so we also need to re-check `i`.
    count = 0
    for j in i:n
        if j in b1 && j in b2
            count += 1
        end
    end
    return count
end
_ray_pair(p1, p2) = false
_ray_pair(::CutoffRayIndex, ::CutoffRayIndex) = true
# Necessary condition for adjacency.
# See Proposition 9 (NC1) of [FP96].
function is_adjacent_nc1(data, i, p1, p2)
    rhs = 1 + _ray_pair(p1, p2) + data.nlines[i] + _bitdot_range(
        tight_halfspace_indices(data, p1),
        tight_halfspace_indices(data, p2),
        i, length(data.halfspaces)
    )
    # `length(data.cutline) - length(data.halfspaces)` is the number of
    # hyperplanes (some may be redundant) in which case this condition is less effective.
    # TODO detect linearity and store the dimension instead of `fulldim`.
    dim = data.fulldim - length(data.nlines) + length(data.halfspaces)
    return dim <= rhs
end
function is_adjacency_breaker(data, i, p, p1, p2)
    return p != p1 && p != p2 && all((i + 1):length(data.halfspaces)) do j
        isin(data, j, p) || !(isin(data, j, p1) && isin(data, j, p2))
    end
end
function isadjacent(data, i::Integer, p1, p2)
    return is_adjacent_nc1(data, i, p1, p2) &&
        # According to Proposition 7 (c) of [FP96], we need to check
        # that there is not other point or ray that is in the same hyperplanes as `p1` and `p2`.
        # If it's the case, it is in hyperplane `i` in particular,
        # we can use it to restrict our attention to points and
        # rays in `pin[i]` and `rin[i]`.
        (_ray_pair(p1, p2) || !any(p -> is_adjacency_breaker(data, i, p, p1, p2), data.pin[i])) &&
        !any(r -> is_adjacency_breaker(data, i, r, p1, p2), data.rin[i])
end
isin(data, i, p) = i in tight_halfspace_indices(data, p)

resized_bitset(data) = sizehint!(BitSet(), length(data.halfspaces))
function add_index!(data, cutoff::Nothing, p::AbstractVector, tight::BitSet)
    push!(data.points, p)
    push!(data.pz, tight)
    return CutoffPointIndex(0, length(data.points))
end
function add_index!(data, cutoff::Integer, p::AbstractVector, tight::BitSet)
    push!(data.cutpoints[cutoff], p)
    push!(data.cutpz[cutoff], tight)
    return CutoffPointIndex(cutoff, length(data.cutpoints[cutoff]))
end
function add_index!(data, cutoff::Nothing, r::Polyhedra.Ray, tight::BitSet)
    push!(data.rays, r)
    push!(data.rz, tight)
    return CutoffRayIndex(0, length(data.rays))
end
function add_index!(data, cutoff::Integer, r::Polyhedra.Ray, tight::BitSet)
    push!(data.cutrays[cutoff], r)
    push!(data.cutrz[cutoff], tight)
    return CutoffRayIndex(cutoff, length(data.cutrays[cutoff]))
end

function add_in!(data, i, index::CutoffPointIndex)
    push!(data.pin[i], index)
end
function add_in!(data, i, index::CutoffRayIndex)
    push!(data.rin[i], index)
end
function set_in!(data, I, index)
    for i in I
        if isin(data, i, index)
            add_in!(data, i, index)
        end
    end
end
function add_element!(data, k, el, tight)
    cutoff = nothing
    for i in reverse(1:k)
        if data.cutline[i] !== nothing
            el = line_project(el, data.cutline[i], data.halfspaces[i])
            # The line is in all halfspaces from `i+1` up so projecting with it does not change it.
            push!(tight, i)
            return add_adjacent_element!(data, i, k, el, data.lineray[i], tight)
        end
        # Could avoid computing `dot` twice between `el` and the halfspace here.
        if !(el in data.halfspaces[i])
            cutoff = i
            break
        end
        if el in hyperplane(data.halfspaces[i])
            push!(tight, i)
        end
    end
    index = add_index!(data, cutoff, el, tight)
    set_in!(data, (index.cutoff + 1):k, index)
    return index
end
function project_onto_affspace(data, offset, el, hyperplanes)
    for i in reverse(eachindex(hyperplanes))
        line = data.cutline[offset + i]
        h = hyperplanes[i]
        if line !== nothing
            el = line_project(el, line, h)
        elseif !(el in h)
            # e.g. 0x1 + 0x2 = 1 or -1 = x1 = 1, ...
            return nothing
        end
    end
    return el
end
function add_adjacent_element!(data, i, k, el, parent, tight)
    index = add_element!(data, i - 1, el, tight)
    add_intersection!(data, index, parent, nothing, i - 1)
    set_in!(data, i:k, index)
    return index
end

function combine(β, p1::AbstractVector, value1, p2::AbstractVector, value2)
    λ = (value2 - β) / (value2 - value1)
    return λ * p1 + (1 - λ) * p2
end
function combine(β, p::AbstractVector, pvalue, r::Polyhedra.Ray, rvalue)
    λ = (β - pvalue) / rvalue
    return p + λ * r
end
combine(β, r::Polyhedra.Ray, rvalue, p::AbstractVector, pvalue) = combine(β, p, pvalue, r, rvalue)
function combine(β, r1::Polyhedra.Ray, value1, r2::Polyhedra.Ray, value2)
    # should take
    # λ = value2 / (value2 - value1)
    @assert 0 <= value2 / (value2 - value1) <= 1
    # By homogeneity we can avoid the division and do
    #newr = value2 * r1 - value1 * r2
    # but this can generate very large numbers (see JuliaPolyhedra/Polyhedra.jl#48)
    # so we still divide
    newr = (value2 * r1 - value1 * r2) / (value2 - value1)
    # In CDD, it does value2 * r1 - value1 * r2 but then it normalize the ray
    # by dividing it by its smallest nonzero entry (see dd_CreateNewRay)
    return Polyhedra.simplify(newr)
end

combine(h, el1, el2) = combine(h.β, el1, h.a ⋅ el1, el2, h.a ⋅ el2)

"""
    add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)

Add the intersection of the elements of indices `idx1` and `idx2` that belong to
the hyperplane corresponding the halfspace of index `hp_idx` (see 3.2 (ii) of
[FP96]) or have inherited adjacency if `hp_idx === nothing` (see 3.2 (i) of
[FP96]). In case of inherited adjacency, `hs_idx` is the halfspace where `idx1`
was created.
If `idx1` and `idx2` are not adjacent or if they are in the hyperplane
corresponding to a halfspace of lower index then the intersection is not added
to avoid adding redundant elements.

[FP96] Fukuda, K. and Prodon, A.
Double description method revisited
*Combinatorics and computer science*, *Springer*, **1996**, 91-111
"""
function add_intersection!(data, idx1, idx2, hp_idx, hs_idx = hp_idx - 1)
    if idx1.cutoff > idx2.cutoff
        return add_intersection!(data, idx2, idx1, hp_idx, hs_idx)
    end
    i = idx2.cutoff
       # Condition (c_k) in [FP96]
    if idx1.cutoff == i ||
        isin(data, i, idx1) ||
        # If it's in both at some lower `j` then we'll
        # add it then otherwise, it will be added twice.
        any(j -> isin(data, j, idx1) && isin(data, j, idx2), (i + 1):hs_idx) ||
        (hp_idx !== nothing && !isadjacent(data, hp_idx, idx1, idx2))
        return
    end
    newel = combine(data.halfspaces[i], data[idx1], data[idx2])
    # `newel` and `idx1` have inherited adjacency, see 3.2 (i) of [FP96].
    # TODO are we sure that they are adjacent ?
    #      or should we check rank or combinatorial check ?
    #      In `DDMethodVariation2` of [FP96], `Adj` is not called
    #      in case `rj, rj'` have inherited adjacency
    tight = tight_halfspace_indices(data, idx1) ∩ tight_halfspace_indices(data, idx2)
    push!(tight, i)
    return add_adjacent_element!(data, i, i, newel, idx1, tight)
end

_shift(el::AbstractVector, line::Line) = el + Polyhedra.coord(line)
_shift(el::Line, line::Line) = el + line
_shift(el::Ray, line::Line) = el + Polyhedra.Ray(Polyhedra.coord(line))
function _λ_proj(r::VStruct, line::Line, h::HRepElement)
    # Line or ray `r`:
    # (r + λ * line) ⋅ h.a == 0
    # λ = -r ⋅ h.a / (line ⋅ h.a)
    return -r ⋅ h.a / (line ⋅ h.a)
end
function _λ_proj(x::AbstractVector, line::Line, h::HRepElement)
    # Point `x`:
    # (x + λ * line) ⋅ h.a == h.β
    # λ = (h.β - x ⋅ h.a) / (line ⋅ h.a)
    return (h.β - x ⋅ h.a) / (line ⋅ h.a)
end
function line_project(el, line, h)
    λ = _λ_proj(el, line, h)
    return Polyhedra.simplify(_shift(el, λ * line))
end
function hline(data, line::Line, i, h)
    value = h.a ⋅ line
    if !Polyhedra.isapproxzero(value)
        if data.cutline[i] === nothing
            # Make `lineray` point inward
            data.cutline[i] = value > 0 ? -line : line
            return true, line
        else
            line = line_project(line, data.cutline[i], h)
        end
    end
    data.nlines[i] += 1
    return false, line
end
# TODO remove solver arg `_`, it is kept to avoid breaking code
function doubledescription(hr::HRepresentation, _ = nothing)
    v = Polyhedra.dualfullspace(hr)
    hps = Polyhedra.lazy_collect(hyperplanes(hr))
    hss = Polyhedra.lazy_collect(halfspaces(hr))
    data = DoubleDescriptionData{pointtype(v), raytype(v), linetype(v)}(fulldim(hr), hps, hss)
    for line in lines(v)
        cut = false
        for i in reverse(eachindex(hps))
            cut, line = hline(data, line, nhalfspaces(hr) + i, hps[i])
            cut && break
        end
        if !cut
            for i in reverse(eachindex(hss))
                cut, line = hline(data, line, i, hss[i])
                cut && break
            end
        end
        if !cut
            push!(data.lines, line)
        end
    end
    # Add line rays after all lines are added so that the rays can be `line_project`ed.
    # We only do that for halfspaces, hyperplanes do not create rays from cutoff lines.
    # We use increasing index order since higher index may need the `lineray` of lower index.
    for i in eachindex(hss)
        line = data.cutline[i]
        if line !== nothing
            ray = Polyhedra.Ray(Polyhedra.coord(line))
            data.lineray[i] = add_element!(data, i - 1, ray, BitSet((i + 1):length(hss)))
        end
    end
    @assert isone(npoints(v))
    # Add the origin
    orig = project_onto_affspace(data, nhalfspaces(hr), first(points(v)), hps)
    if orig !== nothing
        add_element!(data, nhalfspaces(hr), orig, resized_bitset(data))
    end
    for i in reverse(eachindex(hss))
        if data.cutline[i] === nothing && isempty(data.cutpoints[i]) && isempty(data.cutrays[i])
            # Redundant, remove its contribution to allow more pruning with `is_adjacent_nc1`
            for p in data.pin[i]
                delete!(tight_halfspace_indices(data, p), i)
            end
            for r in data.rin[i]
                delete!(tight_halfspace_indices(data, r), i)
            end
        end
        if i > 1
            # Catches new adjacent rays, see 3.2 (ii) of [FP96]
            for p1 in data.pin[i], p2 in data.pin[i]
                if p1.cutoff < p2.cutoff
                    add_intersection!(data, p1, p2, i)
                end
            end
            for p in data.pin[i], r in data.rin[i]
                add_intersection!(data, p, r, i)
            end
        end
        deleteat!(data.cutpoints, i)
        deleteat!(data.cutpz, i)
        if i > 1
            for r1 in data.rin[i], r2 in data.rin[i]
                # We encounter both `r1, r2` and `r2, r1`.
                # Break this symmetry with:
                if r1.cutoff < r2.cutoff
                    add_intersection!(data, r1, r2, i)
                end
            end
        end
        deleteat!(data.cutrays, i)
        deleteat!(data.cutrz, i)
        deleteat!(data.pin, i)
        deleteat!(data.rin, i)
    end
    if isempty(data.points)
        # Empty polyhedron, there may be rays left,
        # Example 1: for 0x_1 + x_2 = -1 ∩ 0x_1 + x_2 = 1, the line (0, 1) is detected as correct
        # Example 2: for 0x_1 + 0x_2 = 1, the lines (1, 0) and (0, 1) are detected as correct
        # but since there is no point, the polyhedron is empty and we should drop all rays/lines
        empty!(data.lines)
        empty!(data.rays)
    end
    return similar(v, data.points, data.lines, data.rays)
end

function doubledescription(v::VRepresentation{T}, solver = nothing; kws...) where {T}
    checkvconsistency(v)
    lv = convert(LiftedVRepresentation{T, Matrix{T}}, v)
    R = -lv.R
    vl = doubledescription(MixedMatHRep{T}(R, zeros(T, size(R, 1)), lv.linset), solver; kws...)
    LiftedHRepresentation{T}(vl.R, vl.Rlinset)
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
# ------------ Simple gradients ------------ #

fn1 = () -> begin
    m = learnable(:m)
    q = trace(:q, Normal(m, 3.0))
    q
end

fn2 = () -> begin
    q = trace(:q, Normal(1.0, 3.0))
    q
end

grad_test_model = () -> begin
    N = 10
    q = trace(:fn1, fn1)
    z = trace(:fn2, fn2)
    z + q
end

@testset "Choice gradients" begin
    params = learnables([(:fn1, :m) => 5.0])
    ret, cl = simulate(params, grad_test_model)
    tg = target([(:fn2, :q)])
    vals, as, cgs = get_choice_gradients(tg, params, cl, 1.0)
    @test haskey(cgs, (:fn2, :q))
    @test as == ()
    @test haskey(vals, (:fn2, :q))
end

@testset "Learnable parameter gradients" begin
    params = learnables([(:fn1, :m) => 5.0])
    ret, cl = simulate(params, grad_test_model)
    vals, pgs = get_learnable_gradients(params, cl, 1.0)
    @test haskey(pgs, (:fn1, :m))
end

ext: jl


text:
using Test

function writepkg(name, precomp::Bool, submod::Bool)
    action = """
        global flag = true
    """

    if submod
        open("$(name)_submod.jl", "w") do io
            println(io, """
                export SubModule
                module SubModule
                    using Colors
                    flag = true
                end
            """)
        end

        action *= """
            include("$(name)_submod.jl")
        """
    end

    open("$name.jl", "w") do io
        println(io, """
__precompile__($precomp)

module $name

using Requires

flag = false

function __init__()
    @require Colors="5ae59095-9a9b-59fe-a467-6f913c188581" begin
        $(action)
    end
end

end
""")
    end
end

@testset "Requires" begin
    mktempdir() do pkgsdir
        local rm_CachedIncludeTest_submod_file
        cd(pkgsdir) do
            npcdir = joinpath("FooNPC", "src")
            mkpath(npcdir)
            cd(npcdir) do
                writepkg("FooNPC", false, false)
            end
            npcdir = joinpath("FooPC", "src")
            mkpath(npcdir)
            cd(npcdir) do
                writepkg("FooPC", true, false)
            end
            npcdir = joinpath("FooSubNPC", "src")
            mkpath(npcdir)
            cd(npcdir) do
                writepkg("FooSubNPC", false, true)
            end
            npcdir = joinpath("FooSubPC", "src")
            mkpath(npcdir)
            cd(npcdir) do
                writepkg("FooSubPC", true, true)
            end
            npcdir = joinpath("CachedIncludeTest", "src")
            mkpath(npcdir)
            cd(npcdir) do
                writepkg("CachedIncludeTest", true, true)
                submod_file = abspath("CachedIncludeTest_submod.jl")
                @test isfile(submod_file)
                rm_CachedIncludeTest_submod_file = ()->rm(submod_file)
            end
        end
        push!(LOAD_PATH, pkgsdir)

        @eval using FooNPC
        @test !FooNPC.flag
        @eval using FooPC
        @test !FooPC.flag
        @eval using FooSubNPC
        @test !(:SubModule in names(FooSubNPC))
        @eval using FooSubPC
        @test !(:SubModule in names(FooSubPC))
        @eval using CachedIncludeTest
        # Test that the content of the file which defines
        # CachedIncludeTest.SubModule is cached by `@require` so it can be used
        # even when the file itself is removed.
        rm_CachedIncludeTest_submod_file()
        @test !(:SubModule in names(CachedIncludeTest))

        @eval using Colors

        @test FooNPC.flag
        @test FooPC.flag
        @test :SubModule in names(FooSubNPC)
        @test FooSubNPC.SubModule.flag
        @test :SubModule in names(FooSubPC)
        @test FooSubPC.SubModule.flag
        @test :SubModule in names(CachedIncludeTest)

        cd(pkgsdir) do
            npcdir = joinpath("FooAfterNPC", "src")
            mkpath(npcdir)
            cd(npcdir) do
                writepkg("FooAfterNPC", false, false)
            end
            pcidr = joinpath("FooAfterPC", "src")
            mkpath(pcidr)
            cd(pcidr) do
                writepkg("FooAfterPC", true, false)
            end
            sanpcdir = joinpath("FooSubAfterNPC", "src")
            mkpath(sanpcdir)
            cd(sanpcdir) do
                writepkg("FooSubAfterNPC", false, true)
            end
            sapcdir = joinpath("FooSubAfterPC", "src")
            mkpath(sapcdir)
            cd(sapcdir) do
                writepkg("FooSubAfterPC", true, true)
            end
        end

        @eval using FooAfterNPC
        @eval using FooAfterPC
        @eval using FooSubAfterNPC
        @eval using FooSubAfterPC

        pop!(LOAD_PATH)

        @test FooAfterNPC.flag
        @test FooAfterPC.flag
        @test :SubModule in names(FooSubAfterNPC)
        @test FooSubAfterNPC.SubModule.flag
        @test :SubModule in names(FooSubAfterPC)
        @test FooSubAfterPC.SubModule.flag
    end
end

module EvalModule end

@testset "Notifications" begin
    push!(LOAD_PATH, joinpath(@__DIR__, "pkgs"))
    @eval using NotifyMe

    mktempdir() do pkgsdir
        ndir = joinpath("NotifyTarget", "src")
        cd(pkgsdir) do
            mkpath(ndir)
            cd(ndir) do
                open("NotifyTarget.jl", "w") do io
                    println(io, """
                    module NotifyTarget
                        using Requires
                        function __init__()
                            @require Example="7876af07-990d-54b4-ab0e-23690620f79a" begin
                                f(x) = 2x
                            end
                        end
                    end
                    """)
                end
            end
        end
        push!(LOAD_PATH, pkgsdir)
        @test isempty(NotifyMe.notified_args)
        @eval using NotifyTarget
        @test isempty(NotifyMe.notified_args)
        @eval using Example
        @test length(NotifyMe.notified_args) == 1
        nargs = NotifyMe.notified_args[1]
        @test nargs[1] == joinpath(pkgsdir, ndir, "NotifyTarget.jl")
        @test nargs[2] == NotifyTarget
        @test nargs[3] == "7876af07-990d-54b4-ab0e-23690620f79a"
        @test nargs[4] == "Example"
        Core.eval(EvalModule, nargs[5])
        @test Base.invokelatest(EvalModule.f, 3) == 6
    end

    pop!(LOAD_PATH)
end


text:
using DitherImage
using Base.Test
using Images
using TestImages

img = testimage("lena_gray_256")
imA = convert(Array{Float64, 2}, data(img))
imA_out = ditherimage(imA)
img_expected = load("lena_dither.png")
imA_expected = convert(Array{Float64, 2}, data(img_expected))

w = size(imA_out, 1)
h = size(imA_out, 2)
@test w == size(imA_expected, 1)
@test h == size(imA_expected, 2)

for i = 1:h
    for j = 1:w
        @test imA_out[w, h] == imA_expected[w, h]
    end
end


text:
using PyCall

rclpy = pyimport("rclpy")
rosNode = pyimport("rclpy.node")
std_msg = pyimport("std_msgs.msg")

@pydef mutable struct TimerNode <: rosNode.Node
    function __init__(self)
        rosNode.Node.__init__(self, "time")
	self.sent = [false, false, false]
        # publisher
        self.publisher_ = self.create_publisher(std_msg.Float64, "time", 10)
        timer_period = 1/500
        dt = 1/500
        self.i = 0
        self.t = 0.0
        function timer_callback(self)
            msg = std_msg.Float64()
            msg.data = self.t
	    # to avoid the communication delay due to JIT compilation
	    if self.sent == [false, false, false]
		    self.sent = [true, false, false]
	    elseif self.sent == [true, false, false]
		    self.sent = [true, true, false]
	    elseif self.sent == [true, true, false]
		    self.sent = [true, true, true]
		    # sleep_time = 10
		    # self.get_logger().info("sleep time: $(sleep_time) [s]. Wait a moment plz...")
		    # sleep(sleep_time)
	    end
            self.publisher_.publish(msg)  # publish
            self.get_logger().info("time: $(self.t)")
            sleep(0.001)  # too fast publishing may yield divergence of simulation due to communication delay
            self.i = self.i + 1
            self.t = self.t + dt
        end
        self.timer = self.create_timer(timer_period, () -> timer_callback(self))
    end
end

function main(args=nothing)
     rclpy.init(args=args)
     println("timer node started")
     timer_node = TimerNode()
     rclpy.spin(timer_node)
     timer_node.destroy_node()
     rclpy.shutdown()
end

main()


text:
"""
```
t = find_threshold(MinimumError(), histogram, edges)
```

Under the assumption that the histogram is a mixture of two Gaussian
distributions the threshold is chosen such that the expected misclassification
error rate is minimised.

# Output

Returns a real number `t` that specifies the threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(MinimumError(), counts[1:end], edges)
```

# References

1. J. Kittler and J. Illingworth, “Minimum error thresholding,” Pattern Recognition, vol. 19, no. 1, pp. 41–47, Jan. 1986. [doi:10.1016/0031-3203(86)90030-0](https://doi.org/10.1016/0031-3203%2886%2990030-0)
2. Q.-Z. Ye and P.-E. Danielsson, “On minimum error thresholding and its implementations,” Pattern Recognition Letters, vol. 7, no. 4, pp. 201–206, Apr. 1988. [doi:10.1016/0167-8655(88)90103-1](https://doi.org/10.1016/0167-8655%2888%2990103-1)
"""
function find_threshold(algorithm::MinimumError, histogram::AbstractArray, edges::AbstractRange)
    frequencies = cumsum(histogram)
    X = cumsum(edges .* histogram)
    X² = cumsum(edges.^2 .* histogram)
    𝔼X₁, 𝔼X₁², 𝔼X₂, 𝔼X₂² = compute_expectations(frequencies, X, X²)
    minJ = typemax(Float64)
    index = firstindex(frequencies)
    for t =  firstindex(frequencies):lastindex(frequencies)-1
        P₁ = max(frequencies[t] / frequencies[end], eps())
        P₂ = 1 - P₁
        # Var(X) = E[X²] - E[X]² = E[X²] - μ²
        σ₁ = sqrt(max(𝔼X₁²[t] - 𝔼X₁[t]^2, eps()))
        σ₂ = sqrt(max(𝔼X₂²[t] - 𝔼X₂[t]^2, eps()))
        J = 1 + 2*(P₁*log(σ₁) + P₂*log(σ₂)) - 2*(P₁*log(P₁)+ P₂*log(P₂))
        if (J < minJ && isfinite(J))
            minJ = J
            index = t
        end
    end
    edges[index]
end

#=
  If the histogram consists of integers then we risk integer overflow when we
  take the cumulative sums. Hence, we explicitlty convert the histogram to
  Float.
=#
function find_threshold(algorithm::MinimumError, histogram::AbstractArray{T}, edges::AbstractRange) where T <: Int
    find_threshold(algorithm, convert.(Float64, histogram), edges)
end

function compute_expectations(frequencies, X, X²)
    dimensions = axes(frequencies)
    𝔼X₁ = zeros(dimensions)
    𝔼X₁² = zeros(dimensions)
    𝔼X₂ = zeros(dimensions)
    𝔼X₂² = zeros(dimensions)
    for i = firstindex(X):lastindex(X)-1
        𝔼X₁[i] = X[i] / frequencies[i]
        𝔼X₁²[i] = X²[i] / frequencies[i]
        𝔼X₂[i] = (X[end] - X[i]) /  (frequencies[end] - frequencies[i])
        𝔼X₂²[i] = (X²[end] - X²[i]) /  (frequencies[end] - frequencies[i])
    end
    𝔼X₁, 𝔼X₁², 𝔼X₂, 𝔼X₂²
end


text:
using Test
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

@testset "Gromacs" begin

  #
  # Tests with Gromacs-XSC trajectory
  #

  dir="./data/Gromacs"
  atoms = readPDB("$dir/system.pdb")
  options = Options(stride=5,seed=321,StableRNG=true,nthreads=1,silent=true)

  # Example 1: protein-EMIM

  # save(R,"$dir/protein_EMI.json")
  R_save = load("$dir/protein_EMI.json")
  protein = Selection(select(atoms,"protein"),nmols=1)
  emi = Selection(select(atoms,"resname EMI"),natomspermol=20)
  traj = Trajectory("$dir/trajectory.xtc",protein,emi)
  R = mddf(traj,options)
  @test isapprox(R,R_save,debug=true)

  # Example 1: EMIM-DCA

  # save(R,"$dir/EMI_DCA.json")
  R_save = load("$dir/EMI_DCA.json")
  emi = Selection(select(atoms,"resname EMI"),natomspermol=20)
  dca = Selection(select(atoms,"resname NC"),natomspermol=5)
  traj = Trajectory("$dir/trajectory.xtc",emi,dca)
  R = mddf(traj,options)
  @test isapprox(R,R_save,debug=true)

  # Example 1: EMIM-EMIM

  # save(R,"$dir/EMI_EMI.json")
  R_save = load("$dir/EMI_EMI.json")
  emi = Selection(select(atoms,"resname EMI"),natomspermol=20)
  traj = Trajectory("$dir/trajectory.xtc",emi)
  R = mddf(traj,options)
  @test isapprox(R,R_save,debug=true)

end


text:
module KTBC
using Test

import Pkg
Pkg.add("ToeplitzMatrices")

import ToeplitzMatrices.Toeplitz
import LinearAlgebra.Tridiagonal

function CreateKTBC(dims=2)
    toeplitzArgs = [2, -1, zeros(dims - 2)...]
    K = Tridiagonal(Toeplitz(toeplitzArgs, toeplitzArgs))
    T = copy(K)
    T[1, 1] = 1
    B = copy(T)
    B[dims, dims] = 1
    C = convert(Array{Float64}, K)
    C[1, dims] = -1
    C[dims, 1] = -1
    return K, T, B, C
end

export CreateKTBC

K, T, B, C = CreateKTBC(3)

@test K == [
  2.0 -1.0 0.0
  -1.0 2.0 -1.0
  0.0 -1.0 2.0
]

@test T == [
  1.0 -1.0 0.0
  -1.0 2.0 -1.0
  0.0 -1.0 2.0
]

@test B == [
  1.0 -1.0 0.0
  -1.0 2.0 -1.0
  0.0 -1.0 1.0
]

@test C == [
  2.0 -1.0 -1.0
  -1.0 2.0 -1.0
  -1.0 -1.0 2.0
]
end


text:
import FireflyModule.Firefly

function arc_distance(f1::Firefly, f2::Firefly)::Int
    """Calculates the arc distance between two firefly paths.
    Example:
        Firefly i -> 1 2 4 6 3 5 7
        Firefly j -> 1 2 4 3 5 7 6
        1-2, 2-4, 3-5, 5-7 arcs in firefly i are also presented
        in firefly j. But 4-6, 6-3 arcs are missing in
        firefly j. so the arc distance between firefly i
        and firefly j is 2 here.
    """
    missmatched_pair_counts = 0
    f1_pairs = zip(f1.path[1:end-1], f1.path[2:end])
    f2_pairs = zip(f2.path[1:end-1], f2.path[2:end])

    for f1_pair in f1_pairs
        if f1_pair ∉ f2_pairs missmatched_pair_counts += 1 end
    end

    return missmatched_pair_counts
end

max_stars_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_issues_repo_name: zahidkizmaz/firefly-algorithm-tsp

max_forks_repo_name: zahidkizmaz/firefly-algorithm-tsp


text:
"""
    Catenary(trans::AffineMap, xmin, xmax, a)

Construct a `Catenary` object. The `AffineMap` brings the catenary into a frame
frame in which the `x` axis is aligned with vertically hanging catenary. In this frame,
the `x` values of the catenary lie between `xmin` and `xmax`, the `y` component is zero,
and `z` is parameterized by `x` as:

z = a * (cosh(x / a) - 1)
"""
struct Catenary{Tfm <: AffineTransformation, T <: Real}
    transform::Tfm # transformation from global frame to quadratic's frame
    lmin::T  # "Longitudinal" coordinate at end of catenary
    lmax::T  # "Longitudinal" coordinate at start of catenary
    a::T # Catenary paramter equivalent to 1/2a (such that ax^2+bx+c)
end

"""
    Catenary(Θ, x, y, z, lmin, lmax, a)

Create a catenary from the parameters stored in the Roames database.

### Arguments
* `Θ` : Angle of catenary in XY plane
* `x`: x position of lowest point on catenary
* `y`: y position of lowest point on catenary
* `z`: z position of lowest point on catenary
* `lmin`: Longitudinal coordinate at start of catenary
* `lmax`: Longitudinal coordinate at end of catenary
* `a`: Catenary paramter equivalent to 1/2a (such that ax^2+bx+c)
"""
function Catenary(Θ::T, x::T, y::T, z::T, lmin::T, lmax::T, a::T) where {T}
    trans = inv(Translation(SVector(x, y, z)) ∘ LinearMap(RotZ(Θ)))
    return Catenary(trans, lmin, lmax, a)
end
@inline Catenary(args...) = Catenary(promote(args...)...)

eltype(::Catenary{T}) where {T} = SVector{3,T}
eltype(::Type{<:Catenary{T}}) where {T} = SVector{3,T}

# Use getindex to get points along the line
function getindex(c::Catenary, x)
    inv(c.transform)(SVector(x, 0, c.a * (cosh(x / c.a) - 1)))
end
length(c::Catenary) = c.lmax - c.lmin
#lastindex(c::Catenary) = c.lmax # seems unsafe with `begin`

getindex(c::Catenary, v::AbstractVector) = map(x -> c[x], v)

# Transformations (which compose into AffineMap's)
# Note: Julia's dispatch system makes it difficult to support all `Transformations`.
(trans::Translation{V})(c::Catenary) where {V} = Catenary(c.transform ∘ inv(trans), c.lmin, c.lmax, c.a)
(trans::LinearMap{M})(c::Catenary) where {M} = Catenary(c.transform ∘ inv(trans), c.lmin, c.lmax, c.a)
(trans::AffineMap{M,V})(c::Catenary) where {M,V} = Catenary(c.transform ∘ inv(trans), c.lmin, c.lmax, c.a)

"""
    Quadratic(c::Catenary)

Create a quadratic which approximates `c`. This is an approximate operation, however
the fitted quadratic is tuned for the relevant segment of the catenary and
`Catenary(Quadratic(c))` should be equal to `c` up to numerical error.
"""
function Quadratic(c::Catenary)
    Quadratic(c[c.lmin], c[0.5*(c.lmin + c.lmax)], c[c.lmax])
end

"""
    Catenary(q::Quadratic)

Create a quadratic which approximates `q`. This is an approximate operation, however
the fitted catenary is tuned for the relevant segment of the quadratic and
`Quadratic(Catenary(q))` should be equal to `q` up to numerical error.
"""
function Catenary(q::Quadratic)
    Catenary(q[0], q[0.5*end], q[end])
end

"""
    Catenary(p_start, p_mid, p_end)

Takes three coordinates for the start, somewhere in the middle, and the end of a catenary
and constructs the appropriate `Catenary` object. A numerical fitting method is used to
fit the catenary parameters.
"""
function Catenary(p_start::StaticVector{3,T}, p_mid::StaticVector{3,T}, p_end::StaticVector{3,T}) where {T}
    dir = SVector(p_end[1] - p_start[1], p_end[2] - p_start[2])
    len = norm(dir) # horizontal length
    @static if VERSION < v"0.7"
        θ = atan2(dir[2], dir[1]) # rotation in x-y plane
    else
        θ = atan(dir[2], dir[1]) # rotation in x-y plane
    end
    tform = inv(AffineMap(RotZ(θ), p_start))

    # Now fit a catenary in the x-z plane. Ignore y value of midpoint on the basis that
    # `Catenary` does not swing (unlike `Quadratic`)
    p_mid_trans = tform(p_mid)
    x_mid = p_mid_trans[1]
    z_mid = p_mid_trans[3]

    x_end = len
    z_end = p_end[3] - p_start[3]

    (a, x0, z0) = fit_catenary_origin_2d(x_mid, z_mid, x_end, z_end)

    # Bring transformation to midpoint
    tform2 = Translation(SVector(-x0, zero(T), -z0)) ∘ tform

    return Catenary(tform2, -x0, -x0+len, a)
end

function fit_catenary_origin_2d(x1::T, z1::T, x2::T, z2::T) where {T}
    # Note: one point is going through the origin

    # First guess - make a quadratic approximation: z = q_a*x^2 + q_b*x + q_c
    q_a = (z2/x2 - z1/x1)/(x2 - x1)
    q_b = (x1*z2/x2 - x2*z1/x1)/(x1 - x2)
    # q_c = 0

    # Find the catenary which matches best at the minimum

    # defend against straight catenaries and degeneracies
    if isnan(q_a) || abs(q_a*max(x1,x2)) < 1e-6 # Second term selects cases where quadratic term changes z by less than 1e-6th of it's length
        # It's either flat or degenerate

        # Choose a bottom which is one million times further than points
        q_b = z2/x2 # slope

        # allow x0 to be close when it's (almost?) flat, but force it to be far when slope is large
        x0 = -1e6*z2 # = -1e6*x2*sign(z2)*abs(q_b)
        a = (q_b == zero(T)) ? inv(2*eps(T)) : -x0/q_b

        # Iteratively improve (as below, but optimizing only `a` w.r.t. `z2`)
        for i = 1:10
            Δz₂ = a*(cosh((x2-x0)/a) - cosh(-x0/a)) - z2
            dΔz₂_da = (cosh((x2-x0)/a) - cosh(-x0/a)) - (((x2-x0)/a)*sinh((x2-x0)/a) - (-x0)/a*sinh(-x0/a))
            a = a - dΔz₂_da \ Δz₂
        end

        z0 = a*(1 - cosh(-x0/a))

        if !isfinite(z0)||!isfinite(x0)||!isfinite(a) # On rare occassions this path proves unstable when the standard one succeeds
            # Continue with standard path
            q_b = (x1*z2/x2 - x2*z1/x1)/(x1 - x2)
        else
            return (a, x0, z0)
        end
    end

    a = inv(2*q_a)
    x0 = -q_b * a # -q_b/(2*q_a)
    # z0 = q_a*x0*x0 + q_b*x0 - a

    # Now iteratively improve upon this guess, fixing z0 = a(1-cosh(-x0/a))
    for i = 1:10
        Δz₁ = a*(cosh((x1-x0)/a) - cosh(-x0/a)) - z1
        Δz₂ = a*(cosh((x2-x0)/a) - cosh(-x0/a)) - z2
        dΔz₁_dx0 = -(sinh((x1-x0)/a) - sinh(-x0/a))
        dΔz₂_dx0 = -(sinh((x2-x0)/a) - sinh(-x0/a))
        dΔz₁_da = (cosh((x1-x0)/a) - cosh(-x0/a)) - (((x1-x0)/a)*sinh((x1-x0)/a) - (-x0)/a*sinh(-x0/a))
        dΔz₂_da = (cosh((x2-x0)/a) - cosh(-x0/a)) - (((x2-x0)/a)*sinh((x2-x0)/a) - (-x0)/a*sinh(-x0/a))

        J = @SMatrix [dΔz₁_dx0 dΔz₁_da; dΔz₂_dx0 dΔz₂_da]
        x = @SVector [Δz₁, Δz₂]
        sol = J \ x
        x0 = x0 - sol[1]
        a = a - sol[2]

        if Δz₁*Δz₁ + Δz₂*Δz₂ < eps(T)
            break
        end
    end

    z0 = a*(1 - cosh(-x0/a))

    return (a, x0, z0)
end


text:
@time @testset "Reachability.options" begin include("unit_options.jl") end
@time @testset "Reachability.solve_continuous" begin include("solve_continuous.jl") end
@time @testset "Reachability.solve_hybrid" begin include("solve_hybrid.jl") end


text:
using Plots, LaTeXStrings, Plots.PlotMeasures, Statistics
using MAT

include("SDOF_Model.jl")

(@isdefined TARGET_FOLDER) ? nothing : TARGET_FOLDER = ""

## Load data

vars = matread(joinpath(@__DIR__, "PE_AD_numerical.mat"))

α              = vars["alphas"]'
N              = vars["N"]'               # number of periods used in the computation
numpers        = vars["npers"]'           # number of periods of the signal
AD_Bathe_teo   = vars["AD_Bathe_teo"]
AD_Bathe_num   = vars["AD_Bathe_num"]
AD_Newmark_teo = vars["AD_Newmark_teo"]
AD_Newmark_num = vars["AD_Newmark_num"]
PE_Bathe       = vars["vectorPETsOctave_Bathe"]
PE_Newmark     = vars["vectorPETsOctave_Newmark"]
PE_Bathe_teo   = vars["PE_Bathe_teo"]
PE_Newmark_teo = vars["PE_Newmark_teo"]

T = 1/2
PE_RFEM_inf, PE_RFEM_sup = compute_period_elongation(T, N, numpers, α)
AD_RFEM_inf = compute_amplitude_decay(T, N, numpers, α)

## Period elongation

function figure_period_elongation()
    fig = plot(xlab=L"\alpha", ylab=L"\% PE", legend=:topleft,
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([0.0, 0.05, 0.10, 0.15, 0.20], [L"0.0", L"0.05", L"0.10", L"0.15", L"0.20"]),
               ytick=([0.0, 2.5, 5.0, 7.5, 10.0], [L"0.0", L"2.5", L"5.0", L"7.5", L"10.0"]),
               xlims=(0.0, 0.22), ylims=(-0.3, 11.0),
               bottom_margin=6mm, left_margin=20mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    plot!(fig, α[1:2:end], PE_Bathe[1:2:end], c=:green, linestyle=:solid, lw=2.5, markersize=6, lab=L"\textrm{Bathe}", marker=:utriangle)
    plot!(fig, α[1:2:end], PE_Newmark[1:2:end], c=:red, linestyle=:solid, lw=2.5, markersize=6, lab=L"\textrm{Newmark}", marker=:circle)
    plot!(fig, α[1:2:end], PE_RFEM_sup[1:2:end], c=:blue, linestyle=:solid, lw=2.5, markersize=6, lab=L"\textrm{Flowpipe}", marker=:square)
    fig
end

fig = figure_period_elongation()
savefig(fig, joinpath(TARGET_FOLDER, "PE_all.pdf"))

## Amplitude decay

function figure_amplitude_decay()
    fig = plot(xlab=L"\alpha", ylab=L"\% AD", legend=:topleft,
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([0.0, 0.05, 0.10, 0.15, 0.20], [L"0.0", L"0.05", L"0.10", L"0.15", L"0.20"]),
               ytick=([0.0, 1.0, 2.0, 3.0], [L"0.0", L"1.0", L"2.0", L"3.0"]),
               xlims=(0.0, 0.22), ylims=(-0.3, 3.2),
               bottom_margin=6mm, left_margin=20mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    plot!(fig, α[1:2:end], AD_Bathe_num[1:2:end], c=:green, linestyle=:solid, lw=2.5, markersize=6, lab=L"\textrm{Bathe}", marker=:utriangle)
    plot!(fig, α[1:2:end], AD_Newmark_num[1:2:end], c=:red, linestyle=:solid, lw=2.5, markersize=6, lab=L"\textrm{Newmark}", marker=:circle)
    plot!(fig, α[1:2:end], AD_RFEM_inf[1:2:end], c=:blue, linestyle=:solid, lw=2.5, markersize=6, lab=L"\textrm{Flowpipe}", marker=:square)
    fig
end
fig = figure_amplitude_decay()
savefig(fig, joinpath(TARGET_FOLDER, "AD_all.pdf"))

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
module Talkon

using Telegram
using Telegram.API
using Serialization
using UnPack
using Dates
using Random

include("database.jl")
include("begin.jl")
include("request.jl")
include("enter.jl")
include("quality.jl")
include("suportive.jl")

export talka, initialize

const DATAFILE = Ref("")

function initialize(datafile = "varTEST.data")
    DATAFILE[] = datafile
    (T, Av, update_id, K, groups) = deserialize(datafile)
    return DataBase(T, Av, update_id, K, groups)
end

function talka(d::DataBase)
    update_id = d.update_id
    av = nothing
    tg = TelegramClient()
    cnt = 0
    @info "Talkon bot started"
    while true
        cnt += 1
        update_id, av = dothing(d, tg, update_id, av)
        d.update_id = update_id
        if cnt == 20
            av = chekeis(d, tg, av)
            cnt = 0
        end
    end
end # end of talka function

function dothing(d::DataBase, tg, update_id, av)
    @unpack T, Av, K, groups = d
    mst = try
        mst=getUpdates(tg, offset = update_id + 1, allowed_updates = ["message", "callback_query"])
    catch
        @warn "getUpdates doesnt work may be no connection?!"
        sleep(2)
        return update_id, av
    end
    if length(mst) > 1
        print("length of the message $(length(mst))")
    elseif length(mst) == 0
        sleep(1)
        return update_id, av
    end
    update_id = mst[end]["update_id"]
    for ms in mst
        if haskey(ms, "callback_query") #messages that have callback string
            id = ms["callback_query"]["from"]["id"]
            av = Av[id]
            av["txt"] = "nothingness"
            av = switcher(d, tg, av, ms[:callback_query][:data])
            #eval(Meta.parse(ms[:callback_query][:data]))
            chat_id = ms["callback_query"]["message"]["chat"]["id"]
            try
                deleteMessage(chat_id = chat_id, message_id = ms["callback_query"]["message"]["message_id"])
            catch
            end
            continue
        elseif haskey(ms["message"], "new_chat_member")
            #if ms["message"]["chat"]["id"] ==-1001715161879# development group. https://t.me/talkon_development
            if ms["message"]["chat"]["id"] == -1001547960563 #message from main test group https://t.me/talka_ukrainai
                if !haskey(Av, ms["message"]["new_chat_member"]["id"]) # new member entered main group
                    try
                    av = newav(d, tg, av, ms) # register and wellcome message for a new member
                    catch
                        println("unsuccesfull registration upon entering main group")
                    end
                end
                continue
            end
            promoteChatMember(chat_id = ms["message"]["chat"]["id"], #new member entered the temporary group
                                user_id = ms["message"]["new_chat_member"]["id"], can_manage_voice_chats = true)
            continue
        elseif !haskey(ms["message"],"text")
            println("unknown type of message")
            continue
        end
        if !haskey(Av, ms["message"]["from"]["id"]) #new member called bot
            av = newav(d, tg, av, ms)
            continue
        end
        av = Av[ms["message"]["from"]["id"]]
        ReplyKeyboardRemove = Dict(:remove_keyboard => true)
        sentms = sendMessage(chat_id = av["id"], text = "entered", reply_markup = ReplyKeyboardRemove)  # sendMessage always require text message,
        deleteMessage(chat_id = av["id"],message_id = sentms["message_id"]) # in this case it has to be deleted
        av["txt"]=replace(ms["message"]["text"],"✓"=>"") # curently tree is represented not with inline keyboard, thus the text message arrives
        if lowercase(av["txt"]) == "namo | додому" || av["txt"] == "/start" || av["txt"] == "/namo"
            tbegin(tg, av)
        #elseif av["txt"] == "/pakviesti" #invite new member, perhaps to help refugees it is not very usefull
        #    invite()
        #elseif av["txt"] == "/taskai"
        #    taskai()
        #elseif av["step"] == "wait_for_invite_code"
        #    registerinvitecode()
        elseif av["step"] == "request"
            trequest(d, tg, av)
        elseif av["step"] == "enter"
            tenter(d, tg, av)
        elseif av["step"] == "requested" || av["step"] == "accepted"
            dealkeis(d, tg, av)
        end
    end # end of for ms
    serialize(DATAFILE[], [T,Av,update_id,K,groups])
    return update_id, av
end # end of dothing function

function switcher(d, tg, av, inp)
    if inp == "tenter()"
        tenter(d, tg, av)
    elseif inp == "trequest()"
        trequest(d, tg, av)
    elseif inp == "valuableyesno(0)"
        valuableyesno(d, tg, av, 0)
    elseif inp == "valuableyesno(1)"
        valuableyesno(d, tg, av, 1)
    elseif inp == "invitecodeyesno(0)"
        invitecodeyesno(0)
    elseif inp == "invitecodeyesno(1)"
        invitecodeyesno(1)
    elseif inp == "closekeis()"
        av = closekeis(d, tg, av, av["requestid"])
    elseif inp == ""
    else
        println("bad command!")
    end

    return av
end

end # module


text:
module CxxWrap

module CxxWrapCore

import Libdl
import MacroTools

export @wrapmodule, @readmodule, @wraptypes, @wrapfunctions, @safe_cfunction, @initcxx,
ConstCxxPtr, ConstCxxRef, CxxRef, CxxPtr,
CppEnum, ConstArray, CxxBool, CxxLong, CxxULong, CxxChar, CxxWchar, CxxUChar, CxxSignedChar, CxxLongLong, CxxULongLong,
ptrunion, gcprotect, gcunprotect, isnull, libcxxwrapversion

const libcxxwrap_version_range = (v"0.7.0",  v"0.8")

using libcxxwrap_julia_jll # for libcxxwrap_julia and libcxxwrap_julia_stl

# These can't be products, since we want to control how and when they are dlopened
for libname in ["jlcxx_containers", "except", "extended", "functions", "hello", "basic_types", "inheritance", "parametric", "pointer_modification", "types"]
  libcxxwrap_julia_name = basename(libcxxwrap_julia_jll.libcxxwrap_julia_path)
  libprefix = startswith(libcxxwrap_julia_name, "lib") ? "lib" : ""
  libext = libcxxwrap_julia_name[findlast('.', libcxxwrap_julia_name):end]
  full_libname = libprefix * libname * libext
  symname = "lib"*libname
  @eval $(Symbol(symname))() = joinpath(dirname(libcxxwrap_julia_jll.libcxxwrap_julia_path), $(full_libname))
end

prefix_path() = dirname(dirname(libcxxwrap_julia_jll.libcxxwrap_julia_path))

libcxxwrapversion() = VersionNumber(unsafe_string(ccall((:cxxwrap_version_string,libcxxwrap_julia), Cstring, ())))

function checkversion()
  jlcxxversion = libcxxwrapversion()
  if !(libcxxwrap_version_range[1] <= jlcxxversion < libcxxwrap_version_range[2])
    error("This version of CxxWrap requires a libcxxwrap-julia in the range $(libcxxwrap_version_range), but version $jlcxxversion was found")
  end
end

# Must also be called during precompile
checkversion()

# Welcome to the C/C++ integer type mess
# See https://en.cppreference.com/w/cpp/language/types and https://en.cppreference.com/w/cpp/types/integer

abstract type CxxSigned <: Signed end
abstract type CxxUnsigned <: Unsigned end

primitive type CxxBool <: CxxUnsigned 8*sizeof(Cuchar) end
const CharSigning = supertype(Cchar) == Signed ? CxxSigned : CxxUnsigned
primitive type CxxChar <: CharSigning 8*sizeof(Cchar) end
const WCharSigning = supertype(Cwchar_t) == Signed ? CxxSigned : CxxUnsigned
primitive type CxxWchar <: WCharSigning 8*sizeof(Cwchar_t) end

function _transform_fundamental_type(fundamental_type_name)
  julianame = "Cxx"
  for part in split(fundamental_type_name)
    if part == "unsigned"
      julianame *= "U"
      continue
    end
    julianame *= titlecase(part)
  end
  return Symbol(julianame)
end

_transform_fixed_type(fixed_type_name) = Symbol(replace(titlecase(fixed_type_name), "int" => "Int")[1:end-2])

# This macro adds the fundamental integer types such as long, which becomes CxxLong
# Names are e.g. CxxInt32 or CxxUInt64
macro add_int_types()
  all_fundamental_types = String[]
  type_sizes = Any[]
  fundamental_types_matched = String[]
  equivalent_types = String[]
  ccall((:get_integer_types,libcxxwrap_julia), Cvoid, (Any,Any,Any,Any), all_fundamental_types, type_sizes, fundamental_types_matched, equivalent_types)
  @assert all(all_fundamental_types .!= "undefined")
  sizedict = Dict(all_fundamental_types .=> type_sizes)
  result = quote end
  for (fundamentaltype, fixedtype) in zip(fundamental_types_matched, equivalent_types)
    push!(result.args, esc(:(const $(_transform_fundamental_type(fundamentaltype)) = $(_transform_fixed_type(fixedtype)))))
  end
  missing_types = setdiff(all_fundamental_types, fundamental_types_matched)
  for missingtype in missing_types
    super = startswith(missingtype, "unsigned") ? CxxUnsigned : CxxSigned
    nbits = sizedict[missingtype]*8
    push!(result.args, esc(:(primitive type $(_transform_fundamental_type(missingtype)) <: $super $nbits end)))
  end
  return result
end
@add_int_types

# Get the equivalent Julia type for a Cxx integer type
@generated julia_int_type(::Type{T}) where {T<:CxxSigned} = Symbol(:Int, 8*sizeof(T))
@generated julia_int_type(::Type{T}) where {T<:CxxUnsigned} = Symbol(:UInt, 8*sizeof(T))

to_julia_int(x::Union{CxxSigned,CxxUnsigned}) = reinterpret(julia_int_type(typeof(x)),x)

const CxxNumber = Union{CxxSigned,CxxUnsigned}
Base.show(io::IO, n::CxxNumber) = show(io, to_julia_int(n))
Base.show(io::IO, b::CxxBool) = show(io, Bool(b))
function Base.promote_rule(::Type{CT}, ::Type{JT}) where {CT <: CxxNumber, JT <: Number}
  if julia_int_type(CT) == JT
    return JT
  end
  return Base.promote_rule(julia_int_type(CT), JT)
end
Base.promote_type(::Type{T}, ::Type{T}) where {T<:CxxNumber} = julia_int_type(T)
Base.promote_rule(::Type{T}, ::Type{T}) where {T<:CxxNumber} = Base.promote_type(T,T)
Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:CxxNumber, T2<:CxxNumber} = Base.promote_rule(julia_int_type(T1), julia_int_type(T2))
Base.AbstractFloat(x::CxxNumber) = Base.AbstractFloat(to_julia_int(x))

# Conversion to and from the equivalent Julia type
Base.convert(::Type{T}, x::Number) where {T<:Union{CxxSigned,CxxUnsigned}} = reinterpret(T, convert(julia_int_type(T), x))
Base.convert(::Type{JT}, x::CT) where {JT<:Number,CT<:Union{CxxSigned,CxxUnsigned}} = convert(JT,reinterpret(julia_int_type(CT), x))
Base.convert(::Type{T}, x::CT) where {T <: Union{CxxWrapCore.CxxSigned, CxxWrapCore.CxxUnsigned}, CT <: Union{CxxWrapCore.CxxSigned, CxxWrapCore.CxxUnsigned}}  = convert(T,reinterpret(julia_int_type(CT), x))

# Convenience constructors
(::Type{T})(x) where {T<:Union{CxxWrapCore.CxxSigned,CxxWrapCore.CxxUnsigned}} = convert(T,x)

Base.flipsign(x::T, y::T) where {T <: CxxSigned} = reinterpret(T, flipsign(to_julia_int(x), to_julia_int(y)))

# Trait type to indicate a type is a C++-wrapped type
struct IsCxxType end
struct IsNormalType end

@inline cpp_trait_type(::Type) = IsNormalType

# Enum type interface
abstract type CppEnum <: Integer end
Base.convert(::Type{T}, x::CppEnum) where {T <: Integer} = T(reinterpret(Int32, x))
Base.convert(::Type{T}, x::Integer) where {T <: CppEnum} = reinterpret(T, Int32(x))
Base.convert(::Type{T}, x::T) where {T <: CppEnum} = x
import Base: +, |
+(a::T, b::T) where {T <: CppEnum} = reinterpret(T, convert(Int32,a) + convert(Int32,b))
|(a::T, b::T) where {T <: CppEnum} = reinterpret(T, convert(Int32,a) | convert(Int32,b))

"""
Base class for smart pointers
"""
abstract type SmartPointer{T} end
@inline cpp_trait_type(::Type{SmartPointer{T}}) where {T} = IsCxxType

Base.show(io::IO, x::SmartPointer) = print(io, "C++ smart pointer of type ", typeof(x))

allocated_type(t::Type) = Any
dereferenced_type(t::Type) = Any

__cxxwrap_smartptr_dereference(p::SmartPointer{T}) where {T} = __cxxwrap_smartptr_dereference(CxxRef(p))
__cxxwrap_smartptr_construct_from_other(t::Type{<:SmartPointer{T}}, p::SmartPointer{T}) where {T} = __cxxwrap_smartptr_construct_from_other(t,CxxRef(p))
__cxxwrap_smartptr_cast_to_base(p::SmartPointer{T}) where {T} = __cxxwrap_smartptr_cast_to_base(CxxRef(p))

function Base.getindex(p::SmartPointer{T}) where {T}
  return __cxxwrap_smartptr_dereference(p)
end

# No conversion if source and target type are identical
Base.convert(::Type{T}, p::T) where {PT,T <: SmartPointer{PT}} = p

# Construct from a related pointer, e.g. a std::weak_ptr from std::shared_ptr
function Base.convert(::Type{T1}, p::T2) where {T, T1 <: SmartPointer{T}, T2 <: SmartPointer{T}}
  return __cxxwrap_smartptr_construct_from_other(T1, p)
end

# upcast to base class
function Base.convert(::Type{T1}, p::T2) where {BaseT,DerivedT, T1 <: SmartPointer{BaseT}, T2 <: SmartPointer{DerivedT}}
  if !(DerivedT <: BaseT)
    error("$DerivedT does not inherit from $BaseT in smart pointer convert")
  end
  # First convert to base type
  base_p = __cxxwrap_smartptr_cast_to_base(p)
  return convert(T1, base_p)
end

struct StrictlyTypedNumber{NumberT}
  value::NumberT
end
function Base.convert(t::Type{<:StrictlyTypedNumber}, n::Number)
  @assert t == StrictlyTypedNumber{typeof(n)}
  return StrictlyTypedNumber{typeof(n)}(n)
end
Base.convert(t::Type{StrictlyTypedNumber{CxxBool}}, b::Bool) = StrictlyTypedNumber{CxxBool}(b)

abstract type CxxBaseRef{T} <: Ref{T} end

struct CxxPtr{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  CxxPtr{T}(x::Ptr) where {T} = new{T}(x)
  CxxPtr{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

struct ConstCxxPtr{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  ConstCxxPtr{T}(x::Ptr) where {T} = new{T}(x)
  ConstCxxPtr{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

struct CxxRef{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  CxxRef{T}(x::Ptr) where {T} = new{T}(x)
  CxxRef{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

struct ConstCxxRef{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  ConstCxxRef{T}(x::Ptr) where {T} = new{T}(x)
  ConstCxxRef{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

_ref_type(::Type{RefT}, ::Type{<:CxxBaseRef{T}}) where {RefT,T} = _ref_type(RefT,T)
_ref_type(::Type{RefT}, ::Type{T}) where {RefT,T} = _ref_type(RefT, T, cpp_trait_type(T))
_ref_type(::Type{RefT}, ::Type{T}, ::Type) where {RefT,T} = error("Manual Cxx Reference creation is only for C++ types")
function _ref_type(::Type{RefT}, ::Type{T}, ::Type{IsCxxType}) where {RefT,T}
  if isconcretetype(T)
    return RefT{supertype(T)}
  end
  return RefT{T}
end
_make_ref(::Type{RefT}, x::T) where {RefT,T} = _ref_type(RefT,T)(x.cpp_object)

CxxPtr(x) = _make_ref(CxxPtr,x)
ConstCxxPtr(x) = _make_ref(ConstCxxPtr,x)
CxxRef(x) = _make_ref(CxxRef,x)
ConstCxxRef(x) = _make_ref(ConstCxxRef,x)

Base.convert(t::Type{<:CxxBaseRef{T}}, x::Ptr{NT}) where {T <: Union{CxxUnsigned, CxxSigned}, NT <: Integer} = t(reinterpret(Ptr{T}, x))

Base.:(==)(a::Union{CxxPtr,ConstCxxPtr}, b::Union{CxxPtr,ConstCxxPtr}) = (a.cpp_object == b.cpp_object)
Base.:(==)(a::CxxBaseRef, b::Ptr) = (a.cpp_object == b)
Base.:(==)(a::Ptr, b::CxxBaseRef) = (b == a)
Base.:(==)(a::Union{CxxRef,ConstCxxRef}, b) = (a[] == b)
Base.:(==)(a, b::Union{CxxRef,ConstCxxRef}) = (b == a)
Base.:(==)(a::Union{CxxRef,ConstCxxRef}, b::Union{CxxRef,ConstCxxRef}) = (a[] == b[])

_deref(p::CxxBaseRef, ::Type) = unsafe_load(p.cpp_object)
_deref(p::CxxBaseRef{T}, ::Type{IsCxxType}) where {T} = dereferenced_type(T)(p.cpp_object)

_store_to_cxxptr(::Any,::Any,::Type) = error("Resetting the value to a C++ pointer or reference is only supported for non-wrapped types")
_store_to_cxxptr(r::Union{ConstCxxPtr{T},ConstCxxRef{T}}, x::T, ::Type{IsNormalType}) where {T} = error("Setting the value of a const reference or pointer is not allowed")
_store_to_cxxptr(r::Union{CxxPtr{T},CxxRef{T}}, x::T, ::Type{IsNormalType}) where {T} = unsafe_store!(r.cpp_object, x)

Base.unsafe_load(p::CxxBaseRef{T}) where {T} = _deref(p, cpp_trait_type(T))
_julia_pointer(p::CxxBaseRef{T}) where {T} = reinterpret(Ptr{julia_int_type(T)}, p.cpp_object)
Base.unsafe_string(p::CxxBaseRef) = unsafe_string(_julia_pointer(p))
Base.unsafe_string(p::CxxBaseRef, len::Integer) = unsafe_string(_julia_pointer(p), len)
Base.getindex(r::CxxBaseRef) = unsafe_load(r)
Base.setindex!(r::CxxBaseRef{T}, x::T) where {T}  = _store_to_cxxptr(r, x, cpp_trait_type(T))

# Delegate iteration to the contained type
Base.iterate(x::CxxBaseRef) = Base.iterate(x[])
Base.iterate(x::CxxBaseRef, state) = Base.iterate(x[], state)
Base.length(x::CxxBaseRef) = Base.length(x[])
Base.size(x::CxxBaseRef, d) = Base.size(x[], d)

# Delegate broadcast operations to the contained type
Base.BroadcastStyle(::Type{<:CxxBaseRef{T}}) where {T} = Base.BroadcastStyle(T)
Base.axes(x::CxxBaseRef) = Base.axes(x[])
Base.broadcastable(x::CxxBaseRef) = Base.broadcastable(x[])

# Delegate indexing to the wrapped type
Base.getindex(x::CxxBaseRef, i::Int) = Base.getindex(x[], i)
Base.setindex!(x::CxxBaseRef, val, i::Int) = Base.setindex!(x[], val, i)

Base.convert(::Type{RT}, p::SmartPointer{T}) where {T, RT <: CxxBaseRef{T}} = p[]
Base.cconvert(::Type{RT}, p::SmartPointer{T}) where {T, RT <: CxxBaseRef{T}} = p[]
function Base.convert(::Type{T1}, p::SmartPointer{DerivedT}) where {BaseT,T1 <: BaseT, DerivedT <: BaseT}
  return cxxupcast(T1, p[])[]
end
Base.convert(to_type::Type{Any}, x::CxxWrapCore.SmartPointer{DerivedT}) where {DerivedT} = x
function Base.convert(to_type::Type{<:Ref{T1}}, p::T2) where {BaseT,DerivedT, T1 <: BaseT, T2 <: SmartPointer{DerivedT}}
  return to_type(convert(T1,p))
end

Base.unsafe_convert(to_type::Type{<:CxxBaseRef}, x) = to_type(x.cpp_object)

# This is defined on the C++ side for each wrapped type
cxxupcast(x) = cxxupcast(CxxRef(x))
cxxupcast(x::CxxRef) = error("No upcast for type $(supertype(typeof(x))). Did you specialize SuperType to enable automatic upcasting?")
function cxxupcast(::Type{T}, x) where {T}
  cxxupcast(T, cxxupcast(x))
end
cxxupcast(::Type{T}, x::CxxBaseRef{T}) where {T} = x

struct ConstArray{T,N} <: AbstractArray{T,N}
  ptr::ConstCxxPtr{T}
  size::NTuple{N,Int}
end

ConstArray(ptr::ConstCxxPtr{T}, args::Vararg{Int,N}) where {T,N} = ConstArray{T,N}(ptr, (args...,))

Base.IndexStyle(::ConstArray) = IndexLinear()
Base.size(arr::ConstArray) = arr.size
Base.getindex(arr::ConstArray, i::Integer) = unsafe_load(arr.ptr.cpp_object, i)

function __delete end
function delete(x)
  __delete(CxxPtr(x))
  x.cpp_object = C_NULL
end

# Encapsulate information about a function
mutable struct CppFunctionInfo
  name::Any
  argument_types::Array{Type,1}
  return_type::Type
  julia_return_type::Type
  function_pointer::Ptr{Cvoid}
  thunk_pointer::Ptr{Cvoid}
  override_module::Module
end

# Interpreted as a constructor for Julia  > 0.5
mutable struct ConstructorFname
  _type::DataType
end

# Interpreted as an operator call overload
mutable struct CallOpOverload
  _type::DataType
end

_parameter_name(p) = string(p)
_parameter_name(p::DataType) = _fully_qualified_name(p)

_type_name(m::Module) = string(nameof(m))
function _type_name(::Type{T}) where {T}
  if T isa UnionAll
    return string(T)
  end
  typename = string(nameof(T))
  if isempty(T.parameters)
    return typename
  end
  return "$typename{$(join(_parameter_name.(T.parameters),","))}"
end

function _fully_qualified_name(root, modules=[parentmodule(root)])
  parent = parentmodule(modules[end])
  if parent == modules[end] || parent == Main
    if root == modules[end] || modules[end] == Main
      pop!(modules)
      if isempty(modules)
        return _type_name(root)
      end
    end
    return join(reverse(nameof.(modules)), ".") * "." * _type_name(root)
  end
  push!(modules, parent)
  return _fully_qualified_name(root, modules)
end

# Provide globally unique names for each distinct method signature
_method_name_string(funcname::ConstructorFname, mod) = "(::$(_fully_qualified_name(funcname._type)))"
_method_name_string(funcname::CallOpOverload, mod) = "(callable::$(_fully_qualified_name(funcname._type)))"
_method_name_string(funcname::Symbol, mod) = "$(_fully_qualified_name(mod)).$funcname"

function Base.print(io::IO, f::CppFunctionInfo)
  print(io, "$(_method_name_string(f.name, f.override_module))($(join(_fully_qualified_name.(f.argument_types), ", ")))::$(_fully_qualified_name(f.julia_return_type))")
end

# Pointers to function and thunk
const FunctionPointers = Tuple{Ptr{Cvoid},Ptr{Cvoid},Bool}

# Store a unique map between methods and their pointer, filled whenever a method is created in a module
# This solves a problem with e.g. vectors of vectors ov vectors of... where it is impossible to predict
# how many times and in which module a method will be defined
# This map is used to update a per-module vector of pointers upon module initialization, so it doesn't slow
# down each function call
const __global_method_map = Dict{String, FunctionPointers}()

function _register_function_pointers(func, precompiling)
  fstring = string(func)
  fptrs = (func.function_pointer, func.thunk_pointer, precompiling)
  if haskey(__global_method_map, fstring)
    existing = __global_method_map[fstring]
    if existing[3] == precompiling
      error("Double registration for method $fstring")
    end
  end
  __global_method_map[fstring] = fptrs
  return (fstring, fptrs)
end

function _get_function_pointer(fstr)
  if !haskey(__global_method_map, fstr)
    error("Method $fstr is undefined, maybe you need to precompile the Julia module?")
  end
  return __global_method_map[fstr]
end

function initialize_cxx_lib()
  ccall((:initialize_cxxwrap,libcxxwrap_julia), Cvoid, (Any, Any), @__MODULE__, CppFunctionInfo)
end

# Must also be called during precompile
initialize_cxx_lib()

function __init__()
  checkversion()
  initialize_cxx_lib()
end

function has_cxx_module(mod::Module)
  r = ccall((:has_cxx_module,libcxxwrap_julia), Cuchar, (Any,), mod)
  return r != 0
end

function register_julia_module(mod::Module, fptr::Ptr{Cvoid})
  ccall((:register_julia_module,libcxxwrap_julia), Cvoid, (Any,Ptr{Cvoid}), mod, fptr)
end

function initialize_julia_module(mod::Module)
  if has_cxx_module(mod) # Happens when not precompiling
    return
  end
  fptr = Libdl.dlsym(Libdl.dlopen(mod.__cxxwrap_sopath, mod.__cxxwrap_flags), mod.__cxxwrap_wrapfunc)
  register_julia_module(mod, fptr)
  funcs = get_module_functions(mod)
  precompiling = false
  for func in funcs
    _register_function_pointers(func, precompiling)
  end
  for (fidx,fstr) in enumerate(mod.__cxxwrap_methodnames)
    mod.__cxxwrap_pointers[fidx] = _get_function_pointer(fstr)
  end
end

function get_module_functions(mod::Module)
  ccall((:get_module_functions,libcxxwrap_julia), Any, (Any,), mod)
end

function bind_constants(m::Module, symbols::Array, values::Array)
  ccall((:bind_module_constants,libcxxwrap_julia), Cvoid, (Any,Any,Any), m, symbols, values)
end

function box_types(mod::Module)
  ccall((:get_box_types,libcxxwrap_julia), Any, (Any,), mod)
end

"""
Protect a variable from garbage collection by adding it to the global array kept by CxxWrap
"""
function gcprotect(x)
  ccall((:gcprotect,libcxxwrap_julia), Cvoid, (Any,), x)
end

"""
Unprotect a variable from garbage collection by removing it from the global array kept by CxxWrap
"""
function gcunprotect(x)
  ccall((:gcunprotect,libcxxwrap_julia), Cvoid, (Any,), x)
end

process_fname(fn::Tuple{<:Any,Module}, julia_mod) = process_fname(fn[1])
function process_fname(fn::Tuple{Symbol,Module}, julia_mod)
  (fname, mod) = fn
  if mod != julia_mod # Adding a method to a function from another module
    return :($mod.$fname)
  end
  return fname # defining a new function in the wrapped module, or adding a method to it
end
process_fname(fn::ConstructorFname) = :(::$(Type{fn._type}))
function process_fname(fn::CallOpOverload)
  return :(arg1::$(fn._type))
end

make_func_declaration(fn, argmap, julia_mod) = :($(process_fname(fn, julia_mod))($(argmap...)))
function make_func_declaration(fn::Tuple{CallOpOverload,Module}, argmap, julia_mod)
  return :($(process_fname(fn, julia_mod))($((argmap[2:end])...)))
end

# By default, no argument overloading happens
argument_overloads(t::Type) = Type[]
@static if Int != Cint
  argument_overloads(t::Type{Cint}) = [Int]
end
@static if UInt != Cuint
  argument_overloads(t::Type{Cuint}) = [UInt, Int]
else
  argument_overloads(t::Type{Cuint}) = [Int]
end
argument_overloads(t::Type{Float64}) = [Int, Irrational]
function argument_overloads(t::Type{Array{AbstractString,1}})
  return [Array{String,1}]
end
argument_overloads(t::Type{Ptr{T}}) where {T <: Number} = [Array{T,1}]

"""
Create a Union containing the type and a smart pointer to any type derived from it
"""
function ptrunion(::Type{T}) where {T}
  result{T2 <: T} = Union{T2, SmartPointer{T2}}
  return result
end

smart_pointer_type(t::Type) = smart_pointer_type(cpp_trait_type(t), t)
smart_pointer_type(::Type{IsNormalType}, t::Type) = t
smart_pointer_type(::Type{IsCxxType}, x::Type{T}) where {T} = ptrunion(x)

function smart_pointer_type(::Type{<:SmartPointer{T}}) where {T}
  result{T2 <: T} = SmartPointer{T2}
  return result
end

map_julia_arg_type(t::Type) = Union{Base.invokelatest(smart_pointer_type,t),argument_overloads(t)...}
map_julia_arg_type(a::Type{StrictlyTypedNumber{T}}) where {T} = T
map_julia_arg_type(a::Type{StrictlyTypedNumber{CxxBool}}) = Union{Bool,CxxBool}
map_julia_arg_type(x::Type{CxxBool}) = Union{Bool,CxxBool}
map_julia_arg_type(x::Type{T}) where {T<:Integer} = map_julia_arg_type(x, Base.invokelatest(cpp_trait_type, T))
map_julia_arg_type(x::Type{<:Integer}, ::Type{IsNormalType}) = Integer
map_julia_arg_type(x::Type{<:Integer}, ::Type{IsCxxType}) = x

const PtrTypes{T} = Union{CxxPtr{T}, Array{T}, CxxRef{T}, Base.RefValue{T}, Ptr{T},T}
const ConstPtrTypes{T} = Union{Ref{T}, Array{T}}

map_julia_arg_type(t::Type{<:CxxBaseRef{T}}) where {T} = map_julia_arg_type(t, Base.invokelatest(cpp_trait_type, T))

map_julia_arg_type(t::Type{ConstCxxRef{T}}, ::Type{IsNormalType}) where {T} = Union{Ref{T},map_julia_arg_type(T)}
map_julia_arg_type(t::Type{ConstCxxPtr{T}}, ::Type{IsNormalType}) where {T} = Union{ConstPtrTypes{T},Ptr{Cvoid}}
map_julia_arg_type(t::Type{CxxRef{T}}, ::Type{IsNormalType}) where {T} = PtrTypes{T}
map_julia_arg_type(t::Type{CxxPtr{T}}, ::Type{IsNormalType}) where {T} = Union{PtrTypes{T},Ptr{Cvoid}}

map_julia_arg_type(t::Type{ConstCxxRef{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = Union{ConstPtrTypes{julia_int_type(T)},map_julia_arg_type(T)}
map_julia_arg_type(t::Type{ConstCxxPtr{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = Union{ConstPtrTypes{julia_int_type(T)},Ptr{Cvoid}}
map_julia_arg_type(t::Type{CxxRef{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = PtrTypes{julia_int_type(T)}
map_julia_arg_type(t::Type{CxxPtr{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = Union{PtrTypes{julia_int_type(T)},Ptr{Cvoid}}

map_julia_arg_type(t::Type{ConstCxxRef{T}}, ::Type{IsCxxType}) where {T} = Union{map_julia_arg_type(T),ConstCxxRef{<:T},CxxRef{<:T}}
map_julia_arg_type(t::Type{ConstCxxPtr{T}}, ::Type{IsCxxType}) where {T} = Union{CxxPtr{<:T},ConstCxxPtr{<:T}, Ptr{Cvoid}}
map_julia_arg_type(t::Type{CxxRef{T}}, ::Type{IsCxxType}) where {T} = Union{map_julia_arg_type(T),CxxRef{<:T}}
map_julia_arg_type(t::Type{CxxPtr{T}}, ::Type{IsCxxType}) where {T} = Union{CxxPtr{<:T},Ptr{Cvoid}}

map_julia_arg_type(t::Type{CxxPtr{CxxChar}}) = Union{PtrTypes{Cchar}, String}
map_julia_arg_type(t::Type{<:Array{T}}) where {T <: Union{CxxSigned,CxxUnsigned}} = Union{t, Array{julia_int_type(T)}}
map_julia_arg_type(t::Type{<:Array{Ptr{T}}}) where {T <: Union{CxxSigned,CxxUnsigned}} = Union{t, Array{Ptr{julia_int_type(T)}}}
map_julia_arg_type(t::Type{ConstCxxPtr{CxxChar}}) = Union{ConstPtrTypes{Cchar}, String}

# names excluded from julia type mapping
const __excluded_names = Set([
      :cxxupcast,
      :__cxxwrap_smartptr_dereference,
      :__cxxwrap_smartptr_construct_from_other,
      :__cxxwrap_smartptr_cast_to_base
])

function Base.cconvert(to_type::Type{<:CxxBaseRef{T}}, x) where {T}
  return cxxconvert(to_type, x, cpp_trait_type(T))
end
function Base.cconvert(to_type::Type{<:CxxBaseRef{T}}, x::Ptr{PT}) where {T<:Integer,PT<:Integer}
  @assert T == PT || julia_int_type(T) == PT
  return to_type(x)
end
@inline Base.cconvert(::Type{T}, v::T) where {T2,T <: CxxBaseRef{T2}} = v
Base.cconvert(to_type::Type{<:CxxBaseRef{T}}, v::CxxBaseRef{T}) where {T} = to_type(v.cpp_object)
Base.unsafe_convert(to_type::Type{<:CxxBaseRef{T}}, v::CxxBaseRef) where {T} = to_type(v.cpp_object)
Base.unsafe_convert(to_type::Type{<:CxxBaseRef}, v::Base.RefValue) = to_type(pointer_from_objref(v))

cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x, ::Type{IsNormalType}) where {T} = Ref{T}(convert(T,x))
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Base.RefValue{T}, ::Type{IsNormalType}) where {T} = x
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Ptr{T}, ::Type{IsNormalType}) where {T} = to_type(x)
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Ptr{Cvoid}, ::Type{IsNormalType}) where {T} = to_type(x)
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Union{Array,String}, ::Type{IsNormalType}) where {T} = to_type(pointer(x))
function cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x, ::Type{IsCxxType}) where {T}
  return convert(T,x)
end

function cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::CxxBaseRef, ::Type{IsCxxType}) where {T}
  return to_type(convert(T,x[]))
end

# Build the expression to wrap the given function
function build_function_expression(func::CppFunctionInfo, funcidx, julia_mod)
  # Arguments and types
  argtypes = func.argument_types
  argsymbols = map((i) -> Symbol(:arg,i[1]), enumerate(argtypes))

  map_c_arg_type(t::Type) = t
  map_c_arg_type(::Type{Array{T,1}}) where {T <: AbstractString} = Any
  map_c_arg_type(::Type{Type}) = Any
  map_c_arg_type(::Type{T}) where {T <: Tuple} = Any
  map_c_arg_type(::Type{ConstArray{T,N}}) where {T,N} = Any
  map_c_arg_type(::Type{T}) where {T<:Union{CxxSigned,CxxUnsigned}} = julia_int_type(T)

  # Builds the return type passed to ccall
  map_c_return_type(t) = t
  map_c_return_type(::Type{T}) where {T <: Tuple} = Any
  map_c_return_type(::Type{ConstArray{T,N}}) where {T,N} = Any
  map_c_return_type(::Type{T}) where {T<:Union{CxxSigned,CxxUnsigned}} = map_c_arg_type(T)

  # Build the types for the ccall argument list
  c_arg_types = map_c_arg_type.(func.argument_types)
  c_return_type = map_c_return_type(func.return_type)

  # Builds the return-type annotation for the Julia function
  map_julia_return_type(t) = t
  map_julia_return_type(::Type{T}) where {T<:Union{CxxSigned,CxxUnsigned}} = map_c_arg_type(T)
  map_julia_return_type(::Type{CxxBool}) = Bool

  # Build the final call expression
  call_exp = quote end
  if func.thunk_pointer == C_NULL
    push!(call_exp.args, :(ccall(__cxxwrap_pointers[$funcidx][1], $c_return_type, ($(c_arg_types...),), $(argsymbols...)))) # Direct pointer call
  else
    push!(call_exp.args, :(ccall(__cxxwrap_pointers[$funcidx][1], $c_return_type, (Ptr{Cvoid}, $(c_arg_types...)), __cxxwrap_pointers[$funcidx][2], $(argsymbols...)))) # use thunk (= std::function)
  end

  function map_julia_arg_type_named(fname, t)
    if fname ∈ __excluded_names
      return t
    end
    return map_julia_arg_type(t)
  end

  # Build an array of arg1::Type1... expressions
  function argmap(signature)
    result = Expr[]
    for (t, s) in zip(signature, argsymbols)
      push!(result, :($s::$(map_julia_arg_type_named(func.name, t))))
    end
    return result
  end

  function_expression = :($(make_func_declaration((func.name,func.override_module), argmap(argtypes), julia_mod))::$(map_julia_return_type(func.julia_return_type)) = $call_exp)
  return function_expression
end

function makereftype(::Type{T}, mod) where {T}
  basename = T.name.name
  refname = Symbol(basename,"Dereferenced")
  tmod = T.name.module
  params = (T.parameters...,)
  if isempty(params)
    # If the type is non-parametric, this should be hit only once
    @assert !isdefined(mod, refname)
    @assert mod == tmod
    return Core.eval(mod, :(struct $refname <: $T cpp_object::Ptr{Cvoid} end; $refname))
  end
  if !isdefined(tmod, refname)
    parameternames = (Symbol(:T,i) for i in 1:length(params))
    expr = :(struct $refname{$(parameternames...)} <: $basename{$(parameternames...)} cpp_object::Ptr{Cvoid} end)
    Core.eval(mod, expr)
  end
  expr = :($tmod.$refname{$(params...)})
  return Core.eval(mod, expr)
end

function wrap_reference_converters(julia_mod)
  boxtypes = box_types(julia_mod)
  for bt in boxtypes
    st = supertype(bt)
    Core.eval(julia_mod, :(@inline $(@__MODULE__).cpp_trait_type(::Type{<:$st}) = $(@__MODULE__).IsCxxType))
    Core.eval(julia_mod, :(Base.convert(t::Type{$st}, x::T) where {T <: $st} = $(cxxupcast)($st,x)))
    Core.eval(julia_mod, :($(@__MODULE__).allocated_type(::Type{$st}) = $bt))
    reftype = makereftype(st, julia_mod)
    Core.eval(julia_mod, :($(@__MODULE__).dereferenced_type(::Type{$st}) = $reftype))
    Core.eval(julia_mod, :(Base.convert(::Type{$st}, x::$bt) = x))
    Core.eval(julia_mod, :(Base.convert(::Type{$st}, x::$reftype) = x))
    Core.eval(julia_mod, :(Base.:(==)(a::Union{CxxRef,ConstCxxRef,$bt}, b::$reftype) = (a.cpp_object == b.cpp_object)))
    Core.eval(julia_mod, :(Base.:(==)(a::$reftype, b::Union{CxxRef,ConstCxxRef,$bt}) = (b == a)))
  end
end

# Wrap functions from the cpp module to the passed julia module
function wrap_functions(functions, julia_mod)
  @assert isempty(julia_mod.__cxxwrap_pointers)
  precompiling = true

  for func in functions
    (fstring,fptrs) = _register_function_pointers(func, precompiling)
    push!(julia_mod.__cxxwrap_methodnames, fstring)
    push!(julia_mod.__cxxwrap_pointers, fptrs)
    funcidx = length(julia_mod.__cxxwrap_pointers)

    Core.eval(julia_mod, build_function_expression(func, funcidx, julia_mod))
  end
end

# Place the types for the module with the name corresponding to the current module name in the current module
function wraptypes(jlmod)
  symbols = Any[]
  values = Any[]
  bind_constants(jlmod, symbols, values)
  for (sym,val) in zip(symbols, values)
    Core.eval(jlmod, :(const $sym = $val))
  end
end

function wrapfunctions(jlmod)
  module_functions = get_module_functions(jlmod)
  wrap_reference_converters(jlmod)
  wrap_functions(module_functions, jlmod)
end

function readmodule(so_path::AbstractString, funcname, m::Module, flags)
  if flags === nothing
    flags = Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND
  end
  Core.eval(m, :(const __cxxwrap_methodnames = String[]))
  Core.eval(m, :(const __cxxwrap_pointers = $(FunctionPointers)[]))
  Core.eval(m, :(const __cxxwrap_sopath = $so_path))
  Core.eval(m, :(const __cxxwrap_wrapfunc = $(QuoteNode(funcname))))
  Core.eval(m, :(const __cxxwrap_flags = $flags))
  fptr = Libdl.dlsym(Libdl.dlopen(so_path, flags), funcname)
  register_julia_module(m, fptr)
end

function wrapmodule(so_path::AbstractString, funcname, m::Module, flags)
  readmodule(so_path, funcname, m, flags)
  wraptypes(m)
  wrapfunctions(m)
end

"""
  @wrapmodule libraryfile [functionname]

Place the functions and types from the C++ lib into the module enclosing this macro call
Calls an entry point named `define_julia_module`, unless another name is specified as
the second argument.
"""
macro wrapmodule(libraryfile, register_func=:(:define_julia_module), flags=:(nothing))
  return :(wrapmodule($(esc(libraryfile)), $(esc(register_func)), $__module__, $(esc(flags))))
end

"""
  @readmodule libraryfile [functionname]

Read a C++ module and associate it with the Julia module enclosing the macro call.
"""
macro readmodule(libraryfile, register_func=:(:define_julia_module), flags=:(nothing))
  return :(readmodule($(esc(libraryfile)), $(esc(register_func)), $__module__, $(esc(flags))))
end

"""
  @wraptypes

Wrap the types defined in the C++ side of the enclosing module. Requires that
`@readmodule` was called first.
"""
macro wraptypes()
  return :(wraptypes($__module__))
end

"""
  @wrapfunctions

Wrap the functions defined in the C++ side of the enclosing module. Requires that
`@readmodule` and `@wraptypes` was called first.
"""
macro wrapfunctions()
  return :(wrapfunctions($__module__))
end

"""
  @initcxx

Initialize the C++ pointer tables in a precompiled module using CxxWrap. Must be called from within
`__init__` in the wrapped module
"""
macro initcxx()
  return :(initialize_julia_module($__module__))
end

struct SafeCFunction
  fptr::Ptr{Cvoid}
  return_type::Type
  argtypes::Array{Type,1}
end

macro safe_cfunction(f, rt, args)
  return esc(:($(@__MODULE__).SafeCFunction(@cfunction($f, $rt, $args), $rt, [$(args.args...)])))
end

isnull(x::CxxBaseRef) = (x.cpp_object == C_NULL)

reference_type_union(::Type{T}) where {T} = reference_type_union(T, Base.invokelatest(cpp_trait_type, T))
reference_type_union(::Type{T}, ::Type{IsNormalType}) where {T} = T
reference_type_union(::Type{T}, ::Type{IsCxxType}) where {T} = Union{T, CxxBaseRef{T}, SmartPointer{T}}
reference_type_union(t::TypeVar) = t

dereference_argument(x) = x
dereference_argument(x::CxxBaseRef) = x[]
dereference_argument(x::SmartPointer) = x[]

macro cxxdereference(f)
  fdict = MacroTools.splitdef(f)

  function maparg(a)
    (argname, argtype, slurp, default) = MacroTools.splitarg(a)
    return MacroTools.combinearg(argname, :($(@__MODULE__).reference_type_union($(argtype))), slurp, default)
  end

  # Adapt the signature
  fdict[:args] .= maparg.(fdict[:args])
  fdict[:kwargs] .= maparg.(fdict[:kwargs])

  # Dereference the arguments
  deref_expr = quote end
  for arg in vcat(fdict[:args], fdict[:kwargs])
    (argname, _, slurp, _) = MacroTools.splitarg(arg)
    if argname === nothing
      continue
    end
    if !slurp
      push!(deref_expr.args, :($argname = $(@__MODULE__).dereference_argument($argname)))
    else
      push!(deref_expr.args, :($argname = (($(@__MODULE__).dereference_argument.($argname))...,)))
    end
  end
  insert!(fdict[:body].args, 1, deref_expr)
  fdict[:body] = MacroTools.flatten(fdict[:body])

  # Reassemble the function
  return esc(MacroTools.combinedef(fdict))
end

export @cxxdereference

end

include("StdLib.jl")

using .CxxWrapCore
using .CxxWrapCore: CxxBaseRef, argument_overloads, SafeCFunction, reference_type_union, dereference_argument, prefix_path

export @wrapmodule, @readmodule, @wraptypes, @wrapfunctions, @safe_cfunction, @initcxx,
ConstCxxPtr, ConstCxxRef, CxxRef, CxxPtr,
CppEnum, ConstArray, CxxBool, CxxLong, CxxULong, CxxChar, CxxWchar, CxxUChar, CxxSignedChar, CxxLongLong, CxxULongLong,
ptrunion, gcprotect, gcunprotect, isnull

using .StdLib: StdVector, StdString, StdWString, StdValArray

export StdVector, StdString, StdWString, StdValArray, @cxxdereference

end # module


text:
using ZChop # For chopping small imaginary parts in ρ
using Distributed

"""
    (t, FI, QFI) = Eff_QFI_PD(Nj, Ntraj, Tfinal, dt; kwargs... )

Evaluate the continuous-time FI and QFI of a final strong measurement for the
estimation of the frequency ω with continuous photo-detection monitoring
of each half-spin particle affected by noise at an angle θ,
with efficiency η using SME (stochastic master equation).

The function returns a tuple `(t, FI, QFI)` containing the time vector and the
vectors containing the FI and average QFI

# Arguments

* `Nj`: number of spins
* `Ntraj`: number of trajectories for the SSE
* `Tfinal`: final time of evolution
* `dt`: timestep of the evolution
* `κ = 1`: the noise coupling
* `θ = 0`: noise angle (0 parallel, π/2 transverse)
* `ω = 0`: local value of the frequency
* `η = 1`: measurement efficiency
"""
function Eff_QFI_PD(Nj::Int64,          # Number of spins
    Ntraj::Int64,                       # Number of trajectories
    Tfinal::Float64,                    # Final time
    dt::Float64;                        # Time step
    κ = 1.,                             # Noise coupling
    θ = 0.,                             # Noise angle
    ω = 0.,                             # Frequency
    η = 1.)                             # Measurement efficiency

    Ntime = Int(floor(Tfinal/dt)) # Number of timesteps
    dimJ = Int(2^Nj)   # Dimension of the corresponding Hilbert space

    # Define an array of noise channels
    cj = sqrt(κ/2) *
        [(cos(θ) * σ_j(:z, j, Nj) + sin(θ) * σ_j(:x, j, Nj)) for j = 1:Nj]

    # ω is the parameter that we want to estimate
    H = ω * σ(:z, Nj) / 2       # Hamiltonian of the spin system
    dH = σ(:z, Nj) / 2          # Derivative of H wrt the parameter ω

    # Kraus-like operator, trajectory-independent part
    M0 = sparse(I -
                1im * H * dt - 0.5 * dt * sum([c'*c for c in cj]) )
    M1 = sqrt(η * dt) * cj

    # Derivative of the Kraus-like operator wrt to ω
    dM0 = - 1im * dH * dt
    dM1 = 0

    # Probability of detecting a single photon (the probability of
    # detecting more than one in dt is negligible)
    # p_PD = η * sum(tr(ρ cj'*cj)) * dt,
    # but each noise operator gives 1/2 * κ * id when squared
    pPD = 0.5 * Nj * η * κ * dt;

    # Initial state of the system
    ψ0 = ghz_state(Nj)
    ρ0 = ψ0 * ψ0'

    t = (1 : Ntime) * dt

    # Run evolution for each trajectory, and build up the average
    # for FI and final strong measurement QFI
    result = @distributed (+) for ktraj = 1 : Ntraj
        ρ = ρ0 # Assign initial state to each trajectory

        # Derivative of ρ wrt the parameter
        # Initial state does not depend on the paramter
        dρ = zero(ρ)
        τ = dρ

        FisherT = zero(t)
        QFisherT = zero(t)

        for jt=1:Ntime
            # Has the photon been detected?
            if (rand() < real(pPD)) # Detected
                # Choose randomly which channel detected the photon
                # (with equal probabiltiy)
                ch = rand(1:Nj)

                new_ρ = M1[ch] * ρ * M1[ch]' ;

                zchop!(new_ρ) # Round off elements smaller than 1e-14
                tr_ρ = real(tr(new_ρ));

                τ = (M1[ch] * (τ * M1[ch]' + ρ * dM1') + dM1 * ρ * M1[ch]') / tr_ρ

            else # Not detected
                new_ρ = M0 * ρ * M0' + (1 - η) * dt * sum([c * ρ * c' for c in cj])
                zchop!(new_ρ)

                tr_ρ = real(tr(new_ρ));

                τ = (M0 * (τ * M0' + ρ * dM0') + dM0 * ρ * M0' +
                       (1 - η)* dt * sum([c* τ * c' for c in cj]))/ tr_ρ;
            end

            zchop!(τ) # Round off elements smaller than 1e-14

            tr_τ = tr(τ)

            # Now we can renormalize ρ and its derivative wrt ω
            ρ = new_ρ / tr_ρ
            dρ = τ - tr_τ * ρ

            # We evaluate the classical FI for the continuous measurement
            FisherT[jt] = real(tr_τ^2)

            # We evaluate the QFI for a final strong measurement done at time t
            QFisherT[jt] = QFI(ρ, dρ)
        end

        # Use the reduction feature of @distributed for
        # (at the end of each cicle, sum the result to result)
        hcat(FisherT, QFisherT)
    end

    return (t, result[:,1] / Ntraj, result[:,2] / Ntraj)
end


text:
const ncclDatatype_t = Cint
datatype(::Type{Int8})    = ncclDatatype_t(0)
datatype(::Type{UInt8})   = ncclDatatype_t(1)
datatype(::Type{Int32})   = ncclDatatype_t(2)
datatype(::Type{UInt32})  = ncclDatatype_t(3)
datatype(::Type{Int64})   = ncclDatatype_t(4)
datatype(::Type{UInt64})  = ncclDatatype_t(5)
datatype(::Type{Float16}) = ncclDatatype_t(6)
datatype(::Type{Float32}) = ncclDatatype_t(7)
datatype(::Type{Float64}) = ncclDatatype_t(8)
datatype(T::DataType) = error("NCCL doesn't support datatype $T")

const ncclRedOp_t = Cint
op(::Type{typeof(+)})   = ncclRedOp_t(0)
op(::Type{typeof(*)})   = ncclRedOp_t(1)
op(::Type{typeof(min)}) = ncclRedOp_t(2)
op(::Type{typeof(max)}) = ncclRedOp_t(3)
op(T::DataType) = error("NCCL doesn't support reduction operator $T")

function allReduce(::F, send::CuArray{T}, recv::CuArray{T}, comm::Communicator, stream=CUDAdrv.CuDefaultStream()) where {F, T}
    @assert size(send) == size(recv)
    @apicall(:ncclAllReduce, (CuPtr{T}, CuPtr{T}, Csize_t, ncclDatatype_t, ncclRedOp_t, ncclComm_t, CuStream_t),
                             send, recv, length(send), datatype(T), op(F), comm, stream)
end


text:
module Inputs

include("PWscf.jl")

end


text:
using PLSRegressor
using Test
using MLJBase
using MLJ
using Statistics
using Random

include("./utils_test.jl")
include("./pls1_test.jl")
include("./pls2_test.jl")
include("./kpls_test.jl")


text:
using SparseArrays
using Distributions
using StatsBase
using GVFHordes
import GVFHordes: AbstractHorde
using MinimalRLCore

ExperienceReplayCuriosity(size, obs_size, obs_type, num_act) =
    DynaExperienceReplay(size,
                         (obs_type, Int, obs_type, Int, Float64, Float64, Float64, Float64, Bool),
                         (obs_size, 1, obs_size, 1, num_act, num_act, 1, 1, 1),
                         (:s, :a, :ns, :na, :actprob, :actprob_prime, :env_r, :rint, :t))

mutable struct AgentER{O, Φ, SC, R} <: AbstractAgent

    behaviour_learner::Learner
    behaviour_demons::Union{AbstractHorde,Nothing} # Round Robin learners have behaviour demons set to nothing
    behaviour_gamma::Float64

    demons::AbstractHorde
    demon_learner::Learner

    last_obs::O
    last_state::Φ

    last_action::Int
    num_actions::Int

    intrinsic_reward::IntrinsicReward
    state_constructor::SC
    use_external_reward::Bool
    exploration::ExplorationStrategy

    random_first_action::Bool

    batch_size::Int
    replay::R

end

function AgentER(horde,
                 behaviour_feature_size::Int,
                 behaviour_learner,
                 behaviour_horde,
                 behaviour_gamma,
                 demon_learner,
                 observation_size::Int,
                 num_actions::Int,
                 intrinsic_reward_type,
                 state_constructor,
                 use_external_reward,
                 exploration_strategy,
                 random_first_action,
                 batch_size,
                 replay_size)

    behaviour_weight_dims = (num_actions, behaviour_feature_size)

    intrinsic_reward = if intrinsic_reward_type == "weight_change"
        #TODO: The intrinsic reward is defined by how the components of the agent are put together. For example, an intrinsic reward
        # could be the model error, which would then require different components that are assembled in the agent
        # Not sure how the construction of intrinisic reward could be abstracted out of this constructor
        WeightChange(demon_learner)
    elseif intrinsic_reward_type == "no_reward"
        NoReward()
    else
        throw(ArgumentError("Not a valid intrinsic reward"))
    end

    replay = ExperienceReplayCuriosity(replay_size, observation_size, Float64, num_actions)

    AgentER(behaviour_learner,
            behaviour_horde,
            behaviour_gamma,

            horde,
            demon_learner,

            zeros(observation_size),
            spzeros(behaviour_feature_size),

            0, # last_action
            num_actions,

            intrinsic_reward,
            state_constructor,

            use_external_reward,
            exploration_strategy,
            random_first_action,
            batch_size,
            replay)

end

function proc_input(agent::AgentER, obs)
    return agent.state_constructor(obs)
end

function get_action(agent::AgentER, state, obs)
    action_probs = if agent.behaviour_learner isa OneDTMazeUtils.RoundRobinPolicy || agent.behaviour_learner isa TwoDGridWorldUtils.RoundRobinPolicy
        get_action_probs(agent.behaviour_learner, state, obs)
    elseif agent.behaviour_learner.update isa TabularRoundRobin
        get_action_probs(agent.behaviour_learner.update, state, obs)
    elseif agent.behaviour_learner isa BaselineUtils.FollowDemon
        agent.behaviour_learner(obs)
    elseif agent.behaviour_learner isa BaselineUtils.RandomDemons
        qs = agent.behaviour_learner(obs)
        agent.exploration(qs)
    else
        qs = agent.behaviour_learner(state)
        agent.exploration(qs)
    end
    action = StatsBase.sample(1:agent.num_actions, Weights(action_probs))
    return action, action_probs
end

function assign_horde!(agent::AgentER, horde)
    agent.demons = horde
end

MinimalRLCore.end!(agent::AgentER, obs, reward, is_terminal) =
    MinimalRLCore.step!(agent, obs, reward, is_terminal)

function MinimalRLCore.start!(agent::AgentER, obs, args...)
    next_state = proc_input(agent, obs)
    #Always exploring starts
    step!(agent.exploration)
    next_action = if agent.random_first_action
        sample(1:agent.num_actions, Weights(ones(agent.num_actions)))
    else
        a, probs = get_action(agent,next_state, obs)
        a
    end

    # NOTE: IS THIS RIGHT, seems wierd!?
    _, discounts, _ = get(agent.demons; state_t = obs, action_t = next_action, state_tp1 = obs, action_tp1 = next_action)

    agent.last_state = next_state
    agent.last_action = next_action
    agent.last_obs = obs
    zero_eligibility_traces!(agent.demon_learner)

    return next_action
end

function MinimalRLCore.step!(agent::AgentER, obs, r, is_terminal, args...)
    next_state = proc_input(agent, obs)
    step!(agent.exploration)
    next_action, next_action_probs = get_action(agent, next_state, obs)

    #=
    Dealing with replay
    =#
    let # (:s, :a, :ns, :na, :actprob, :r, :t)
        ap = get_behaviour_pis(agent, agent.last_state, agent.last_obs)
        ap_prime = get_behaviour_pis(agent, next_state, obs)
        r_int = update_reward!(agent.intrinsic_reward, agent)
        total_reward = agent.use_external_reward ? r_int + r : r_int

        exp = (agent.last_obs, agent.last_action, obs, next_action, ap, ap_prime, r, total_reward, is_terminal)
        # @show typeof(exp)
        # @show typeof(agent.replay)
        push!(agent.replay, exp)
    end

    update_demons!(agent.demon_learner,
                   agent,
                   agent.last_obs,
                   obs,
                   agent.last_state,
                   agent.last_action,
                   next_state,
                   next_action,
                   is_terminal,
                   r)

    #get intrinssic reward

    update_behaviour!(agent)
                      # agent.last_obs,
                      # obs,
                      # agent.last_state,
                      # agent.last_action,
                      # next_state,
                      # next_action,
                      # is_terminal,
                      # total_reward)

    agent.last_state = next_state
    agent.last_action = next_action
    agent.last_obs = obs
    return next_action
end

get_behaviour_pis(agent::AgentER, state, obs) =
    get_action(agent, state, obs)[2]

function μ_π(agent::AgentER, obs)
    state = proc_input(agent, obs)
    get_behaviour_pis(agent,state,obs)
end

function update_demons!(::QLearner, agent::AgentER, args...)#, obs, next_obs, state, action, next_state, next_action, is_terminal, env_reward)

    exp = StatsBase.sample(agent.replay, agent.batch_size)
    # (:s, :a, :ns, :na, :actprob, :actprob_prime, :env_r, :rint, :t))
    for traj in exp
        state = proc_input(agent, traj.s)
        next_state = proc_input(agent, traj.ns)
        update!(agent.demon_learner,
                agent.demons,
                traj.s,
                traj.ns,
                state,
                traj.a,
                next_state,
                traj.na,
                is_terminal,
                (s, obs) -> obs === traj.s ? traj.actprob : traj.actprob_prime, #get_behaviour_pis(agent, state, obs),
                traj.env_r)
    end

end

function update_demons!(::SRLearner, agent::AgentER, obs, next_obs, state, action, next_state, next_action, is_terminal, env_reward)

    # update reward part

    update_rew_part!(agent.demon_learner.update,
                     agent.demon_learner,
                     agent.demons,
                     obs,
                     next_obs,
                     state,
                     action,
                     next_state,
                     next_action,
                     is_terminal,
                     (state, obs) -> get_behaviour_pis(agent, state, obs),
                     env_reward)

    # update SR part
    exp = StatsBase.sample(agent.replay, agent.batch_size)
    # (:s, :a, :ns, :na, :actprob, :actprob_prime, :env_r, :rint, :t))
    for traj in exp
        state = proc_input(agent, traj.s)
        next_state = proc_input(agent, traj.ns)
        update_sr_part!(agent.demon_learner.update,
                        agent.demon_learner,
                        agent.demons,
                        traj.s,
                        traj.ns,
                        state,
                        traj.a,
                        next_state,
                        traj.na,
                        is_terminal,
                        (s, obs) -> obs === traj.s ? traj.actprob : traj.actprob_prime,
                        traj.env_r)
    end

end

function update_behaviour!(agent::AgentER) #, obs, next_obs, state, action, next_state, next_action, is_terminal, reward)

    # behaviour_pis = get_behaviour_pis(agent, state, obs)
    # next_behaviour_pis = get_behaviour_pis(agent, next_state, next_obs)

    #NOTE: Different call than demon updates as the reward and environment pseudotermination function
    #
    exp = StatsBase.sample(agent.replay, agent.batch_size)

    for traj in exp
        state = proc_input(agent, traj.s)
        next_state = proc_input(agent, traj.ns)
        update!(agent.behaviour_learner,
                agent.behaviour_demons,
                traj.s,
                traj.ns,
                state,
                traj.a,
                next_state,
                traj.na,
                is_terminal,
                (s, obs) -> obs === traj.s ? traj.actprob : traj.actprob_prime, #get_behaviour_pis(agent, state, obs),
                traj.rint)

    end
end

function get_demon_prediction(agent::AgentER, obs, action)
    ϕ = agent.state_constructor(obs)
    agent.demon_learner(ϕ, action)
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity