text:
# VSCode specific
# ---------------

nodocuemnt_error(uri, data=nothing) =
    return JSONRPC.JSONRPCError(-32099, "document $(uri) requested but not present in the JLS", data)

function mismatched_version_error(uri, doc, params, msg, data=nothing)
    return JSONRPC.JSONRPCError(
        -32099,
        "version mismatch in $(msg) request for $(uri): JLS $(doc._version), client: $(params.version)",
        data
    )
end

# lookup
# ------

traverse_by_name(f, cache = SymbolServer.stdlibs) = traverse_store!.(f, values(cache))

traverse_store!(_, _) = return
traverse_store!(f, store::SymbolServer.EnvStore) = traverse_store!.(f, values(store))
function traverse_store!(f, store::SymbolServer.ModuleStore)
    for (sym, val) in store.vals
        f(sym, val)
        traverse_store!(f, val)
    end
end

# misc
# ----

function uri2filepath(uri::AbstractString)
    parsed_uri = try
        URIParser.URI(uri)
    catch
        throw(LSUriConversionFailure("Cannot parse `$uri`."))
    end

    if parsed_uri.scheme !== "file"
        return nothing
    end

    path_unescaped = URIParser.unescape(parsed_uri.path)
    host_unescaped = URIParser.unescape(parsed_uri.host)

    value = ""

    if host_unescaped != "" && length(path_unescaped) > 1
        # unc path: file://shares/c$/far/boo
        value = "//$host_unescaped$path_unescaped"
    elseif length(path_unescaped) >= 3 &&
            path_unescaped[1] == '/' &&
            isascii(path_unescaped[2]) && isletter(path_unescaped[2]) &&
            path_unescaped[3] == ':'
        # windows drive letter: file:///c:/far/boo
        value = lowercase(path_unescaped[2]) * path_unescaped[3:end]
    else
        # other path
        value = path_unescaped
    end

    if Sys.iswindows()
        value = replace(value, '/' => '\\')
    end

    value = normpath(value)

    return value
end

function filepath2uri(file::String)
    isabspath(file) || throw(LSRelativePath("Relative path `$file` is not valid."))
    if Sys.iswindows()
        file = normpath(file)
        file = replace(file, "\\" => "/")
        file = URIParser.escape(file)
        file = replace(file, "%2F" => "/")
        if startswith(file, "//")
            # UNC path \\foo\bar\foobar
            return string("file://", file[3:end])
        else
            # windows drive letter path
            return string("file:///", file)
        end
    else
        file = normpath(file)
        file = URIParser.escape(file)
        file = replace(file, "%2F" => "/")
        return string("file://", file)
    end
end

function escape_uri(uri::AbstractString)
    if !startswith(uri, "file://") # escaping only file URI
        return uri
    end
    escaped_uri = uri[8:end] |> URIParser.unescape |> URIParser.escape
    return string("file://", replace(escaped_uri, "%2F" => "/"))
end

function should_file_be_linted(uri, server)
    !server.runlinter && return false

    uri_path = uri2filepath(uri)

    if length(server.workspaceFolders) == 0 || uri_path === nothing
        return false
    else
        return any(i -> startswith(uri_path, i), server.workspaceFolders)
    end
end

# CompletionItemKind(t) = t in [:String, :AbstractString] ? 1 :
#                                 t == :Function ? 3 :
#                                 t == :DataType ? 7 :
#                                 t == :Module ? 9 : 6

# SymbolKind(t) = t in [:String, :AbstractString] ? 15 :
#                         t == :Function ? 12 :
#                         t == :DataType ? 5 :
#                         t == :Module ? 2 :
#                         t == :Bool ? 17 : 13

# Find location of default datatype constructor
const DefaultTypeConstructorLoc = let def = first(methods(Int))
    Base.find_source_file(string(def.file)), def.line
end

# TODO I believe this will also remove files from documents that were added
# not because they are part of the workspace, but by either StaticLint or
# the include follow logic.
function remove_workspace_files(root, server)
    for (uri, doc) in getdocuments_pair(server)
        # We first check whether the doc still exists on the server
        # because a previous loop iteration could have deleted it
        # via dependency removal of files
        hasdocument(server, uri) || continue
        fpath = getpath(doc)
        isempty(fpath) && continue
        get_open_in_editor(doc) && continue
        # If the file is in any other workspace folder, don't delete it
        any(folder -> startswith(fpath, folder), server.workspaceFolders) && continue
            deletedocument!(server, uri)
        end
    end

function Base.getindex(server::LanguageServerInstance, r::Regex)
    out = []
    for (uri, doc) in getdocuments_pair(server)
        occursin(r, uri._uri) && push!(out, doc)
    end
    return out
end

function _offset_unitrange(r::UnitRange{Int}, first=true)
    return r.start - 1:r.stop
end

function get_toks(doc, offset)
    ts = CSTParser.Tokenize.tokenize(get_text(doc))
    ppt = CSTParser.Tokens.RawToken(CSTParser.Tokens.ERROR, (0, 0), (0, 0), 1, 0, CSTParser.Tokens.NO_ERR, false, false)
    pt = CSTParser.Tokens.RawToken(CSTParser.Tokens.ERROR, (0, 0), (0, 0), 1, 0, CSTParser.Tokens.NO_ERR, false, false)
    t = CSTParser.Tokenize.Lexers.next_token(ts)

    prevpos = -1 # TODO: remove.
    while t.kind != CSTParser.Tokenize.Tokens.ENDMARKER
        if t.startbyte === prevpos # TODO: remove.
            throw(LSInfiniteLoop("Loop did not progress between iterations.")) # TODO: remove.
        else # TODO: remove.
            prevpos = t.startbyte # TODO: remove.
        end # TODO: remove.

        if t.startbyte < offset <= t.endbyte + 1
            break
        end
        ppt = pt
        pt = t
        t = CSTParser.Tokenize.Lexers.next_token(ts)
    end
    return ppt, pt, t
end

function isvalidjlfile(path)
    endswith(path, ".jl")
end

function get_expr(x, offset, pos=0, ignorewhitespace=false)
    if pos > offset
        return nothing
    end
    if length(x) > 0 && headof(x) !== :NONSTDIDENTIFIER
        for a in x
            if pos < offset <= (pos + a.fullspan)
                return get_expr(a, offset, pos, ignorewhitespace)
            end
            pos += a.fullspan
        end
    elseif pos == 0
        return x
    elseif (pos < offset <= (pos + x.fullspan))
        ignorewhitespace && pos + x.span < offset && return nothing
        return x
    end
end

# like get_expr, but only returns a expr if offset is not on the edge of its span
function get_expr_or_parent(x, offset, pos=0)
    if pos > offset
        return nothing, pos
    end
    ppos = pos
    if length(x) > 0 && headof(x) !== :NONSTDIDENTIFIER
        for a in x
            if pos < offset <= (pos + a.fullspan)
                if pos < offset < (pos + a.span)
                    return get_expr_or_parent(a, offset, pos)
                else
                    return x, ppos
                end
            end
            pos += a.fullspan
        end
    elseif pos == 0
        return x, pos
    elseif (pos < offset <= (pos + x.fullspan))
        if pos + x.span < offset
            return x.parent, ppos
        end
        return x, pos
    end
    return nothing, pos
end

function get_expr(x, offset::UnitRange{Int}, pos=0, ignorewhitespace=false)
    if all(pos .> offset)
        return nothing
    end
    if length(x) > 0 && headof(x) !== :NONSTDIDENTIFIER
        for a in x
            if all(pos .< offset .<= (pos + a.fullspan))
                return get_expr(a, offset, pos, ignorewhitespace)
            end
            pos += a.fullspan
        end
    elseif pos == 0
        return x
    elseif all(pos .< offset .<= (pos + x.fullspan))
        ignorewhitespace && all(pos + x.span .< offset) && return nothing
        return x
    end
    pos -= x.fullspan
    if all(pos .< offset .<= (pos + x.fullspan))
        ignorewhitespace && all(pos + x.span .< offset) && return nothing
        return x
    end
end

function get_expr1(x, offset, pos=0)
    if length(x) == 0 || headof(x) === :NONSTDIDENTIFIER
        if pos <= offset <= pos + x.span
            return x
        else
            return nothing
        end
    else
        for i = 1:length(x)
            arg = x[i]
            if pos < offset < (pos + arg.span) # def within span
                return get_expr1(arg, offset, pos)
            elseif arg.span == arg.fullspan
                if offset == pos
                    if i == 1
                        return get_expr1(arg, offset, pos)
                    elseif headof(x[i - 1]) === :IDENTIFIER
                        return get_expr1(x[i - 1], offset, pos)
                    else
                        return get_expr1(arg, offset, pos)
                    end
                elseif i == length(x) # offset == pos + arg.fullspan
                    return get_expr1(arg, offset, pos)
                end
            else
                if offset == pos
                    if i == 1
                        return get_expr1(arg, offset, pos)
                    elseif headof(x[i - 1]) === :IDENTIFIER
                        return get_expr1(x[i - 1], offset, pos)
                    else
                        return get_expr1(arg, offset, pos)
                    end
                elseif offset == pos + arg.span
                    return get_expr1(arg, offset, pos)
                elseif offset == pos + arg.fullspan
                elseif pos + arg.span < offset < pos + arg.fullspan
                    return nothing
                end
            end
            pos += arg.fullspan
        end
        return nothing
    end
end

function get_identifier(x, offset, pos=0)
    if pos > offset
        return nothing
    end
    if length(x) > 0
        for a in x
            if pos <= offset <= (pos + a.span)
                return get_identifier(a, offset, pos)
            end
            pos += a.fullspan
        end
    elseif headof(x) === :IDENTIFIER && (pos <= offset <= (pos + x.span)) || pos == 0
        return x
    end
end

if VERSION < v"1.1" || Sys.iswindows() && VERSION < v"1.3"
    _splitdir_nodrive(path::String) = _splitdir_nodrive("", path)
    function _splitdir_nodrive(a::String, b::String)
        m = match(Base.Filesystem.path_dir_splitter, b)
        m === nothing && return (a, b)
        a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
        a, String(m.captures[3])
    end
    splitpath(p::AbstractString) = splitpath(String(p))

    function splitpath(p::String)
        drive, p = _splitdrive(p)
        out = String[]
        isempty(p) && (pushfirst!(out, p))  # "" means the current directory.
        while !isempty(p)
            dir, base = _splitdir_nodrive(p)
            dir == p && (pushfirst!(out, dir); break)  # Reached root node.
            if !isempty(base)  # Skip trailing '/' in basename
                pushfirst!(out, base)
            end
            p = dir
        end
        if !isempty(drive)  # Tack the drive back on to the first element.
            out[1] = drive * out[1]  # Note that length(out) is always >= 1.
        end
        return out
    end
    _path_separator    = "\\"
    _path_separator_re = r"[/\\]+"
    function _pathsep(paths::AbstractString...)
        for path in paths
            m = match(_path_separator_re, String(path))
            m !== nothing && return m.match[1:1]
        end
        return _path_separator
    end
    function joinpath(a::String, b::String)
        isabspath(b) && return b
        A, a = _splitdrive(a)
        B, b = _splitdrive(b)
        !isempty(B) && A != B && return string(B,b)
        C = isempty(B) ? A : B
        isempty(a)                              ? string(C,b) :
        occursin(_path_separator_re, a[end:end]) ? string(C,a,b) :
                                                  string(C,a,_pathsep(a,b),b)
    end
    joinpath(a::AbstractString, b::AbstractString) = joinpath(String(a), String(b))
    joinpath(a, b, c, paths...) = joinpath(joinpath(a, b), c, paths...)
end

@static if Sys.iswindows() && VERSION < v"1.3"
    function _splitdir_nodrive(a::String, b::String)
        m = match(r"^(.*?)([/\\]+)([^/\\]*)$", b)
        m === nothing && return (a, b)
        a = string(a, isempty(m.captures[1]) ? m.captures[2][1] : m.captures[1])
        a, String(m.captures[3])
    end
    function _dirname(path::String)
        m = match(r"^([^\\]+:|\\\\[^\\]+\\[^\\]+|\\\\\?\\UNC\\[^\\]+\\[^\\]+|\\\\\?\\[^\\]+:|)(.*)$"s, path)
        m === nothing && return ""
        a, b = String(m.captures[1]), String(m.captures[2])
        _splitdir_nodrive(a, b)[1]
    end
    function _splitdrive(path::String)
        m = match(r"^([^\\]+:|\\\\[^\\]+\\[^\\]+|\\\\\?\\UNC\\[^\\]+\\[^\\]+|\\\\\?\\[^\\]+:|)(.*)$"s, path)
        m === nothing && return "", path
        String(m.captures[1]), String(m.captures[2])
    end
    function _splitdir(path::String)
        a, b = _splitdrive(path)
        _splitdir_nodrive(a, b)
    end
else
    _dirname = dirname
    _splitdir = splitdir
    _splitdrive = splitdrive
end

function valid_id(s::String)
    !isempty(s) && all(i == 1 ? Base.is_id_start_char(c) : Base.is_id_char(c) for (i, c) in enumerate(s))
end

function sanitize_docstring(doc::String)
    doc = replace(doc, "```jldoctest" => "```julia")
    doc = replace(doc, "\n#" => "\n###")
    return doc
end

function parent_file(x::EXPR)
    if parentof(x) isa EXPR
        return parent_file(parentof(x))
    elseif parentof(x) === nothing && StaticLint.haserror(x) && StaticLint.errorof(x) isa Document
        return x.meta.error
    else
        return nothing
    end
end

function resolve_op_ref(x::EXPR, env)
    StaticLint.hasref(x) && return true
    !CSTParser.isoperator(x) && return false
    pf = parent_file(x)
    pf === nothing && return false
    scope = StaticLint.retrieve_scope(x)
    scope === nothing && return false

    return op_resolve_up_scopes(x, CSTParser.str_value(x), scope, env)
end

function op_resolve_up_scopes(x, mn, scope, env)
    scope isa StaticLint.Scope || return false
    if StaticLint.scopehasbinding(scope, mn)
        StaticLint.setref!(x, scope.names[mn])
        return true
    elseif scope.modules isa Dict && length(scope.modules) > 0
        for (_, m) in scope.modules
            if m isa SymbolServer.ModuleStore && StaticLint.isexportedby(Symbol(mn), m)
                StaticLint.setref!(x, StaticLint.maybe_lookup(m[Symbol(mn)], env))
                return true
            elseif m isa StaticLint.Scope && StaticLint.scopehasbinding(m, mn)
                StaticLint.setref!(x, StaticLint.maybe_lookup(m.names[mn], env))
                return true
            end
        end
    end
    CSTParser.defines_module(scope.expr) || !(StaticLint.parentof(scope) isa StaticLint.Scope) && return false
    return op_resolve_up_scopes(x, mn, StaticLint.parentof(scope), env)
end

function is_in_target_dir_of_package(pkgpath, target)
    try # Safe failure - attempts to read disc.
        spaths = splitpath(pkgpath)
        if (i = findfirst(==(target), spaths)) !== nothing && "src" in readdir(joinpath(spaths[1:i - 1]...))
            return true
        end
        return false
    catch
        return false
    end
end


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

using Agents
using InteractiveDynamics
using GLMakie
using LinearAlgebra

include("Utils.jl")
include("Parameters.jl")

import .Parameters as P
import .Utils as U

#=
FUNCTIONS FOR DRAWING ON THE PLOT
=#
draw_line!(pos1, pos2; kwargs...) = lines!([pos1[1], pos2[1]], [pos1[2], pos2[2]]; kwargs...)

function draw_road_boundaries!(road::HorizontalRoad)
    draw_line!(
        top_boundary(road)...;
        color=P.ROAD_BOUNDARY_COLOR,
        linewidth=P.ROAD_BOUNDARY_LW
    )
    draw_line!(
        bottom_boundary(road)...;
        color=P.ROAD_BOUNDARY_COLOR,
        linewidth=P.ROAD_BOUNDARY_LW
    )
    return nothing
end

function draw_cross_path!(cp::HcrossPath)
    topy = top_boundary(cp.road)[1][2]
    boty = bottom_boundary(cp.road)[1][2]
    hw = P.CROSS_PATH_WIDTH /2
    line1 = ((cp.xpos - hw, topy), (cp.xpos - hw, boty))
    line2 = ((cp.xpos + hw, topy), (cp.xpos + hw, boty))
    for line in [line1, line2]
        draw_line!(
            line...;
            color=P.CP_COLOR,
            linewidth=P.CP_LW
        )
    end
end

function draw_cross_path!(cp::VcrossPath)
    leftx = left_boundary(cp.road)[1][1]
    rightx = right_boundary(cp.road)[1][1]
    hw = P.CROSS_PATH_WIDTH /2
    line1 = ((leftx, cp.ypos - hw), (rightx, cp.ypos - hw))
    line2 = ((leftx, cp.ypos + hw), (rightx, cp.ypos + hw))
    for line in [line1, line2]
        draw_line!(
            line...;
            color=P.CP_COLOR,
            linewidth=P.CP_LW
        )
    end
end

function draw_road_boundaries!(road::VerticalRoad)
    draw_line!(
        right_boundary(road)...;
        color=P.ROAD_BOUNDARY_COLOR,
        linewidth=P.ROAD_BOUNDARY_LW
    )
    draw_line!(
        left_boundary(road)...;
        color=P.ROAD_BOUNDARY_COLOR,
        linewidth=P.ROAD_BOUNDARY_LW
    )
    return nothing
end

function draw_signal!(road::Road)
    draw_signal!(road.signal)
    return nothing
end

function draw_signal!(road::Union{TwoWayHroad, TwoWayVroad, TwoWayIntersectingRoads})
    foreach(draw_signal!, signals(road))
    return nothing
end

function draw_signal!(s::Signal)
    s.active || return nothing
    if s.plot === missing
        signalPt = Point2f(s.pos...)
        plot = scatter!(signalPt, color=s.state, markersize=P.SIGNAL_MS)
        s.plot = plot
    else
        processSignalState!(s)
        s.plot.color = s.state
    end
    return nothing
end

function draw_lane_markers!(road::HorizontalRoad)
    last_boundary = top_boundary(road)
    for _ in 2:road.numLanes
        last_boundary = U.reduce_y(last_boundary, road.laneWidth)
        draw_line!(last_boundary..., color=P.LANE_MARKER_COLOR, linestyle=P.LANE_MARKER_LS)
    end
    return nothing
end

function draw_lane_markers!(road::VerticalRoad)
    last_boundary = right_boundary(road)
    for _ in 2:road.numLanes
        last_boundary = U.reduce_x(last_boundary, road.laneWidth)
        draw_line!(last_boundary..., color=P.LANE_MARKER_COLOR, linestyle=P.LANE_MARKER_LS)
    end
    return nothing
end

function draw_pedestrian_walkway!(road::HorizontalRoad)
    pedWay = nothing
    if road.startPos[1] < road.endPos[1]
        pedWay = U.increase_y(top_boundary(road), P.PEDESTRIAN_WALKWAY_WIDTH)
    else
        pedWay = U.reduce_y(bottom_boundary(road), P.PEDESTRIAN_WALKWAY_WIDTH)
    end
    draw_line!(pedWay..., color=P.PW_COLOR)
end

function draw_pedestrian_walkway!(road::VerticalRoad)
    pedWay = nothing
    if road.startPos[2] > road.endPos[2]
        pedWay = U.increase_x(right_boundary(road), P.PEDESTRIAN_WALKWAY_WIDTH)
    else
        pedWay = U.reduce_x(left_boundary(road), P.PEDESTRIAN_WALKWAY_WIDTH)
    end
    draw_line!(pedWay..., color=P.PW_COLOR)
end

function drawRoad!(road::Road)
    draw_road_boundaries!(road)
    draw_lane_markers!(road)
    draw_signal!(road)
    draw_pedestrian_walkway!(road)
    return nothing
end

function drawMedian!(road::TwoWayHroad)
    startPt = (startXPos(road), road.Ypos)
    endPt = (endXPos(road), road.Ypos)
    draw_line!(startPt, endPt, color=P.MEDIAN_COLOR, linewidth=road.medianWidth)
    return nothing
end

function drawMedian!(road::TwoWayVroad)
    startPt = (road.Xpos, startYPos(road))
    endPt = (road.Xpos, endYPos(road))
    draw_line!(startPt, endPt, color=P.MEDIAN_COLOR, linewidth=road.medianWidth)
    return nothing
end

function drawRoad!(road::TwoWayHroad)
    drawRoad!(road.L2Rroad)
    drawMedian!(road)
    drawRoad!(road.R2Lroad)
    return nothing
end

function drawRoad!(road::TwoWayVroad)
    drawRoad!(road.B2Troad)
    drawMedian!(road)
    drawRoad!(road.T2Broad)
    return nothing
end

function drawRoad!(interRoad::TwoWayIntersectingRoads)
    drawRoad!(interRoad.leftRoadSeg)
    drawRoad!(interRoad.topRoadSeg)
    drawRoad!(interRoad.rightRoadSeg)
    drawRoad!(interRoad.bottomRoadSeg)
    return nothing
end

ext: jl


text:
using ReachabilityAnalysis, JLD2
using ReachabilityAnalysis: add_dimension

LazySets.set_ztol(Float64, 1e-15)
ISS_path = joinpath(@__DIR__, "ISS.jld2")

@load ISS_path C
const C3 = C[3, :] # variable y₃
const C3_ext = vcat(C3, fill(0.0, 3))

function ISSF01()
    @load ISS_path A B

    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.]);  # input set
    X0 = BallInf(zeros(size(A, 1)), 0.0001)  # -0.0001 <= xi <= 0.0001 for all i
    prob_ISSF01 = @ivp(x' = A*x + B*u, x(0) ∈ X0, u ∈ U, x ∈ Universe(270))
end

function ISSC01()
    @load ISS_path A B

    Aext = add_dimension(A, 3)
    Aext[1:270, 271:273] = B

    U = Hyperrectangle(low=[0.0, 0.8, 0.9], high=[0.1, 1., 1.])  # input set
    X0 = BallInf(zeros(size(A, 1)), 0.0001)  # -0.0001 <= xi <= 0.0001 for all i
    X0 = X0 * U
    prob_ISSC01 = @ivp(x' = Aext*x, x(0) ∈ X0)
end


text:
"""
    ProjectFile(; version=v"0.1.0")

Creates a `Project.toml`.

## Keyword Arguments
- `version::VersionNumber`: The initial version of created packages.
"""
@with_kw_noshow struct ProjectFile <: Plugin
    version::VersionNumber = v"0.1.0"
end

# Other plugins like Tests will modify this file.
priority(::ProjectFile, ::typeof(hook)) = typemax(Int) - 5

function hook(p::ProjectFile, t::Template, pkg_dir::AbstractString)
    toml = Dict(
        "name" => basename(pkg_dir),
        "uuid" => string(uuid4()),
        "authors" => t.authors,
        "version" => string(p.version),
        "compat" => Dict("julia" => compat_version(t.julia)),
    )
    open(io -> TOML.print(io, toml), joinpath(pkg_dir, "Project.toml"), "w")
end

"""
    compat_version(v::VersionNumber) -> String

Format a `VersionNumber` to exclude trailing zero components.
"""
function compat_version(v::VersionNumber)
    return if v.patch == 0 && v.minor == 0
        "$(v.major)"
    elseif v.patch == 0
        "$(v.major).$(v.minor)"
    else
        "$(v.major).$(v.minor).$(v.patch)"
    end
end


text:
push!(LOAD_PATH,"./src/")

using Documenter, Slack

makedocs(sitename = "Slack.jl")


text:
using Documenter, Sads

makedocs(;
    modules=[Sads],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/kandread/Sads.jl/blob/{commit}{path}#L{line}",
    sitename="Sads.jl",
    authors="Kostas Andreadis",
    assets=[],
)

deploydocs(;
    repo="github.com/kandread/Sads.jl",
)


text:
import Polyhedra.fulldim
export EntropyCone, polymatroidcone, redundant, getinequalities, getextremerays, tight!

# Entropy Cone

abstract type AbstractEntropyCone{T<:Real} end

Polyhedra.fulldim(h::AbstractEntropyCone) = Polyhedra.fulldim(h.poly)

mutable struct EntropyCone{T<:Real} <: AbstractEntropyCone{T}
    n::Int
    poly::Polyhedron{T}

    function EntropyCone{T}(n::Int, p::Polyhedron{T}) where T
        if ntodim(n) != Polyhedra.fulldim(p)
            error("The number of variables does not match the dimension of the polyhedron")
        end
        new{T}(n, p)
    end
end
EntropyCone(n::Int, p::Polyhedron{T}) where {T} = EntropyCone{T}(n, p)

#EntropyCone{T<:AbstractFloat}(n::Int, A::AbstractMatrix{T}) = EntropyCone{Float64}(n, AbstractMatrix{Float64}(A), BitSet([]))
#EntropyCone{T<:Real}(n::Int, A::AbstractMatrix{T}) = EntropyCone{Rational{BigInt}}(n, AbstractMatrix{Rational{BigInt}}(A), BitSet())

#Base.getindex{T<:Real}(H::EntropyCone{T}, i) = DualEntropy(H.n, H.A[i,:], i in H.equalities) # FIXME

Base.copy(h::EntropyCone{T}) where {T<:Real} = EntropyCone{T}(h.n, copy(h.poly))

function indset(h::EntropyCone)
    indset(h.n)
end

function getinequalities(h::EntropyCone)
    removeredundantinequalities!(h.poly)
    ine = getinequalities(h.poly)
    if sum(abs(ine.b)) > 0
        error("Error: b is not zero-valued.")
    end
    [DualEntropy(ine.A[i,:], i in ine.linset) for i in 1:size(ine.A, 1)]
end

function getextremerays(h::EntropyCone)
    removeredundantgenerators!(h.poly)
    ext = SimpleVRepresentation(getgenerators(h.poly))
    if size(ext.V, 1) > 0
        error("Error: There are vertices.")
    end
    [PrimalEntropy(ext.R[i,:]) for i in 1:size(ext.R, 1)]
end

function Base.intersect!(H::AbstractEntropyCone, h::AbstractDualEntropy{L}) where {L}
    @assert Polyhedra.fulldim(H) == Polyhedra.fulldim(h)
    if H.n != h.n
        error("The dimension of the cone and entropy differ")
    end
    intersect!(H.poly, hrep(h))
end
function Base.intersect!(H::EntropyCone, hs::Vector{<:DualEntropy{L}}) where {L}
    @assert all(Polyhedra.fulldim(H) == Polyhedra.fulldim(h) for h in hs)
    intersect!(H.poly, hrep(hs))
end

function Base.intersect!(h::AbstractEntropyCone, hr::HRepresentation)
    h.poly = intersect(h.poly, hr)
end
function Base.intersect!(h1::AbstractEntropyCone, h2::AbstractEntropyCone)
    if h1.n != h2.n
        error("The dimension for the cones differ")
    end
    intersect!(h1.poly, h2.poly)
end

function Base.intersect(h1::AbstractEntropyCone, h2::AbstractEntropyCone)
    if h1.n != h2.n
        error("The dimension for the cones differ")
    end
    typeof(h1)(h1.n, intersect(h1.poly, h2.poly))
end

function polymatroidcone(::Type{T}, n::Integer, lib = nothing, minimal = true) where T
    # 2^n-1           nonnegative   inequalities H(S) >= 0
    # n*2^(n-1)-n     nondecreasing inequalities H(S) >= H(T) https://oeis.org/A058877
    # n*(n+1)*2^(n-2) submodular    inequalities              https://oeis.org/A001788

    # Actually, nonnegative is not required and nondecreasing only for H([n]) and H([n] \ i)
    n_nonnegative   = minimal ? 0 : 2^n-1
    n_nondecreasing = minimal ? n : n*2^(n-1)-n
    n_submodular    = 0
    if n >= 3
        n_submodular  = (n-1)*n*2^(n-3)
    elseif n == 2
        n_submodular  = 1
    end
    offset_nonnegative   = 0
    offset_nondecreasing = n_nonnegative
    offset_submodular    = n_nonnegative + n_nondecreasing
    cur_nonnegative   = 1
    cur_nondecreasing = 1
    cur_submodular    = 1
    HT = HalfSpace{T, SparseVector{T, Int}}
    hss = Vector{HT}(undef, n_nonnegative + n_nondecreasing + n_submodular)
    for j = 1:n
        for k = (j+1):n
            hss[offset_submodular+cur_submodular] = HRepElement(submodular(n, set(j), set(k)))
            cur_submodular += 1
        end
    end
    for I = indset(n)
        if !minimal
            hss[offset_nonnegative+cur_nonnegative] = HRepElement(nonnegative(n, I))
            cur_nonnegative += 1
        end
        for j = 1:n
            if !myin(j, I)
                if !minimal || card(I) == n-1
                    hss[offset_nondecreasing+cur_nondecreasing] = HRepElement(nondecreasing(n, I, set(j)))
                    cur_nondecreasing += 1
                end
                for k = (j+1):n
                    if !myin(k, I)
                        hss[offset_submodular+cur_submodular] = HRepElement(submodular(n, set(j), set(k), I))
                        cur_submodular += 1
                    end
                end
            end
        end
    end
    @assert cur_nonnegative == n_nonnegative+1
    @assert cur_nondecreasing == n_nondecreasing+1
    @assert cur_submodular == n_submodular+1
    h = hrep(hss)
    if lib === nothing
        p = polyhedron(h)
    else
        p = polyhedron(h, lib)
    end
    EntropyCone(n, p)
end
polymatroidcone(n::Integer, lib = nothing, minimal = true) = polymatroidcone(Int, n, lib, minimal)

function tight!(h::EntropyCone)
    tightness = DualEntropy{Int}[setequality(nondecreasing(h.n, setdiff(fullset(h.n), set(i)), set(i))) for i in 1:h.n]
    intersect!(h, tightness)
end


text:
export Metadata

# subtypes are expected to have a 'ref::API.LLVMMetadataRef' field
abstract type Metadata end

Base.unsafe_convert(::Type{API.LLVMMetadataRef}, md::Metadata) = md.ref

const metadata_kinds = Vector{Type}(fill(Nothing, typemax(API.LLVMMetadataKind)+1))
function identify(::Type{Metadata}, ref::API.LLVMMetadataRef)
    kind = API.LLVMGetMetadataKind(ref)
    typ = @inbounds metadata_kinds[kind+1]
    typ === Nothing && error("Unknown metadata kind $kind")
    return typ
end
function register(T::Type{<:Metadata}, kind::API.LLVMMetadataKind)
    metadata_kinds[kind+1] = T
end

function refcheck(::Type{T}, ref::API.LLVMMetadataRef) where T<:Metadata
    ref==C_NULL && throw(UndefRefError())
    if Base.JLOptions().debug_level >= 2
        T′ = identify(Metadata, ref)
        if T != T′
            error("invalid conversion of $T′ metadata reference to $T")
        end
    end
end

# Construct a concretely typed metadata object from an abstract metadata ref
function Metadata(ref::API.LLVMMetadataRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(Metadata, ref)
    return T(ref)
end

function Base.show(io::IO, md::Metadata)
    output = unsafe_message(API.LLVMExtraPrintMetadataToString(md))
    print(io, output)
end

## metadata as value

# this is for interfacing with (older) APIs that accept a Value*, not a Metadata*

@checked struct MetadataAsValue <: Value
    ref::API.LLVMValueRef
end
register(MetadataAsValue, API.LLVMMetadataAsValueValueKind)

# NOTE: this can be used to both pack e.g. metadata as a value, and to extract the
#       value from an ValueAsMetadata, so we don't type-assert narrowly here
Value(md::Metadata; ctx::Context) = Value(API.LLVMMetadataAsValue2(ctx, md))

Base.convert(T::Type{<:Value}, val::Metadata) = Value(val)::T

# NOTE: we can't do this automatically, as we can't query the context of metadata...
#       add wrappers to do so? would also simplify, e.g., `string(::MDString)`

## value as metadata

abstract type ValueAsMetadata <: Metadata end

@checked struct ConstantAsMetadata <: ValueAsMetadata
    ref::API.LLVMMetadataRef
end
register(ConstantAsMetadata, API.LLVMConstantAsMetadataMetadataKind)

@checked struct LocalAsMetadata <: ValueAsMetadata
    ref::API.LLVMMetadataRef
end
register(LocalAsMetadata, API.LLVMLocalAsMetadataMetadataKind)

# NOTE: this can be used to both pack e.g. constants as metadata, and to extract the
#       metadata from an MetadataAsValue, so we don't type-assert narrowly here
Metadata(val::Value) = Metadata(API.LLVMValueAsMetadata(val))

Base.convert(T::Type{<:Metadata}, val::Value) = Metadata(val)::T

## values

export MDString

@checked struct MDString <: Metadata
    ref::API.LLVMMetadataRef
end
register(MDString, API.LLVMMDStringMetadataKind)

MDString(val::String; ctx::Context) =
    MDString(API.LLVMMDStringInContext2(ctx, val, length(val)))

function Base.string(md::MDString)
    len = Ref{Cuint}()
    ptr = API.LLVMExtraGetMDString2(md, len)
    return unsafe_string(convert(Ptr{Int8}, ptr), len[])
end

## nodes

export MDNode, operands

abstract type MDNode <: Metadata end

function operands(md::MDNode)
    nops = API.LLVMExtraGetMDNodeNumOperands2(md)
    ops = Vector{API.LLVMMetadataRef}(undef, nops)
    API.LLVMExtraGetMDNodeOperands2(md, ops)
    return [op == C_NULL ? nothing : Metadata(op) for op in ops]
end

# TODO: setindex?
function replace_operand(md::MDNode, i, new::Metadata)
    API.LLVMReplaceMDNodeOperandWith(md, i-1, new)
end

## tuples

export MDTuple

@checked struct MDTuple <: MDNode
    ref::API.LLVMMetadataRef
end
register(MDTuple, API.LLVMMDTupleMetadataKind)

# MDTuples are commonly referred to as MDNodes, so keep that name
MDNode(mds::Vector{<:Metadata}; ctx::Context) =
    MDTuple(API.LLVMMDNodeInContext2(ctx, mds, length(mds)))
MDNode(vals::Vector; ctx::Context) =
    MDNode(convert(Vector{Metadata}, vals); ctx)

# we support passing `nothing`, but convert it to a non-exported `MDNull` instance
# so that we can keep everything as a subtype of `Metadata`
struct MDNull <: Metadata end
Base.convert(::Type{Metadata}, ::Nothing) = MDNull()
Base.unsafe_convert(::Type{API.LLVMMetadataRef}, md::MDNull) =
    convert(API.LLVMMetadataRef, C_NULL)


text:
"""
All is used instead of MIToS 1.0 "all" or "*", because it's possible to dispatch on it.
"""
struct All end

_get_function_name(str::String)::String = split(str,'.')[end]

"""
`get_n_words{T <: Union{ASCIIString, UTF8String}}(line::T, n::Int)`
It returns a `Vector{T}` with the first `n` (possibles) words/fields (delimited
by space or tab). If there is more than `n` words, the last word
returned contains the finals words and the delimiters. The length of the
returned vector is `n` or less (if the number of words is less than `n`).
This is used for parsing the Stockholm format.

```jldoctest
julia> using MIToS.Utils

julia> get_n_words("#=GR O31698/18-71 SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH", 3)
3-element Array{String,1}:
 "#=GR"
 "O31698/18-71"
 "SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH"

```
"""
function get_n_words(line::String, n::Int)
    if isempty(line)
        return String[]
    end
    words = Array{String}(undef, n)
    N  = 1
    last_spaces = 0:0
    while true
        if N == n
            @inbounds words[N] = line[(last(last_spaces)+1):end]
            break
        end
        spaces = something(findnext(r"[ |\t]+", line, last(last_spaces)+1), 0:-1)
        if first(spaces) == 0
            @inbounds words[N] = line[(last(last_spaces)+1):end]
            break
        end
        @inbounds words[N] = line[(last(last_spaces)+1):(first(spaces)-1)]
        last_spaces = spaces
        N += 1
    end
    if N != n
      resize!(words, N)
    end
    words
end

"""
`hascoordinates(id)`
It returns `true` if `id`/sequence name has the format: **UniProt/start-end**
(i.e. O83071/192-246)
"""
function hascoordinates(id)
    occursin(r"^\w+/\d+-\d+$", id)
end

"""
Selects the first element of the vector. This is useful for unpacking one element vectors.
Throws a warning if there are more elements. `element_name` is *element* by default,
but the name can be changed using the second argument.
"""
function select_element(vector::Array{T,1}, element_name::String="element") where T
    len = length(vector)
    if len == 0
        throw(ErrorException("There is not $element_name"))
    elseif len != 1
        warn("There are more than one ($len) $element_name using the first.")
    end
    @inbounds return(vector[1])
end

"""
Returns a vector with the `part` ("upper" or "lower") of the square matrix `mat`.
The `diagonal` is not included by default.
"""
function matrix2list(mat::AbstractMatrix{T}; part="upper", diagonal::Bool=false) where T
    nrow, ncol = size(mat)
    if nrow != ncol
        throw(ErrorException("Should be a square matrix"))
    end
    if diagonal
        d = 0
        N = div((ncol * ncol) + ncol, 2)
    else
        d = 1
        N = div((ncol * ncol) - ncol, 2)
    end
    list = Array{T}(undef, N)
    k = 1
    if part=="upper"
        for i in 1:(ncol-d)
            for j in (i+d):ncol
                list[k] = mat[i, j]
                k += 1
            end
        end
    elseif part=="lower"
        for j in 1:(ncol-d)
            for i in (j+d):ncol
                list[k] = mat[i, j]
                k += 1
            end
        end
    else
        throw(ErrorException("part should be \"upper\" or \"lower\""))
    end
    list
end

"""
Returns a square symmetric matrix from the vector `vec`. `side` is the number of
rows/columns. The `diagonal` is not included by default, set to `true` if there are
diagonal elements in the list.
"""
function list2matrix(vec::AbstractVector{T}, side::Int; diagonal::Bool=false) where T
    d = diagonal ? 0 : 1
    mat = zeros(T, side, side)
    k = 1
    for i in 1:(side-d)
        for j in (i+d):side
            value = vec[k]
            mat[i, j] = value
            mat[j, i] = value
            k += 1
        end
    end
    mat
end

"""
It checks if a PDB code has the correct format.
"""
check_pdbcode(pdbcode::String) = occursin(r"^\w{4}$", pdbcode)

"""
Getter for the `array` field of `NamedArray`s
"""
getarray(x::NamedArray) = x.array


text:
#==============================================================================#
# AWSCore.jl
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

__precompile__()

module AWSCore

export AWSException, AWSConfig, aws_config, AWSRequest, post_request, do_request

using Retry
using SymDict
using XMLDict

using Compat
import Compat: String

typealias AWSConfig SymbolDict
typealias AWSRequest SymbolDict

include("http.jl")
include("AWSException.jl")
include("AWSCredentials.jl")
include("names.jl")
include("mime.jl")

#------------------------------------------------------------------------------#
# Configuration.
#------------------------------------------------------------------------------#

function aws_config(;creds=AWSCredentials(),
                     region=get(ENV, "AWS_DEFAULT_REGION", "us-east-1"),
                     args...)
    @SymDict(creds, region, args...)
end

#------------------------------------------------------------------------------#
# AWSRequest to Request.jl conversion.
#------------------------------------------------------------------------------#

# Construct a HTTP POST request dictionary for "servce" and "query"...
#
# e.g.
# aws = Dict(:creds  => AWSCredentials(),
#            :region => "ap-southeast-2")
#
# post_request(aws, "sdb", "2009-04-15", StrDict("Action" => "ListDomains"))
#
# Dict{Symbol, Any}(
#     :creds    => creds::AWSCredentials
#     :verb     => "POST"
#     :url      => "http://sdb.ap-southeast-2.amazonaws.com/"
#     :headers  => Dict("Content-Type" =>
#                       "application/x-www-form-urlencoded; charset=utf-8)
#     :content  => "Version=2009-04-15&ContentType=JSON&Action=ListDomains"
#     :resource => "/"
#     :region   => "ap-southeast-2"
#     :service  => "sdb"
# )

function post_request(aws::AWSConfig,
                      service::String,
                      version::String,
                      query::Dict)

    resource = get(aws, :resource, "/")
    url = aws_endpoint(service, aws[:region]) * resource
    if version != ""
        query["Version"] = version
    end
    headers = Dict("Content-Type" =>
                   "application/x-www-form-urlencoded; charset=utf-8")
    content = format_query_str(query)

    @SymDict(verb = "POST", service, resource, url, headers, query, content,
             aws...)
end

# Convert AWSRequest dictionary into Requests.Request (Requests.jl)

function Request(r::AWSRequest)
    Request(r[:verb], r[:resource], r[:headers], r[:content], URI(r[:url]))
end

# Call http_request for AWSRequest.

function http_request(request::AWSRequest, args...)
    http_request(Request(request), get(request, :return_stream, false))
end

# Pretty-print AWSRequest dictionary.

function dump_aws_request(r::AWSRequest)

    action = r[:verb]
    name = r[:resource]
    if haskey(r, :query) && haskey(r[:query], "Action")
        action = r[:query]["Action"]
    end
    if haskey(r, :query)
        for k in keys(r[:query])
            if ismatch(r"[^.]Name$", k)
                name *= " "
                name *= r[:query][k]
            end
        end
    end
    println("$(r[:service]).$action $name")
end

#------------------------------------------------------------------------------#
# AWSRequest retry loop
#------------------------------------------------------------------------------#

include("sign.jl")

function do_request(r::AWSRequest)

    response = nothing

    # Try request 3 times to deal with possible Redirect and ExiredToken...
    @repeat 3 try

        # Default headers...
        if !haskey(r, :headers)
            r[:headers] = Dict()
        end
        r[:headers]["User-Agent"] = "JuliaAWS.jl/0.0.0"
        r[:headers]["Host"]       = URI(r[:url]).host

        # Load local system credentials if needed...
        if !haskey(r, :creds) || r[:creds].token == "ExpiredToken"
            r[:creds] = AWSCredentials()
        end

        # Use credentials to sign request...
        sign!(r)

        if debug_level > 0
            dump_aws_request(r)
        end

        # Send the request...
        response = http_request(r)

    catch e

        # Handle HTTP Redirect...
        @retry if http_status(e) in [301, 302, 307] &&
                  haskey(headers(e), "Location")
            r[:url] = headers(e)["Location"]
        end

        e = AWSException(e)

        if debug_level > 0
            println("Warning: AWSCore.do_request() exception: $(typeof(e))")
        end

        # Handle expired signature...
        @retry if ismatch(r"Signature expired", e.message) end

        # Handle ExpiredToken...
        @retry if typeof(e) == ExpiredToken
            r[:creds].token = "ExpiredToken"
        end
    end

    # If there is no reponse data, return raw response object...
    if typeof(response) != Response || length(response.data) < 1
        return response
    end

    # Return raw data if there is no mimetype...
    if isnull(mimetype(response))
        return response.data
    end

    # Parse response data according to mimetype...
    mime = get(mimetype(response))

    if ismatch(r"/xml$", mime)
        return parse_xml(String(response.data))
    end

    if ismatch(r"/x-amz-json-1.0$", mime)
        return JSON.parse(String(response.data))
    end

    if ismatch(r"json$", mime)
        info = JSON.parse(String(response.data))
        @protected try
            action = r[:query]["Action"]
            info = info[action * "Response"]
            info = info[action * "Result"]
        catch e
            @ignore if typeof(e) == KeyError end
        end
        return info
    end

    if ismatch(r"^text/", mime)
        return String(response.data)
    end

    # Return raw data by default...
    return response.data
end

global debug_level = 0

function set_debug_level(n)
    global debug_level = n
end

end # module AWSCore

#==============================================================================#
# End of file.
#==============================================================================#


text:
using Test;

using AdalmPluto;

@testset "libIIO_jl" begin
    include("test_toplevel.jl");
    include("test_scan.jl");
    include("test_context.jl");
    include("test_device.jl");
    include("test_channel.jl");
    include("test_buffer.jl");
    # #  include("test_debug.jl");
end

@testset "AdalmPluto" begin
    #  openPluto
    #  openPluto (txCfg::ChannelCfg, …)
    global pluto = openPluto(Int64(960e6), Int64(2.8e6), 64; bufferSize=UInt64(1024), bandwidth=Int(25e6));
    # Functions called by openPluto :
    #   - updateGain! (pluto::PlutoSDR, …)
    #   - updateGainMode! (pluto::PlutoSDR, …)
    #   - updateEffectiveCfg! (trx::Union, …)

    # Verification of the correct initalization of the fields
    ## PlutoSDR.PlutoTx.txWrapper
    @test pluto.tx.iio.tx     != C_NULL;
    @test pluto.tx.iio.tx0_i  != C_NULL;
    @test pluto.tx.iio.tx0_q  != C_NULL;
    @test pluto.tx.iio.chn    != C_NULL;
    @test pluto.tx.iio.chn_lo != C_NULL;
    ## PlutoSDR.PlutoRx.rxWrapper
    @test pluto.rx.iio.rx     != C_NULL;
    @test pluto.rx.iio.rx0_i  != C_NULL;
    @test pluto.rx.iio.rx0_q  != C_NULL;
    @test pluto.rx.iio.chn    != C_NULL;
    @test pluto.rx.iio.chn_lo != C_NULL;

    ## PlutoSDR.PlutoTx.ChannelCfg
    @test pluto.tx.cfg.rfport       == "A";
    @test pluto.tx.cfg.bandwidth    == 25000000;
    @test pluto.tx.cfg.samplingRate == Int64(2.8e6);
    @test pluto.tx.cfg.carrierFreq  == Int64(960e6);
    ## PlutoSDR.PlutoRx.ChannelCfg
    @test pluto.rx.cfg.rfport       == "A_BALANCED";
    @test pluto.rx.cfg.bandwidth    == 25000000;
    @test pluto.rx.cfg.samplingRate == Int64(2.8e6);
    @test pluto.rx.cfg.carrierFreq  == Int64(960e6);

    ## PlutoSDR.PlutoTx.IIO_Buffer
    @test pluto.tx.buf.C_ptr         != C_NULL;
    @test pluto.tx.buf.C_size        == 1024;
    @test pluto.tx.buf.C_sample_size > 0;
    @test pluto.tx.buf.C_first       != C_NULL;
    @test pluto.tx.buf.C_last        != C_NULL;
    @test pluto.tx.buf.C_step        > 0;
    ## PlutoSDR.PlutoTx.IIO_Buffer
    @test pluto.rx.buf.C_ptr         != C_NULL;
    @test pluto.rx.buf.C_size        == 1024;
    @test pluto.rx.buf.C_sample_size > 0;
    @test pluto.rx.buf.C_first       != C_NULL;
    @test pluto.rx.buf.C_last        != C_NULL;
    @test pluto.rx.buf.C_step        > 0;

    ## PlutoSDR.PlutoTx
    @test isapprox(pluto.tx.effectiveSamplingRate, pluto.tx.cfg.samplingRate; atol=5);
    @test isapprox(pluto.tx.effectiveCarrierFreq, pluto.tx.cfg.carrierFreq; atol=5);
    ## PlutoSDR.PlutoRx
    @test isapprox(pluto.rx.effectiveSamplingRate, pluto.rx.cfg.samplingRate; atol=5);
    @test isapprox(pluto.rx.effectiveCarrierFreq, pluto.rx.cfg.carrierFreq; atol=5);

    #  updateBandwidth! (pluto, …)
    @test updateBandwidth!(pluto, Int64(20e6); doLog=false) == 0;
    @test pluto.tx.cfg.bandwidth == Int64(20e6);
    @test C_iio_channel_attr_read_longlong(pluto.tx.iio.chn, "rf_bandwidth") == (0, Int64(20e6));
    @test pluto.rx.cfg.bandwidth == Int64(20e6);
    @test C_iio_channel_attr_read_longlong(pluto.rx.iio.chn, "rf_bandwidth") == (0, Int64(20e6));

    #  updateSamplingRate! (pluto::PlutoSDR, …)
    @test updateSamplingRate!(pluto, Int64(5e6); doLog=false) == 0;
    # Tx
    @test pluto.tx.cfg.samplingRate == Int64(5e6);
    @test isapprox(pluto.tx.effectiveSamplingRate, pluto.tx.cfg.samplingRate; atol=5);
    # Tx
    @test pluto.rx.cfg.samplingRate == Int64(5e6);
    @test isapprox(pluto.rx.effectiveSamplingRate, pluto.rx.cfg.samplingRate; atol=5);

    #  updateCarrierFreq! (pluto::PlutoSDR, …)
    @test updateCarrierFreq!(pluto, Int64(1050.5e6); doLog=false) == 0;
    # Tx
    @test pluto.tx.cfg.carrierFreq == Int64(1050.5e6);
    @test isapprox(pluto.tx.effectiveCarrierFreq, pluto.tx.cfg.carrierFreq; atol=5);
    # Rx
    @test pluto.rx.cfg.carrierFreq == Int64(1050.5e6);
    @test isapprox(pluto.rx.effectiveCarrierFreq, pluto.rx.cfg.carrierFreq; atol=5);

    sig = zeros(ComplexF32, 3000);
    #  recv (pluto::PlutoSDR, …)
    #  recv! (sig::Array, …)
    @test recv!(sig, pluto) == length(sig);
    @test sig[1] != 0 && sig[end] != 0; # very unlikely to have those equal 0
    @test pluto.rx.buf.nb_samples == 3072 - length(sig);

    # Testing gain
    gain = getGain(pluto)
    @test gain == 64

    # Send
    buffer = zeros(Complex{Float32},1024)
    send(pluto,buffer)

    #  Base.close (pluto::PlutoSDR, …)
    close(pluto);
    @test pluto.released = true;
    @test pluto.tx.released = true;
    @test pluto.rx.released = true;
end


text:
#
# Copyright 2013 - 2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
# Copyright 2015 Oliver Heimlich (oheim@posteo.de)
#
# Original author: Marco Nehmeier (unit tests in libieeep1788)
# Converted into portable ITL format by Oliver Heimlich with minor corrections.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
#Language imports

#Test library imports
using FactCheck

#Arithmetic library imports
using ValidatedNumerics

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :narrow)

facts("minimal_inf_test") do
    @fact infimum(∅) --> Inf
    @fact infimum(Interval(-Inf, Inf)) --> -Inf
    @fact infimum(Interval(1.0, 2.0)) --> 1.0
    @fact infimum(Interval(-3.0, -2.0)) --> -3.0
    @fact infimum(Interval(-Inf, 2.0)) --> -Inf
    @fact infimum(Interval(-Inf, 0.0)) --> -Inf
    @fact infimum(Interval(-Inf, -0.0)) --> -Inf
    @fact infimum(Interval(-2.0, Inf)) --> -2.0
    @fact infimum(Interval(0.0, Inf)) --> -0.0
    @fact infimum(Interval(-0.0, Inf)) --> -0.0
    @fact infimum(Interval(-0.0, 0.0)) --> -0.0
    @fact infimum(Interval(0.0, -0.0)) --> -0.0
    @fact infimum(Interval(0.0, 0.0)) --> -0.0
    @fact infimum(Interval(-0.0, -0.0)) --> -0.0
end

facts("minimal_inf_dec_test") do
    @fact infimum(DecoratedInterval(∅, trv)) --> Inf
    @fact infimum(DecoratedInterval(Interval(-Inf, Inf), def)) --> -Inf
    @fact infimum(DecoratedInterval(Interval(1.0, 2.0), com)) --> 1.0
    @fact infimum(DecoratedInterval(Interval(-3.0, -2.0), trv)) --> -3.0
    @fact infimum(DecoratedInterval(Interval(-Inf, 2.0), dac)) --> -Inf
    @fact infimum(DecoratedInterval(Interval(-Inf, 0.0), def)) --> -Inf
    @fact infimum(DecoratedInterval(Interval(-Inf, -0.0), trv)) --> -Inf
    @fact infimum(DecoratedInterval(Interval(-2.0, Inf), trv)) --> -2.0
    @fact infimum(DecoratedInterval(Interval(0.0, Inf), def)) --> -0.0
    @fact infimum(DecoratedInterval(Interval(-0.0, Inf), trv)) --> -0.0
    @fact infimum(DecoratedInterval(Interval(-0.0, 0.0), dac)) --> -0.0
    @fact infimum(DecoratedInterval(Interval(0.0, -0.0), trv)) --> -0.0
    @fact infimum(DecoratedInterval(Interval(0.0, 0.0), trv)) --> -0.0
    @fact infimum(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> -0.0
end

facts("minimal_sup_test") do
    @fact supremum(∅) --> -Inf
    @fact supremum(Interval(-Inf, Inf)) --> Inf
    @fact supremum(Interval(1.0, 2.0)) --> 2.0
    @fact supremum(Interval(-3.0, -2.0)) --> -2.0
    @fact supremum(Interval(-Inf, 2.0)) --> 2.0
    @fact supremum(Interval(-Inf, 0.0)) --> 0.0
    @fact supremum(Interval(-Inf, -0.0)) --> 0.0
    @fact supremum(Interval(-2.0, Inf)) --> Inf
    @fact supremum(Interval(0.0, Inf)) --> Inf
    @fact supremum(Interval(-0.0, Inf)) --> Inf
    @fact supremum(Interval(-0.0, 0.0)) --> 0.0
    @fact supremum(Interval(0.0, -0.0)) --> 0.0
    @fact supremum(Interval(0.0, 0.0)) --> 0.0
    @fact supremum(Interval(-0.0, -0.0)) --> 0.0
end

facts("minimal_sup_dec_test") do
    @fact supremum(DecoratedInterval(∅, trv)) --> -Inf
    @fact supremum(DecoratedInterval(Interval(-Inf, Inf), def)) --> Inf
    @fact supremum(DecoratedInterval(Interval(1.0, 2.0), com)) --> 2.0
    @fact supremum(DecoratedInterval(Interval(-3.0, -2.0), trv)) --> -2.0
    @fact supremum(DecoratedInterval(Interval(-Inf, 2.0), dac)) --> 2.0
    @fact supremum(DecoratedInterval(Interval(-Inf, 0.0), def)) --> 0.0
    @fact supremum(DecoratedInterval(Interval(-Inf, -0.0), trv)) --> 0.0
    @fact supremum(DecoratedInterval(Interval(-2.0, Inf), trv)) --> Inf
    @fact supremum(DecoratedInterval(Interval(0.0, Inf), def)) --> Inf
    @fact supremum(DecoratedInterval(Interval(-0.0, Inf), trv)) --> Inf
    @fact supremum(DecoratedInterval(Interval(-0.0, 0.0), dac)) --> +0.0
    @fact supremum(DecoratedInterval(Interval(0.0, -0.0), trv)) --> +0.0
    @fact supremum(DecoratedInterval(Interval(0.0, 0.0), trv)) --> +0.0
    @fact supremum(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> +0.0
end

facts("minimal_mid_test") do
    @fact mid(Interval(-Inf, Inf)) --> 0.0
    @fact mid(Interval(-0x1.fffffffffffffp1023, +0x1.fffffffffffffp1023)) --> 0.0
    @fact mid(Interval(0.0, 2.0)) --> 1.0
    @fact mid(Interval(2.0, 2.0)) --> 2.0
    @fact mid(Interval(-2.0, 2.0)) --> 0.0
    @fact mid(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022)) --> 0.0
    @fact mid(Interval(-0x0.0000000000001p-1022, 0x0.0000000000002p-1022)) --> 0.0
end

facts("minimal_mid_dec_test") do
    @fact mid(DecoratedInterval(Interval(-Inf, Inf), def)) --> 0.0
    @fact mid(DecoratedInterval(Interval(-0x1.fffffffffffffp1023, +0x1.fffffffffffffp1023), trv)) --> 0.0
    @fact mid(DecoratedInterval(Interval(0.0, 2.0), com)) --> 1.0
    @fact mid(DecoratedInterval(Interval(2.0, 2.0), dac)) --> 2.0
    @fact mid(DecoratedInterval(Interval(-2.0, 2.0), trv)) --> 0.0
    @fact mid(DecoratedInterval(Interval(-0x0.0000000000002p-1022, 0x0.0000000000001p-1022), trv)) --> 0.0
    @fact mid(DecoratedInterval(Interval(-0x0.0000000000001p-1022, 0x0.0000000000002p-1022), trv)) --> 0.0
end

facts("minimal_rad_test") do
    @fact radius(Interval(0.0, 2.0)) --> 1.0
    @fact radius(Interval(2.0, 2.0)) --> 0.0
    @fact radius(Interval(-Inf, Inf)) --> Inf
    @fact radius(Interval(0.0, Inf)) --> Inf
    @fact radius(Interval(-Inf, 1.2)) --> Inf
end

facts("minimal_rad_dec_test") do
    @fact radius(DecoratedInterval(Interval(0.0, 2.0), trv)) --> 1.0
    @fact radius(DecoratedInterval(Interval(2.0, 2.0), com)) --> 0.0
    @fact radius(DecoratedInterval(Interval(-Inf, Inf), trv)) --> Inf
    @fact radius(DecoratedInterval(Interval(0.0, Inf), def)) --> Inf
    @fact radius(DecoratedInterval(Interval(-Inf, 1.2), trv)) --> Inf
end

facts("minimal_wid_test") do
    @fact diam(Interval(2.0, 2.0)) --> 0.0
    @fact diam(Interval(1.0, 2.0)) --> 1.0
    @fact diam(Interval(1.0, Inf)) --> Inf
    @fact diam(Interval(-Inf, 2.0)) --> Inf
    @fact diam(Interval(-Inf, Inf)) --> Inf
end

facts("minimal_wid_dec_test") do
    @fact diam(DecoratedInterval(Interval(2.0, 2.0), com)) --> 0.0
    @fact diam(DecoratedInterval(Interval(1.0, 2.0), trv)) --> 1.0
    @fact diam(DecoratedInterval(Interval(1.0, Inf), trv)) --> Inf
    @fact diam(DecoratedInterval(Interval(-Inf, 2.0), def)) --> Inf
    @fact diam(DecoratedInterval(Interval(-Inf, Inf), trv)) --> Inf
end

facts("minimal_mag_test") do
    @fact mag(Interval(1.0, 2.0)) --> 2.0
    @fact mag(Interval(-4.0, 2.0)) --> 4.0
    @fact mag(Interval(-Inf, 2.0)) --> Inf
    @fact mag(Interval(1.0, Inf)) --> Inf
    @fact mag(Interval(-Inf, Inf)) --> Inf
    @fact mag(Interval(-0.0, 0.0)) --> 0.0
    @fact mag(Interval(-0.0, -0.0)) --> 0.0
end

facts("minimal_mag_dec_test") do
    @fact mag(DecoratedInterval(Interval(1.0, 2.0), com)) --> 2.0
    @fact mag(DecoratedInterval(Interval(-4.0, 2.0), trv)) --> 4.0
    @fact mag(DecoratedInterval(Interval(-Inf, 2.0), trv)) --> Inf
    @fact mag(DecoratedInterval(Interval(1.0, Inf), def)) --> Inf
    @fact mag(DecoratedInterval(Interval(-Inf, Inf), trv)) --> Inf
    @fact mag(DecoratedInterval(Interval(-0.0, 0.0), trv)) --> 0.0
    @fact mag(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> 0.0
end

facts("minimal_mig_test") do
    @fact mig(Interval(1.0, 2.0)) --> 1.0
    @fact mig(Interval(-4.0, 2.0)) --> 0.0
    @fact mig(Interval(-4.0, -2.0)) --> 2.0
    @fact mig(Interval(-Inf, 2.0)) --> 0.0
    @fact mig(Interval(-Inf, -2.0)) --> 2.0
    @fact mig(Interval(-1.0, Inf)) --> 0.0
    @fact mig(Interval(1.0, Inf)) --> 1.0
    @fact mig(Interval(-Inf, Inf)) --> 0.0
    @fact mig(Interval(-0.0, 0.0)) --> 0.0
    @fact mig(Interval(-0.0, -0.0)) --> 0.0
end

facts("minimal_mig_dec_test") do
    @fact mig(DecoratedInterval(Interval(1.0, 2.0), com)) --> 1.0
    @fact mig(DecoratedInterval(Interval(-4.0, 2.0), trv)) --> 0.0
    @fact mig(DecoratedInterval(Interval(-4.0, -2.0), trv)) --> 2.0
    @fact mig(DecoratedInterval(Interval(-Inf, 2.0), def)) --> 0.0
    @fact mig(DecoratedInterval(Interval(-Inf, -2.0), trv)) --> 2.0
    @fact mig(DecoratedInterval(Interval(-1.0, Inf), trv)) --> 0.0
    @fact mig(DecoratedInterval(Interval(1.0, Inf), trv)) --> 1.0
    @fact mig(DecoratedInterval(Interval(-Inf, Inf), trv)) --> 0.0
    @fact mig(DecoratedInterval(Interval(-0.0, 0.0), trv)) --> 0.0
    @fact mig(DecoratedInterval(Interval(-0.0, -0.0), trv)) --> 0.0
end
# FactCheck.exitstatus()


text:
import Test

Test.@testset "PredictMDFull.jl" begin
    import InteractiveUtils # stdlib
    import Pkg # stdlib
    import Test # stdlib

    @debug(string("Julia depot paths: "), Base.DEPOT_PATH)
    @debug(string("Julia load paths: "), Base.LOAD_PATH)

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        @debug(string("Julia version debug: ",))
        InteractiveUtils.versioninfo(verbose=true)
    end

    @debug(string("Attempting to import PredictMDFull...",))
    import PredictMDFull
    @debug(string("Successfully imported PredictMDFull.",))
    @debug(string("PredictMDFull version: "),PredictMDFull.version(),)
    @debug(string("PredictMDFull package directory: "),PredictMDFull.package_directory(),)

    PredictMDFull.import_all()

    @debug(string("Julia depot paths: "), Base.DEPOT_PATH)
    @debug(string("Julia load paths: "), Base.LOAD_PATH)

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        if ispath(Base.active_project())
            println(logger.stream, "# Location of test environment Project.toml: \"$(Base.active_project())\"")
            println(logger.stream, "# Beginning of test environment Project.toml")
            println(logger.stream, read(Base.active_project(), String))
            println(logger.stream, "# End of test environment Project.toml")
        else
            println(logger.stream, "# File \"$(Base.active_project())\" does not exist")
        end
        if ispath(joinpath(dirname(Base.active_project()), "Manifest.toml"))
            println(logger.stream, "# Location of test environment Manifest.toml: \"$(joinpath(dirname(Base.active_project()), "Manifest.toml"))\"")
            println(logger.stream, "# Beginning of test environment Manifest.toml")
            println(logger.stream, read(joinpath(dirname(Base.active_project()), "Manifest.toml"),String))
            println(logger.stream, "# End of test environment Manifest.toml")
        else
            println(logger.stream, "# File \"$(joinpath(dirname(Base.active_project()), "Manifest.toml"))\" does not exist")
        end
    end

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        if ispath(PredictMDFull.package_directory("Project.toml"))
            println(logger.stream, "# Location of PredictMDFull package Project.toml: \"$(PredictMDFull.package_directory("Project.toml"))\"")
            println(logger.stream, "# Beginning of PredictMDFull package Project.toml")
            println(logger.stream, read(PredictMDFull.package_directory("Project.toml"), String))
            println(logger.stream, "# End of PredictMDFull package Project.toml")
        else
            println(logger.stream, "# File \"$(PredictMDFull.package_directory("Project.toml"))\" does not exist")
        end
        if ispath(PredictMDFull.package_directory("Manifest.toml"))
            println(logger.stream, "# Location of PredictMDFull package Manifest.toml: \"$(PredictMDFull.package_directory("Manifest.toml"))\"")
            println(logger.stream, "# Beginning of PredictMDFull package Manifest.toml")
            println(logger.stream, read(PredictMDFull.package_directory("Manifest.toml"),String))
            println(logger.stream, "# End of PredictMDFull package Manifest.toml")
        else
            println(logger.stream, "# File \"$(PredictMDFull.package_directory("Manifest.toml"))\" does not exist")
        end
    end

    Test.@testset "Unit tests           " begin
        testmodulea_filename::String = joinpath("TestModuleA", "TestModuleA.jl")
        testmoduleb_filename::String  = joinpath(
            "TestModuleB", "directory1", "directory2", "directory3",
            "directory4", "directory5", "TestModuleB.jl",
            )
        testmodulec_filename::String  = joinpath(mktempdir(), "TestModuleC.jl")
        rm(testmodulec_filename; force = true, recursive = true)
        open(testmodulec_filename, "w") do io
            write(io, "module TestModuleC end")
        end
        include(testmodulea_filename)
        include(testmoduleb_filename)
        include(testmodulec_filename)
        include(joinpath("test_package_directory.jl"))
        include(joinpath("test_package_list.jl"))
        include(joinpath("test_registry_url_list.jl"))
        include(joinpath("test_version.jl"))
    end

    Test.@testset "Import PredictMDExtra" begin
        import PredictMDExtra
    end

    Test.@testset "Import PredictMD     " begin
        import PredictMD
    end

    Test.@testset "Test import_all()" begin
        include(joinpath("test_import_all.jl"))
    end

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        if ispath(Base.active_project())
            println(logger.stream, "# Location of test environment Project.toml: \"$(Base.active_project())\"")
            println(logger.stream, "# Beginning of test environment Project.toml")
            println(logger.stream, read(Base.active_project(), String))
            println(logger.stream, "# End of test environment Project.toml")
        else
            println(logger.stream, "# File \"$(Base.active_project())\" does not exist")
        end
        if ispath(joinpath(dirname(Base.active_project()), "Manifest.toml"))
            println(logger.stream, "# Location of test environment Manifest.toml: \"$(joinpath(dirname(Base.active_project()), "Manifest.toml"))\"")
            println(logger.stream, "# Beginning of test environment Manifest.toml")
            println(logger.stream, read(joinpath(dirname(Base.active_project()), "Manifest.toml"),String))
            println(logger.stream, "# End of test environment Manifest.toml")
        else
            println(logger.stream, "# File \"$(joinpath(dirname(Base.active_project()), "Manifest.toml"))\" does not exist")
        end
    end

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        if ispath(PredictMD.package_directory("Project.toml"))
            println(logger.stream, "# Location of PredictMD package Project.toml: \"$(PredictMD.package_directory("Project.toml"))\"")
            println(logger.stream, "# Beginning of PredictMD package Project.toml")
            println(logger.stream, read(PredictMD.package_directory("Project.toml"), String))
            println(logger.stream, "# End of PredictMD package Project.toml")
        else
            println(logger.stream, "# File \"$(PredictMD.package_directory("Project.toml"))\" does not exist")
        end
        if ispath(PredictMD.package_directory("Manifest.toml"))
            println(logger.stream, "# Location of PredictMD package Manifest.toml: \"$(PredictMD.package_directory("Manifest.toml"))\"")
            println(logger.stream, "# Beginning of PredictMD package Manifest.toml")
            println(logger.stream, read(PredictMD.package_directory("Manifest.toml"),String))
            println(logger.stream, "# End of PredictMD package Manifest.toml")
        else
            println(logger.stream, "# File \"$(PredictMD.package_directory("Manifest.toml"))\" does not exist")
        end
    end

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        if ispath(PredictMDExtra.package_directory("Project.toml"))
            println(logger.stream, "# Location of PredictMDExtra package Project.toml: \"$(PredictMDExtra.package_directory("Project.toml"))\"")
            println(logger.stream, "# Beginning of PredictMDExtra package Project.toml")
            println(logger.stream, read(PredictMDExtra.package_directory("Project.toml"), String))
            println(logger.stream, "# End of PredictMDExtra package Project.toml")
        else
            println(logger.stream, "# File \"$(PredictMDExtra.package_directory("Project.toml"))\" does not exist")
        end
        if ispath(PredictMDExtra.package_directory("Manifest.toml"))
            println(logger.stream, "# Location of PredictMDExtra package Manifest.toml: \"$(PredictMDExtra.package_directory("Manifest.toml"))\"")
            println(logger.stream, "# Beginning of PredictMDExtra package Manifest.toml")
            println(logger.stream, read(PredictMDExtra.package_directory("Manifest.toml"),String))
            println(logger.stream, "# End of PredictMDExtra package Manifest.toml")
        else
            println(logger.stream, "# File \"$(PredictMDExtra.package_directory("Manifest.toml"))\" does not exist")
        end
    end

    logger = Base.CoreLogging.current_logger_for_env(Base.CoreLogging.Debug, Symbol(splitext(basename(something(@__FILE__, "nothing")))[1]), something(@__MODULE__, "nothing"))
    if !isnothing(logger)
        if ispath(PredictMDFull.package_directory("Project.toml"))
            println(logger.stream, "# Location of PredictMDFull package Project.toml: \"$(PredictMDFull.package_directory("Project.toml"))\"")
            println(logger.stream, "# Beginning of PredictMDFull package Project.toml")
            println(logger.stream, read(PredictMDFull.package_directory("Project.toml"), String))
            println(logger.stream, "# End of PredictMDFull package Project.toml")
        else
            println(logger.stream, "# File \"$(PredictMDFull.package_directory("Project.toml"))\" does not exist")
        end
        if ispath(PredictMDFull.package_directory("Manifest.toml"))
            println(logger.stream, "# Location of PredictMDFull package Manifest.toml: \"$(PredictMDFull.package_directory("Manifest.toml"))\"")
            println(logger.stream, "# Beginning of PredictMDFull package Manifest.toml")
            println(logger.stream, read(PredictMDFull.package_directory("Manifest.toml"),String))
            println(logger.stream, "# End of PredictMDFull package Manifest.toml")
        else
            println(logger.stream, "# File \"$(PredictMDFull.package_directory("Manifest.toml"))\" does not exist")
        end
    end
end


text:
@doc raw"""

    struct VelocityParameters{T <: SSDFloat}

Values needed to parametrize the longitudinal drift velocity of electrons or hole
along a crystal axis as a function of the electric field strength.

## Background information

The parameterization for the longitudinal drift velocity, ``v_l``, as a function of the electric
field strength, ``E``, was proposed by [D.M. Caughey and R.E. Thomas](https://ieeexplore.ieee.org/document/1448053)
and later expanded by [L. Mihailescu et al.](https://www.sciencedirect.com/science/article/pii/S0168900299012863):
```math
v_l = \frac{\mu_0 E}{(1 + (E/E_0 )^{\beta})^{1/ \beta}} - \mu_{n} E.
```
with the four parameters, ``\mu_0``, ``\beta``, ``E_0`` and ``\mu_n``, which are different
for electrons and holes and for the different crystal axes.

!!! note
    The parameter ``\mu_n`` accounts for the Gunn effects for electrons and should be 0 for holes.

## Fields
* `mu0::T`: Parameter ``\mu_0`` in the parameterization shown above.
* `beta::T`: Parameter ``\beta`` in the parameterization shown above.
* `E0::T`: Parameter ``E_0`` in the parameterization shown above.
* `mun::T`: Parameter ``\mu_n`` in the parameterization shown above.

"""
struct VelocityParameters{T <: SSDFloat}
    mu0::T
    beta::T
    E0::T
    mun::T
end

@doc raw"""

    struct CarrierParameters{T <: SSDFloat}

Parameters needed to describe the electron or hole drift along the <100> and <111> axes.

## Fields
* `axis100::VelocityParameters{T}`: Parameters to describe the charge drift along the <100> axis.
* `axis111::VelocityParameters{T}`: Parameters to describe the charge drift along the <111> axis.

See also [`VelocityParameters`](@ref).
"""
struct CarrierParameters{T <: SSDFloat}
    axis100::VelocityParameters{T}
    axis111::VelocityParameters{T}
end

struct ADLParameters{T <: SSDFloat}
    ml_inv::T
    mt_inv::T
    Γ0::T
    Γ1::T
    Γ2::T
end

# Temperature models
include("TemperatureModels/TemperatureModels.jl")

# Electron model parametrization from [3]
@fastmath function γj(j::Integer, crystal_orientation::SMatrix{3,3,T,9}, γ0::SMatrix{3,3,T,9}, ::Type{HPGe})::SMatrix{3,3,T,9} where {T <: SSDFloat}
    tmp::T = 2 / 3
    a::T = acos(sqrt(tmp))
    Rx::SMatrix{3,3,T,9} = RotX{T}(a)
    b::T = (j - 1) * T(π) / 2 - π / 4
    Rzj::SMatrix{3,3,T,9} = RotZ{T}(-b)
    Rj = Rx * Rzj * crystal_orientation
    transpose(Rj) * γ0 * Rj
end

@fastmath function γj(j::Integer, crystal_orientation::SMatrix{3,3,T,9}, γ0::SMatrix{3,3,T,9}, ::Type{Si})::SMatrix{3,3,T,9} where {T <: SSDFloat}
    # needs to be updated!
    b::T = (j == 1 ? 0 : π / 2)
    Rzj::SMatrix{3,3,T,9} = (j == 3 ? SMatrix{3,3,T,9}(1, 0, 0, 0, 0, 1, 0, -1, 0)  : I) * SMatrix{3,3,T,9}(cos(b), -sin(b), 0, sin(b), cos(b), 0, 0, 0, 1)
    Rj = Rzj * crystal_orientation
    transpose(Rj) * γ0 * Rj
end

@fastmath function setup_γj(crystal_orientation::SMatrix{3,3,T,9}, p::ADLParameters{T}, material::Type{HPGe})::SVector{4, SMatrix{3,3,T,9}} where {T <: SSDFloat}
    γ0 = SMatrix{3,3,T,9}(p.mt_inv, 0, 0, 0, p.ml_inv, 0, 0, 0, p.mt_inv)
    SVector{4, SMatrix{3,3,T,9}}([γj(i, crystal_orientation, γ0, material) for i in Base.OneTo(4)]...)
end

@fastmath function setup_γj(crystal_orientation::SMatrix{3,3,T,9}, p::ADLParameters{T}, material::Type{Si})::SVector{3, SMatrix{3,3,T,9}} where {T <: SSDFloat}
    γ0 = SMatrix{3,3,T,9}(p.mt_inv, 0, 0, 0, p.ml_inv, 0, 0, 0, p.mt_inv)
    SVector{3, SMatrix{3,3,T,9}}([γj(i, crystal_orientation, γ0, material) for i in Base.OneTo(3)]...)
end

@fastmath function ADLParameters{T}(ml::T, mt::T, ::Type{HPGe})::ADLParameters{T} where {T}
    ml_inv::T = inv(ml)
    mt_inv::T = inv(mt)
    tmp::T = sqrt(8*mt_inv + ml_inv) / 3

    Γ0::T = inv(sqrt((2 * mt_inv + ml_inv) / 3))
    Γ1::T = (-4 * sqrt(ml_inv) - 4/3 * tmp) / (sqrt(ml_inv) - tmp)^2
    Γ2::T = Γ1 * (sqrt(ml_inv) + 3 * tmp) / (-4)
    ADLParameters{T}(ml_inv, mt_inv, Γ0, Γ1, Γ2)
end

@fastmath function ADLParameters{T}(ml::T, mt::T, ::Type{Si})::ADLParameters{T} where {T}
    ml_inv::T = inv(ml)
    mt_inv::T = inv(mt)

    Γ0::T = inv(sqrt((2 * mt_inv + ml_inv) / 3))
    Γ1::T = (-3 * sqrt(ml_inv) - 3/2 * sqrt(mt_inv) )/(sqrt(ml_inv) - sqrt(mt_inv))^2
    Γ2::T = Γ1 * (sqrt(ml_inv) + 2 * sqrt(mt_inv)) / (-3)
    ADLParameters{T}(ml_inv, mt_inv, Γ0, Γ1, Γ2)
end

# Longitudinal drift velocity formula
@fastmath function Vl(Emag::T, params::VelocityParameters{T})::T where {T <: SSDFloat}
    params.mu0 * Emag / (1 + (Emag / params.E0)^params.beta)^(1 / params.beta) - params.mun * Emag
end

@doc raw"""
    ADLChargeDriftModel{T <: SSDFloat, M <: AbstractDriftMaterial, N, TM <: AbstractTemperatureModel{T}} <: AbstractChargeDriftModel{T}

Charge drift model for electrons and holes based on the AGATA Detector Library.
Find a detailed description of the calculations in [ADL Charge Drift Model](@ref).

## Fields
- `electrons::CarrierParameters{T}`: Parameters to describe the electron drift along the <100> and <111> axes.
- `holes::CarrierParameters{T}`: Parameters to describe the hole drift along the <100> and <111> axes.
- `crystal_orientation::SMatrix{3,3,T,9}`: Rotation matrix that transforms the global coordinate system to the crystal coordinate system given by the <100>, <010> and <001> axes of the crystal.
- `γ::SVector{N,SMatrix{3,3,T,9}}`: Reciprocal mass tensors to the `N` valleys of the conduction band.
- `parameters::ADLParameters{T}`: Parameters needed for the calculation of the electron drift velocity.
- `temperaturemodel::TM`: Models to scale the resulting drift velocities with respect to temperature

See also [`CarrierParameters`](@ref).
"""
struct ADLChargeDriftModel{T <: SSDFloat, M <: AbstractDriftMaterial, N, TM <: AbstractTemperatureModel{T}} <: AbstractChargeDriftModel{T}
    electrons::CarrierParameters{T}
    holes::CarrierParameters{T}
    crystal_orientation::SMatrix{3,3,T,9}
    γ::SVector{N,SMatrix{3,3,T,9}}
    parameters::ADLParameters{T}
    temperaturemodel::TM
end

function ADLChargeDriftModel{T,M,N,TM}(chargedriftmodel::ADLChargeDriftModel{<:Any,M,N,TM})::ADLChargeDriftModel{T,M,N,TM} where {T <: SSDFloat, M, N, TM}
    cdmf64 = chargedriftmodel
    e100 = VelocityParameters{T}(cdmf64.electrons.axis100.mu0, cdmf64.electrons.axis100.beta, cdmf64.electrons.axis100.E0, cdmf64.electrons.axis100.mun)
    e111 = VelocityParameters{T}(cdmf64.electrons.axis111.mu0, cdmf64.electrons.axis111.beta, cdmf64.electrons.axis111.E0, cdmf64.electrons.axis111.mun)
    h100 = VelocityParameters{T}(cdmf64.holes.axis100.mu0, cdmf64.holes.axis100.beta, cdmf64.holes.axis100.E0, cdmf64.holes.axis100.mun)
    h111 = VelocityParameters{T}(cdmf64.holes.axis111.mu0, cdmf64.holes.axis111.beta, cdmf64.holes.axis111.E0, cdmf64.holes.axis111.mun)
    electrons = CarrierParameters{T}(e100, e111)
    holes     = CarrierParameters{T}(h100, h111)
    crystal_orientation = SMatrix{3,3,T,9}(cdmf64.crystal_orientation)
    γ = Vector{SMatrix{3,3,T,9}}( cdmf64.γ )
    parameters = ADLParameters{T}(cdmf64.parameters.ml_inv, cdmf64.parameters.mt_inv, cdmf64.parameters.Γ0_inv, cdmf64.parameters.Γ1, cdmf64.parameters.Γ2)
    temperaturemodel::AbstractTemperatureModel{T} = cdmf64.temperaturemodel
    ADLChargeDriftModel{T,M,N,TM}(electrons, holes, crystal_orientation, γ, parameters, temperaturemodel)
end

ADLChargeDriftModel{T}(args...; kwargs...) where {T <: SSDFloat} = ADLChargeDriftModel(args..., T=T, kwargs...)
function ADLChargeDriftModel(configfilename::Union{Missing, AbstractString} = missing; T::Type=Float32, material::Type{<:AbstractDriftMaterial} = HPGe,
                             temperature::Union{Missing, Real}= missing, phi110::Union{Missing, Real} = missing)::ADLChargeDriftModel{T}

    if ismissing(configfilename) configfilename = joinpath(get_path_to_example_config_files(), "ADLChargeDriftModel/drift_velocity_config.yaml") end
    if !ismissing(temperature) temperature = T(temperature) end  #if you give the temperature it will be used, otherwise read from config file

    config = parse_config_file(configfilename)

    #mu0 in m^2 / ( V * s )
    #beta dimensionless
    #E0 in V / m
    #mun in m^2 / ( V * s )

    e100mu0::T  = config["drift"]["velocity"]["parameters"]["e100"]["mu0"]
    e100beta::T = config["drift"]["velocity"]["parameters"]["e100"]["beta"]
    e100E0::T   = config["drift"]["velocity"]["parameters"]["e100"]["E0"]
    e100mun::T  = config["drift"]["velocity"]["parameters"]["e100"]["mun"]
    e111mu0::T  = config["drift"]["velocity"]["parameters"]["e111"]["mu0"]
    e111beta::T = config["drift"]["velocity"]["parameters"]["e111"]["beta"]
    e111E0::T   = config["drift"]["velocity"]["parameters"]["e111"]["E0"]
    e111mun::T  = config["drift"]["velocity"]["parameters"]["e111"]["mun"]
    h100mu0::T  = config["drift"]["velocity"]["parameters"]["h100"]["mu0"]
    h100beta::T = config["drift"]["velocity"]["parameters"]["h100"]["beta"]
    h100E0::T   = config["drift"]["velocity"]["parameters"]["h100"]["E0"]
    h111mu0::T  = config["drift"]["velocity"]["parameters"]["h111"]["mu0"]
    h111beta::T = config["drift"]["velocity"]["parameters"]["h111"]["beta"]
    h111E0::T   = config["drift"]["velocity"]["parameters"]["h111"]["E0"]

    e100 = VelocityParameters{T}(e100mu0, e100beta, e100E0, e100mun)
    e111 = VelocityParameters{T}(e111mu0, e111beta, e111E0, e111mun)
    h100 = VelocityParameters{T}(h100mu0, h100beta, h100E0, 0)
    h111 = VelocityParameters{T}(h111mu0, h111beta, h111E0, 0)
    electrons = CarrierParameters{T}(e100, e111)
    holes     = CarrierParameters{T}(h100, h111)

    if "material" in keys(config)
        config_material::String = config["material"]
        if config_material == "HPGe"
            material = HPGe
        elseif config_material == "Si"
            material = Si
        else
            @warn "Material \"$(config_material)\" not supported for ADLChargeDriftModel.\nSupported materials are \"Si\" and \"HPGe\".\nUsing \"$(Symbol(material))\" as default."
        end
    end

    parameters = if "masses" in keys(config["drift"])
        ml = T(config["drift"]["masses"]["ml"])
        mt = T(config["drift"]["masses"]["mt"])
        ADLParameters{T}(ml, mt, material)
    else
        ml = T(material_properties[Symbol(material)].ml)
        mt = T(material_properties[Symbol(material)].mt)
        ADLParameters{T}(ml, mt, material)
    end

    crystal_orientation::SMatrix{3,3,T,9} = if ismissing(phi110)
        if "phi110" in keys(config)
            RotZ{T}(- π/4 - config["phi110"] )
        else
            transpose(parse_rotation_matrix(T, config, u"rad")) # replace u"rad" with the units in the config file
        end
    else
        RotZ{T}(- π/4 - phi110)
    end

    γ = setup_γj(crystal_orientation, parameters, material)

    if "temperature_dependence" in keys(config)
        if "model" in keys(config["temperature_dependence"])
            model::String = config["temperature_dependence"]["model"]
            if model == "Linear"
                temperaturemodel = LinearModel{T}(config, temperature = temperature)
            elseif model == "PowerLaw"
                temperaturemodel = PowerLawModel{T}(config, temperature = temperature)
            elseif model == "Boltzmann"
                temperaturemodel = BoltzmannModel{T}(config, temperature = temperature)
            else
                temperaturemodel = VacuumModel{T}(config)
                println("Config File does not suit any of the predefined temperature models. The drift velocity will not be rescaled.")
            end
        else
            temperaturemodel = VacuumModel{T}(config)
            println("No temperature model specified. The drift velocity will not be rescaled.")
        end
    else
        temperaturemodel = VacuumModel{T}(config)
        # println("No temperature dependence found in Config File. The drift velocity will not be rescaled.")
    end
    return ADLChargeDriftModel{T,material,length(γ),typeof(temperaturemodel)}(electrons, holes, crystal_orientation, γ, parameters, temperaturemodel)
end

@inline function _get_AE_and_RE(cdm::ADLChargeDriftModel{T, HPGe}, V100e::T, V111e::T)::Tuple{T,T} where{T <: SSDFloat}
    edmp::ADLParameters{T} = cdm.parameters
    AE::T = edmp.Γ0 * V100e
    RE::T = edmp.Γ1 * V111e / AE + edmp.Γ2
    AE, RE
end

@inline function _get_AE_and_RE(cdm::ADLChargeDriftModel{T, Si}, V100e::T, V111e::T)::Tuple{T,T} where{T <: SSDFloat}
    edmp::ADLParameters{T} = cdm.parameters
    AE::T = edmp.Γ0 * V111e
    RE::T = edmp.Γ1 * V100e / AE + edmp.Γ2
    AE, RE
end

@fastmath function getVe(fv::SVector{3, T}, cdm::ADLChargeDriftModel{T,M,N}, Emag_threshold::T = T(1e-5))::SVector{3, T} where {T <: SSDFloat, M, N}
    @inbounds begin
        Emag::T = norm(fv)
        Emag_inv::T = inv(Emag)

        if Emag < Emag_threshold return SVector{3,T}(0, 0, 0) end

        f::NTuple{4,T} = scale_to_given_temperature(cdm.temperaturemodel)
        f100e::T = f[1]
        f111e::T = f[2]
        V100e::T = Vl(Emag, cdm.electrons.axis100) * f100e
        V111e::T = Vl(Emag, cdm.electrons.axis111) * f111e

        AE::T, RE::T = _get_AE_and_RE(cdm, V100e, V111e)

        E0 = SVector{3, T}(fv * Emag_inv)
        oneOverSqrtEγE::SVector{N,T} = broadcast(γ -> T(1/sqrt(γ * E0 ⋅ E0)), cdm.γ)
        sumOneOverSqrtEγE_inv::T = inv(sum(oneOverSqrtEγE))
        N_inv::T = T(1/N)

        g0::SVector{3,T} = @SVector T[0,0,0]
        for j in eachindex(cdm.γ)
            NiOverNj::T = RE * (oneOverSqrtEγE[j] * sumOneOverSqrtEγE_inv - N_inv) + N_inv
            g0 += cdm.γ[j] * E0 * NiOverNj * oneOverSqrtEγE[j]
        end

        return g0 * -AE
    end
end

# Hole model parametrization from [1] equations (22)-(26), adjusted
@fastmath Λ(vrel::T) where {T} = T(0.75 * (1.0 - vrel))

@fastmath function Ω(vrel::T, ::Type{HPGe})::T where {T}
    p1::T = -0.29711
    p2::T = -1.12082
    p3::T = 3.83929
    p4::T = 4.80825
    x::T = 1 - vrel
    p1 * x + p2 * x^2 + p3 * x^3 + p4 * x^4
end

@fastmath function Ω(vrel::T, ::Type{Si})::T where {T}
    p1::T = -0.30565
    p2::T = -1.19650
    p3::T = 4.69001
    p4::T = -7.00635
    x::T = 1 - vrel
    p1 * x + p2 * x^2 + p3 * x^3 + p4 * x^4
end

@fastmath function getVh(fv::SVector{3,T}, cdm::ADLChargeDriftModel{T, M}, Emag_threshold::T = T(1e-5))::SVector{3,T} where {T <: SSDFloat, M}
    @inbounds begin
        Emag::T = norm(fv)
        Emag_inv::T = inv(Emag)

        if Emag < Emag_threshold return SVector{3,T}(0, 0, 0) end

        f::NTuple{4,T} = scale_to_given_temperature(cdm.temperaturemodel)
        f100h::T = f[3]
        f111h::T = f[4]

        V100h::T = Vl(Emag, cdm.holes.axis100) * f100h
        V111h::T = Vl(Emag, cdm.holes.axis111) * f111h

        tmp::SVector{3,T} = cdm.crystal_orientation * fv

        vrel::T = V111h / V100h
        Λvrel::T = Λ(vrel)
        Ωvrel::T = Ω(vrel, M)
        θ0::T = acos(tmp[3] / Emag)
        φ0::T = atan(tmp[2], tmp[1])
        sθ0::T, cθ0::T = sincos(θ0)
        s2θ0::T, c2θ0::T = sincos(2*θ0)
        s2φ0::T, c2φ0::T = sincos(2*φ0)

        vr::T = V100h * ( 1 - Λvrel * (sθ0^4 * s2φ0^2 + s2θ0^2) )
        vΩ::T = V100h * Ωvrel * 2 * (sθ0^3 * cθ0 * s2φ0^2 + s2θ0 * c2θ0)
        vφ::T = V100h * Ωvrel * 2 * sθ0^3 * s2φ0 * c2φ0

        Ry::RotY{T} = RotY{T}(θ0)
        Rz::SMatrix{3,3,T,9} = transpose(cdm.crystal_orientation) * RotZ{T}(φ0)
        Rz * (Ry * @SVector T[vΩ, vφ, vr])
    end
end


text:
struct Wrapper{NT}
   x::NT
end

@generated function syms(w::Wrapper{NamedTuple{SYMS, TT}}
                         ) where {SYMS, TT}
   code = "(" * prod( (":" * string(sym) * ",") for sym in SYMS ) * ")"
   quote
      $(Meta.parse(code))
   end
end

nt = (a = 1, b = 2, c = 5)
w = Wrapper(nt)

syms(w)


text:
@testset "Special Cases" begin

phase = ones(3,3,3)
unwrap(phase)

point = ones(1,1,1)
@test_throws AssertionError unwrap(point)

## weight test
function weight_test(w, out)
    w = reshape(w, length(w), 1, 1)
    @test UInt8.(out) == ROMEO.calculateweights(w)[1,:,1,1]
end
weight_test([0.1, 0.2 + 2pi, 0.3, 0.4], [30, 7, 30, 0]) # phase linearity penalty (30) at borders
weight_test([0.1, 0.2 + 2pi, 0.3, NaN], [30, 157, 0, 0]) # 157 to voxel bordering to NaN voxel, 0 to NaN voxel

## NaN test
@test nan_test(unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]), [0.1, 0.2, 0.3, 0.4])
@test nan_test(unwrap([0.1, 0.2 + 2pi, 0.3, NaN]), [0.1, 0.2, 0.3, NaN])
@test nan_test(unwrap([0.1, 0.2 + 2pi, 0.3, 0.4]; mag=[1, 1, 1, NaN]), [0.1, 0.2, 0.3, 0.4])

## test dimensions
function dim_test(a; keyargs...)
    ca = copy(a)
    ua = unwrap(a; keyargs...)
    @test a != ua
    @test size(a) == size(ua)
    @test a == ca
    unwrap!(a; keyargs...)
    @test a != ca
    @test a == ua
end

dim_test(10rand(50))
dim_test(10rand(50,20))
dim_test(10rand(15,20,10))
dim_test(10rand(7,9,10,3); TEs=ones(3))

end


text:
module MDSim

export init, simulate!

using Statistics: mean

mutable struct MDSystem{DT<:AbstractFloat}
    l::DT
    h::DT
    U::DT
    K::DT
    T::DT
    P::DT
    N::Integer
    f::Matrix{DT}
    r::Matrix{DT}
    v::Matrix{DT}
    function MDSystem(N::Integer, T₀::DT, h::DT, l::DT) where {DT<:AbstractFloat}
        r = rand(DT, 2, N) * l
        v = rand(DT, 2, N) .- 0.5
        v *= √(T₀ ÷ mean(v .^ 2))
        f = Matrix{DT}(undef, 2, N)
        return new{DT}(l, h, DT(NaN), DT(NaN), DT(NaN), DT(NaN), N, f, r, v)
    end
end

function init(; N::Integer, T₀::DT, h::DT, l::DT) where {DT<:AbstractFloat}
    sys = MDSystem(N, T₀, h, l)
    alignleft!(sys)
    update_forces!(sys)
    update_phase!(sys)
    check_boundary!(sys)
    update_potential!(sys)
    update_kinetic!(sys)
    update_temperature!(sys)
    update_pressure!(sys)
    return sys
end

function simulate!(sys::MDSystem)
    update_phase!(sys)
    check_boundary!(sys)
    update_potential!(sys)
    update_kinetic!(sys)
    update_temperature!(sys)
    update_pressure!(sys)
end

function alignleft!(sys::MDSystem)    #Thanks to my dearest friend for his beautiful method in setting the initial conditions SLH Hero
    xcount = ceil(Integer, (2 * sys.N)^(1 / 2) / 2)
    ycount = 2 * xcount

    xs = collect(range(0.1 * sys.l / 2, 0.9 * sys.l / 2, length = xcount))
    sys.r[1, :] .= repeat(xs, outer = ycount)[1:sys.N]

    ys = collect(range(0.1 * sys.l, 0.9 * sys.l, length = ycount))
    sys.r[2, :] .= repeat(ys, inner = xcount, outer = 1)[1:sys.N]
end

function update_forces!(sys::MDSystem)
    for i ∈ 1:sys.N
        sys.f[:, i] *= 0
        for j ∈ setdiff(1:sys.N, i)
            rmirror = sys.r[:, j]
            if sys.r[1, i] - rmirror[1] > sys.l / 2
                rmirror[1] += sys.l
            elseif sys.r[1, i] - rmirror[1] < -sys.l / 2
                rmirror[1] -= sys.l
            end
            if sys.r[2, i] - rmirror[2] > sys.l / 2
                rmirror[2] += sys.l
            elseif sys.r[2, i] - rmirror[2] < -sys.l / 2
                rmirror[2] -= sys.l
            end
            Δr = √sum((sys.r[:, i] - rmirror) .^ 2)
            sys.f[:, i] += 48 * (((Δr)^-14) - 0.5 * ((Δr)^-8)) * (sys.r[:, i] - rmirror)
        end
    end
end

function update_force!(sys::MDSystem, i::Integer)
    sys.f[:, i] *= 0
    for j ∈ setdiff(1:sys.N, i)
        rmirror = sys.r[:, j]
        if sys.r[1, i] - rmirror[1] > sys.l / 2
            rmirror[1] += sys.l
        elseif sys.r[1, i] - rmirror[1] < -sys.l / 2
            rmirror[1] -= sys.l
        end
        if sys.r[2, i] - rmirror[2] > sys.l / 2
            rmirror[2] += sys.l
        elseif sys.r[2, i] - rmirror[2] < -sys.l / 2
            rmirror[2] -= sys.l
        end
        Δr = √sum((sys.r[:, i] - rmirror) .^ 2)
        sys.f[:, i] += 48 * (((Δr)^-14) - 0.5 * ((Δr)^-8)) * (sys.r[:, i] - rmirror)
    end
end

"""
important note:
    update_forces!() will update the force of the whole system,
    but update_force! only updates the selected particle.
"""

function update_phase!(sys::MDSystem)
    for i ∈ 1:sys.N
        f₁ = sys.f[:, i]
        sys.r[:, i] += sys.h * sys.v[:, i] + 0.5 * sys.h * sys.h * f₁
        update_force!(sys, i)
        f₂ = sys.f[:, i]
        sys.v[:, i] += sys.h * (f₁ + f₂) / 2
    end
end

function check_boundary!(sys::MDSystem)
    sys.r = (sys.r .+ sys.l) .% sys.l
end

function update_potential!(sys::MDSystem)
    Σrterm = 0.0

    for i ∈ 1:(sys.N-1)
        for j ∈ i+1:sys.N
            rmirror = sys.r[:, j]
            if sys.r[1, i] - rmirror[1] > sys.l / 2
                rmirror[1] += sys.l
            elseif sys.r[1, i] - rmirror[1] < -sys.l / 2
                rmirror[1] -= sys.l
            end
            if sys.r[2, i] - rmirror[2] > sys.l / 2
                rmirror[2] += sys.l
            elseif sys.r[2, i] - rmirror[2] < -sys.l / 2
                rmirror[2] -= sys.l
            end
            Δr = √sum((sys.r[:, i] - rmirror) .^ 2)
            Σrterm += ((Δr)^-12) - ((Δr)^-6)
        end
    end

    U = 4 * Σrterm

    sys.U = U
end

function update_kinetic!(sys::MDSystem)
    sys.K = sum(sys.v .^ 2) / 2
end

function update_temperature!(sys::MDSystem)
    sys.T = sum(sys.v .^ 2) / (2 * sys.N)
end

function update_pressure!(sys::MDSystem)
    Σv² = sum(sys.v .^ 2)
    Σrterm = 0.0

    for i ∈ 1:(sys.N-1)
        for j ∈ i+1:sys.N
            rmirror = sys.r[:, j]
            if sys.r[1, i] - rmirror[1] > sys.l / 2
                rmirror[1] += sys.l
            elseif sys.r[1, i] - rmirror[1] < -sys.l / 2
                rmirror[1] -= sys.l
            end
            if sys.r[2, i] - rmirror[2] > sys.l / 2
                rmirror[2] += sys.l
            elseif sys.r[2, i] - rmirror[2] < -sys.l / 2
                rmirror[2] -= sys.l
            end
            Δr = √sum((sys.r[:, i] - rmirror) .^ 2)
            Σrterm += ((Δr)^-12) - 0.5 * ((Δr)^-6)
        end
    end
    sys.P = ((48 * Σrterm) + Σv²) / (2 * sys.l * sys.l)
end

end


text:
#ДАНО: Робот - Робот - в произвольной клетке ограниченного прямоугольного поля
#РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так: нижний ряд - полностью, следующий - весь, за исключением одной последней клетки на Востоке, следующий - за исключением двух последних клеток на Востоке, и т.д.
function begin!(r)
    for side in (West,Sud)
        while !isborder(r,side)
            move!(r,side)
        end
    end
end

function step(r)
    num = 0
    begin!(r)
    while !isborder(r,Ost)
        move!(r,Ost)
        num+=1
    end
    return num
    begin!(r)
end

function craska!(r)
    a = step(r)
    begin!(r)
    while a!=0
        for _ in 1:a
            putmarker!(r)
            move!(r,Ost)
        end
        putmarker!(r)
        a = a - 2
        move!(r,Nord)
        while !isborder(r,West)
            move!(r,West)
        end
    end
end


text:
#
#  This example is based on demo_m1.m from LAYPACK 1.0 and includes
#  many of the same comments.
#
#  MODEL REDUCTION BY THE ALGORITHMS LRSRM AND DSPMR. THE GOAL IS TO
#  GENERATE A REDUCED SYSTEM OF VERY SMALL ORDER.
#
#  This demo program shows the use of 'lp_lradi' and 'lp_para'

# -----------------------------------------------------------------------
# Generate test problem
# -----------------------------------------------------------------------
#
# This is an artificial test problem of a system, whose Bode plot shows
# "spires".
workspace()
include("../src/LargeMatrixEquations.jl")
using LargeMatrixEquations

A = speye(408,408); B = ones(408,1); C = ones(408,1)
A[1:2,1:2] = [-.01 -200; 200 .001]
A[3:4,3:4] = [-.2 -300; 300 -.1]
A[5:6,5:6] = [-.02 -500; 500 0]
A[7:8,7:8] = [-.01 -520; 520 -.01]
A[9:408,9:408] = diagm(-(1:400))

println("Problem dimensions:")

n = size(A,1);   # problem order (number of states)
m = size(B,2);  # number of inputs
q = size(C,1);   # number of outputs

# -----------------------------------------------------------------------
# Initialization/generation of data structures used in user-supplied
# functions and computation of ADI shift parameters
# -----------------------------------------------------------------------
#
# Note that A is a tridiagonal matrix. No preprocessing needs to be done.

println("Parameters for heuristic algorithm which computes ADI parameters:")
l0 = 10;   # desired number of distinct shift parameters
kp = 30;   # number of steps of Arnoldi process w.r.t. A
km = 15;   # number of steps of Arnoldi process w.r.t. inv(A)

b0 = ones(n,1);   # This is just one way to choose the Arnoldi start
                  # vector.
Bf = []#zeros(408,2)#
Kf = []
# computation of ADI shift parameters
p = lp_para(A,Bf,Kf,l0,kp,km,b0)

println("Actual number of ADI shift parameters:")
l = length(p)

println("ADI shift parameters:")
p

# -----------------------------------------------------------------------
# Solution of Lyapunov equations A*X+X*A' = -B*B' and
# A'*X+X*A = -C'*C
# -----------------------------------------------------------------------

println("Parameters for stopping criteria in LRCF-ADI iteration:")
max_it = 40;   # (will stop the iteration)
min_res = 1e-100;   # (avoided, but the residual history is shown)
with_rs = :N;   # (avoided)
min_in = 0;   # (avoided)

zk = :Z
rc = :R
Bf = ones(408,2)#[]
Kf = zeros(408,2)#[]
K = ones(408,408)
info = 2

println("... solving A*XB+XB*A'' = - B*B''...")
tp = :B

ZB,flag,res =lp_lradi(A,B,p=p,max_it=200)

X=ZB*ZB'
println(norm(A*X+X*A'+B*B'))

max_stars_repo_name: garrettthomaskth/LargeMatrixEquations

max_issues_repo_name: garrettthomaskth/LargeMatrixEquations

max_forks_repo_name: garrettthomaskth/LargeMatrixEquations


text:
using EasyConfig
using StructTypes
using OrderedCollections: OrderedDict
using Test

@testset "Constructors" begin
    @test Config() == Config()
    c = Config()
    c.x = 1
    @test Config(x=1) == Config(:x => 1) == Config(Dict(:x => 1)) == Config((; x=1)) == c
end

@testset "Set/get property/field" begin
    c = Config()
    c.one."two"["three"][:four] = 5
    @test c.one.two.three.four == 5
    c."test" == Config()
    @test propertynames(c) == [:one, :test]
end

@testset "Other base methods" begin
    c = Config()

    @test get(c, "test", 1) == 1
    @test c."test" == Config()

    @test get!(c, "test2", 1) == 1
    @test c."test2" == 1
    delete!(c, "test")
    @test propertynames(c) == [:test2]

    @test keys(Config(x = 1)) == keys(OrderedDict(:x => 1))
    @test all(values(Config(x = 1)) .== values(OrderedDict(:x => 1)))
    @test pairs(Config(x = 1)) == pairs(OrderedDict(:x => 1))
    @test empty!(Config(x=1)) == Config()
end

@testset "from NamedTuple" begin
    nt = (x=1, y=(x=1, y=(x=1, y=(x=1,y=2))))
    c = Config(nt)
    @test c.y.y.y.y == 2
end

@testset "from Dict" begin
    d = Dict(:x => Dict(:x => Dict(:x => Dict(:x => 2))))
    c = Config(d)
    @test c.x.x.x.x == 2
end

@testset "merge!" begin
    a = Config(x = 1, y = Config(x = 1))
    b = Config(x = 5, y = Config(x=5, z="hi"))
    c = merge(a, b)
    merge!(a, b)
    @test a == c
    @test a.x == 5
    @test a.y.x == 5
    @test a.y.z == "hi"
end

@testset "isempty/delete_empty!" begin
    c = Config()
    c.x.x.x
    @test isempty(c) == true
end

@testset "StructTypes" begin
    @test isempty(StructTypes.keyvaluepairs(Config()))
    c = Config(x=1,y=2)
    @test StructTypes.keyvaluepairs(c)[:x] == 1
    @test StructTypes.keyvaluepairs(c)[:y] == 2
end


text:
#Clearing Error Function
function  z(x::Matrix{Int64}, q::Vector{Int64}, p::Vector{Int64})
   # S is the number of students, N is the number of courses
   # x is the SxN allocation matrix
   # q is the Nx1 course capacity vector
   # p is the Nx1 price vector
   N = length(q)
   zer = zeros(N)        #Create empty Nx1 vector of clearing errors
   numstud = transpose(sum(x, 1)) #Create Nx1 vector with the number of students enrolled in each course.
   for i in 1:N        #For every course
       if p[i] > 0     #As long as the price in this course is positive.
           zer[i] = numstud[i] - q[i]
       else
           zer[i] = max(numstud[i] - q[i], 0)
       end
   end
   return zer
end

#TEST
x = [1 0 1;
     1 0 1;
     1 0 1;
     0 0 1;
     1 1 0;
     1 1 1]
q = [4;3;3]
N = 3
p = [30; 0; 50]

z(x, q, p, N)

max_stars_repo_name: ScPo-CompEcon/CourseMatch

max_issues_repo_name: ScPo-CompEcon/CourseMatch

max_forks_repo_name: ScPo-CompEcon/CourseMatch


text:
# Autogenerated wrapper script for gperf_jll for x86_64-linux-musl
export gperf

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `gperf`
const gperf_splitpath = ["bin", "gperf"]

# This will be filled out by __init__() for all products, as it must be done at runtime
gperf_path = ""

# gperf-specific global declaration
function gperf(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(gperf_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"gperf")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    global gperf_path = normpath(joinpath(artifact_dir, gperf_splitpath...))

    push!(PATH_list, dirname(gperf_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')

end  # __init__()