text:
"""
    YASStyle()

Formatting style based on [YASGuide](https://github.com/jrevels/YASGuide)
and [JuliaFormatter#198](https://github.com/domluna/JuliaFormatter.jl/issues/198).

Configurable options with different defaults to [`DefaultStyle`](@ref) are:
- `always_for_in` = true
- `whitespace_ops_in_indices` = true
- `remove_extra_newlines` = true
- `import_to_using` = true
- `pipe_to_function_call` = true
- `short_to_long_function_def` = true
- `always_use_return` = true
- `whitespace_in_kwargs` = false
"""
struct YASStyle <: AbstractStyle
    innerstyle::Union{Nothing,AbstractStyle}
end
YASStyle() = YASStyle(nothing)
@inline getstyle(s::YASStyle) = s.innerstyle === nothing ? s : s.innerstyle

function options(style::YASStyle)
    return (;
        always_for_in = true,
        whitespace_ops_in_indices = true,
        remove_extra_newlines = true,
        import_to_using = true,
        pipe_to_function_call = true,
        short_to_long_function_def = true,
        always_use_return = true,
        whitespace_in_kwargs = false,
    )
end

function is_binaryop_nestable(::YASStyle, cst::CSTParser.EXPR)
    (CSTParser.defines_function(cst) || is_assignment(cst)) && return false
    (cst[2].val == "=>" || cst[2].val == "->") && return false
    return true
end

function p_import(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = p_import(DefaultStyle(style), cst, s)
    idx = findfirst(n -> n.typ === PLACEHOLDER, t.nodes)
    if idx !== nothing && is_colon(t[idx-1])
        t[idx] = Whitespace(1)
    end
    t
end

function p_using(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_import(ys, cst, s)
    t.typ = Using
    t
end

function p_export(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_import(ys, cst, s)
    t.typ = Export
    t
end

function p_curly(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Curly, cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if CSTParser.is_comma(a) && i == length(cst) - 1
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(0), s)
        elseif is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif i > 1 && is_opener(cst[i-1])
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

function p_braces(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_curly(ys, cst, s)
    t.typ = Braces
    t
end

function p_tuple(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(TupleN, cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = if is_binary(a) && a[2].val == "="
            p_kw(style, a, s)
        else
            pretty(style, a, s)
        end

        if CSTParser.is_comma(a) && i + 1 == length(cst)
            if length(cst.args) == 1
                add_node!(t, n, s, join_lines = true)
            elseif !is_closer(cst[i+1])
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Placeholder(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif i > 1 && is_opener(cst[i-1])
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

function p_tuple(ys::YASStyle, nodes::Vector{CSTParser.EXPR}, s::State)
    style = getstyle(ys)
    t = FST(TupleN, nspaces(s))
    for (i, a) in enumerate(nodes)
        n = pretty(style, a, s)
        if CSTParser.is_comma(a) && i + 1 == length(nodes)
            if length(nodes) == 1
                add_node!(t, n, s, join_lines = true)
            elseif !is_closer(nodes[i+1])
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Placeholder(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(nodes) && !is_punc(nodes[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif i > 1 && is_opener(nodes[i-1])
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

# Brackets
function p_invisbrackets(
    ys::YASStyle,
    cst::CSTParser.EXPR,
    s::State;
    nonest = false,
    nospace = false,
)
    style = getstyle(ys)
    t = FST(Brackets, cst, nspaces(s))

    if is_block(cst[2]) || (cst[2].head === :generator && is_block(cst[2][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(
        t,
        pretty(style, cst[1], s),
        s,
        join_lines = true,
        override_join_lines_based_on_source = true,
    )

    n = if cst[2].head === :block
        pretty(style, cst[2], s, from_quote = true)
    elseif is_opcall(cst[2])
        pretty(style, cst[2], s, nonest = nonest, nospace = nospace)
    else
        pretty(style, cst[2], s)
    end
    add_node!(t, n, s, join_lines = true, override_join_lines_based_on_source = true)

    add_node!(
        t,
        pretty(style, cst[3], s),
        s,
        join_lines = true,
        override_join_lines_based_on_source = true,
    )

    t
end

function p_call(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Call, cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if CSTParser.is_comma(a) && i + 1 == length(cst)
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif i > 1 && is_opener(cst[i-1])
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    if !parent_is(cst, is_function_or_macro_def)
        separate_kwargs_with_semicolon!(t)
    end
    t
end
function p_vect(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_call(ys, cst, s)
    t.typ = Vect
    t
end

function p_vcat(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Vcat, nspaces(s))
    st = cst.head === :vcat ? 1 : 2

    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif !is_closer(a) && i > st
            join_lines = i == st + 1 ? true : t.endline == n.startline
            if join_lines && i != st + 1
                add_node!(t, Placeholder(1), s)
            end

            add_node!(
                t,
                n,
                s;
                join_lines = join_lines,
                override_join_lines_based_on_source = i == st + 1,
            )
            if i != length(cst) - 1
                has_semicolon(s.doc, n.startline) &&
                    add_node!(t, InverseTrailingSemicolon(), s)
            end
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end
function p_typedvcat(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_vcat(ys, cst, s)
    t.typ = TypedVcat
    t
end

function p_ref(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(RefN, cst, nspaces(s))
    nospace = !s.opts.whitespace_ops_in_indices

    for (i, a) in enumerate(cst)
        n = if is_binary(a) || is_chain(a) || a.head === :brackets || a.head === :comparison
            pretty(style, a, s, nonest = true, nospace = nospace)
        else
            pretty(style, a, s)
        end

        if CSTParser.is_comma(a) && i + 1 == length(cst)
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif i > 1 && is_opener(cst[i-1])
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

function p_comprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Comprehension, cst, nspaces(s))

    if is_block(cst[2]) || (cst[2].head === :Generator && is_block(cst[2][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(t, pretty(style, cst[1], s), s, join_lines = true)
    add_node!(
        t,
        pretty(style, cst[2], s),
        s,
        join_lines = true,
        override_join_lines_based_on_source = true,
    )
    add_node!(
        t,
        pretty(style, cst[3], s),
        s,
        join_lines = true,
        override_join_lines_based_on_source = true,
    )
    t
end

function p_typedcomprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(TypedComprehension, cst, nspaces(s))

    if is_block(cst[3]) || (cst[3].head === :generator && is_block(cst[3][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(t, pretty(style, cst[1], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[2], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[3], s), s, join_lines = true)
    add_node!(t, pretty(style, cst[4], s), s, join_lines = true)
    t
end

function p_macrocall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(MacroCall, cst, nspaces(s))

    args = get_args(cst)
    nest = length(args) > 0 && !(length(args) == 1 && unnestable_node(args[1]))
    has_closer = is_closer(cst[end])

    !has_closer && (t.typ = MacroBlock)

    # same as CSTParser.Call but whitespace sensitive
    for (i, a) in enumerate(cst)
        if CSTParser.is_nothing(a)
            s.offset += a.fullspan
            continue
        end

        n = pretty(style, a, s)
        if CSTParser.ismacroname(a)
            if has_closer || length(args) == 0
                add_node!(t, n, s, join_lines = true)
            else
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Whitespace(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif is_closer(n)
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif i > 1 && is_opener(cst[i-1])
            add_node!(
                t,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif t.typ === MacroBlock
            if has_closer
                add_node!(t, n, s, join_lines = true)
                if i < length(cst) - 1 && cst[i+1].head != :parameters
                    add_node!(t, Whitespace(1), s)
                end
            else
                padding = is_block(n) ? 0 : -1
                add_node!(t, n, s, join_lines = true, max_padding = padding)
                i < length(cst) && add_node!(t, Whitespace(1), s)
            end
        else
            if has_closer
                add_node!(t, n, s, join_lines = true)
            else
                padding = is_block(n) ? 0 : -1
                add_node!(t, n, s, join_lines = true, max_padding = padding)
            end
        end
    end
    # move placement of @ to the end
    #
    # @Module.macro -> Module.@macro
    t[1] = move_at_sign_to_the_end(t[1], s)
    t
end

function p_whereopcall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Where, cst, nspaces(s))

    add_node!(t, pretty(style, cst[1], s), s)
    add_node!(t, Whitespace(1), s)
    add_node!(t, pretty(style, cst[2], s), s, join_lines = true)
    add_node!(t, Whitespace(1), s)

    curly_ctx =
        cst.parent.head === :curly ||
        cst[3].head === :curly ||
        cst[3].head === :bracescat ||
        cst[3].head === :parameters

    add_braces = !curly_ctx && !CSTParser.is_lbrace(cst[3])

    bc = curly_ctx ? t : FST(BracesCat, nspaces(s))

    if add_braces
        brace = FST(PUNCTUATION, -1, t.endline, t.endline, "{")
        add_node!(bc, brace, s, join_lines = true)
    end

    for i = 3:length(cst)
        a = cst[i]
        n = is_binary(a) ? pretty(style, a, s, nospace = true) : pretty(style, a, s)

        if CSTParser.is_comma(a) && i == length(cst) - 1
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(bc, n, s, join_lines = true)
            add_node!(bc, Placeholder(0), s)
        elseif is_closer(n)
            add_node!(
                bc,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        elseif is_opener(cst[i-1])
            add_node!(
                bc,
                n,
                s,
                join_lines = true,
                override_join_lines_based_on_source = true,
            )
        else
            add_node!(bc, n, s, join_lines = true)
        end
    end

    if add_braces
        brace = FST(PUNCTUATION, -1, t.endline, t.endline, "}")
        add_node!(bc, brace, s, join_lines = true)
    end

    !curly_ctx && add_node!(t, bc, s, join_lines = true)

    t
end

function p_generator(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    style = getstyle(ys)
    t = FST(Generator, cst, nspaces(s))
    has_for_kw = false

    for (i, a) in enumerate(cst)
        n = pretty(style, a, s)
        if CSTParser.iskeyword(a)
            if !has_for_kw && a.head === :FOR
                has_for_kw = true
            end

            incomp =
                parent_is(a, is_iterable, ignore = n -> is_gen(n) || n.head === :brackets)

            if is_block(cst[i-1])
                add_node!(t, Newline(), s)
            elseif incomp
                add_node!(t, Placeholder(1), s)
            else
                add_node!(t, Whitespace(1), s)
            end
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Whitespace(1), s)

            if !is_gen(cst[i+1])
                tupargs = CSTParser.EXPR[]
                for j = i+1:length(cst)
                    push!(tupargs, cst[j])
                end

                # verify this is not another for loop
                any(b -> b.head === :FOR, tupargs) && continue

                tup = p_tuple(style, tupargs, s)
                add_node!(t, tup, s, join_lines = true)

                if has_for_kw
                    for nn in tup.nodes
                        eq_to_in_normalization!(nn, s.opts.always_for_in)
                    end
                end
                break
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end

        has_for_kw && eq_to_in_normalization!(n, s.opts.always_for_in)
    end

    t
end

function p_filter(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_generator(ys, cst, s)
    t.typ = Filter
    t
end

function p_flatten(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = p_generator(ys, cst, s)
    t.typ = Flatten
    t
end


text:
using ClickHouse: ClickHouseSock, CHSettings, is_connected,
                is_busy, chwrite, chread, has_temporary_tables, ClientInfo,
                VarUInt, write_packet, read_packet, @using_socket, ClientHello

@testset "guarded" begin
    sock = ClickHouseSock(PipeBuffer())
    ClickHouse.@guarded sock sock.busy = true
    @test sock.busy == true
    @test ClickHouse.@guarded(sock, sock.busy) == true
    name = "test test"
    ClickHouse.@guarded sock begin
        sock.server_name = name
    end
    @test sock.server_name == "test test"
end
@testset "busy" begin
    sock = ClickHouseSock(nothing)

    try
        @using_socket sock begin
            sleep(1)
        end
        @test false
    catch e
        @test e.msg == "ClickHouseSock not connected"
    end

    sock.io = PipeBuffer()

    a = @async @using_socket sock begin
        sleep(1)
    end
    sleep(0.2)
    @test is_busy(sock)
    try
        @using_socket sock begin
            sleep(1)
        end
        @test false
    catch e
        @test e.msg == "ClickHouseSock is busy"
    end
    wait(a)
    @test !is_busy(sock)
end

@testset "ch structs" begin
    c_info = ClientInfo(
        0,
        "user",
        "aaa-aaa",
        ":0",
        1,
        "osu",
        "host",
        "name",
        10,
        2,
        23331,
        "quota",
        10
    )

    sock = ClickHouseSock(PipeBuffer())
    chwrite(sock, c_info)
    res = chread(sock, ClientInfo)

    @test res == ClientInfo(
        0,
        "user",
        "aaa-aaa",
        ":0",
        1,
        "osu",
        "host",
        "name",
        10,
        2,
        23331,
        "",
        0
    )

    sock.server_rev = ClickHouse.has_quota_key_rev()
    chwrite(sock, c_info)
    res = chread(sock, ClientInfo)

    @test res == ClientInfo(
        0,
        "user",
        "aaa-aaa",
        ":0",
        1,
        "osu",
        "host",
        "name",
        10,
        2,
        23331,
        "quota",
        0
    )

    sock.server_rev = ClickHouse.has_version_patch_rev()
    chwrite(sock, c_info)
    res = chread(sock, ClientInfo)

    @test res == ClientInfo(
        0,
        "user",
        "aaa-aaa",
        ":0",
        1,
        "osu",
        "host",
        "name",
        10,
        2,
        23331,
        "quota",
        10
    )
end
@testset "ch packets" begin
    c_hello = ClientHello(
        "ddddd",
        10,
        10,
        10,
        "db",
        "us",
        "pass"
    )

    sock = ClickHouseSock(PipeBuffer())
    write_packet(sock, c_hello)
    res = read_packet(sock, ClickHouse.ClientCodes)
    @test res isa ClientHello

    sock = ClickHouseSock(PipeBuffer())
    write_packet(sock, ClickHouse.ServerPong())
    res = read_packet(sock, ClickHouse.ServerCodes)
    @test res isa ClickHouse.ServerPong
end

ext: jl


text:
using Test
using jInv.Mesh

h = rand(6)
msh = getTensorMesh3D(h,h,h)
m = rand(msh.nc)

Ne, Qe, activeEdges = getEdgeConstraints(msh);
Nf,Qf = getFaceConstraints(msh);

Curl = getCurlMatrix(msh)
Msig = getEdgeMassMatrix(msh, m)
Mmu = getFaceMassMatrix(msh, m)

CurlC = Qf  * Curl * Ne
MsigC = Ne' * Msig * Ne
MmuC  = Nf' * Mmu * Nf

@test Curl==CurlC
@test Msig==MsigC
@test Mmu==MmuC


text:
"""
	gmtconnect(cmd0::String="", arg1=[], kwargs...)

Connect individual lines whose end points match within tolerance

Full option list at [`gmtconnect`](http://gmt.soest.hawaii.edu/doc/latest/gmtconnect.html)

Parameters
----------

- **C** : **closed** : -- Str or [] --        Flags = [closed]

    Write all the closed polygons to closed [gmtgmtconnect_closed.txt] and return all other
    segments as they are. No gmtconnection takes place.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/gmtconnect.html#c)
- **D** : **dump** : -- Str or [] --   Flags = [template]

    For multiple segment data, dump each segment to a separate output file
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/gmtconnect.html#d)
- **L** : **linkfile** : -- Str or [] --      Flags = [linkfile]

    Writes the link information to the specified file [gmtgmtconnect_link.txt].
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/gmtconnect.html#l)
- **Q** : **list_file** : -- Str or [] --      Flags =  [listfile]

    Used with **D** to write a list file with the names of the individual output files.
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/gmtconnect.html#q)
- **T** : **tolerance ** : -- List or Str --    Flags = [cutoff[unit][/nn_dist]]

    Specifies the separation tolerance in the data coordinate units [0]; append distance unit.
    If two lines has end-points that are closer than this cutoff they will be joined.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmtconnect.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtconnect(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("gmtconnect", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'C', d, [:C :closed])
	cmd = add_opt(cmd, 'D', d, [:D :dump])
	cmd = add_opt(cmd, 'L', d, [:L :linkfile])
	cmd = add_opt(cmd, 'Q', d, [:Q :list_file])
	cmd = add_opt(cmd, 'T', d, [:T :tolerance ])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "gmtconnect", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtconnect(arg1=[], cmd0::String=""; kw...) = gmtconnect(cmd0, arg1; kw...)


text:
using Distributions: truncated

## Extra distributions ##

@dist shifted_poisson(rate, offset::Int) = poisson(rate) + offset
@dist delta(val) = uniform_discrete(val, val)

struct TruncatedNormal <: Gen.Distribution{Float64} end

"""
    trunc_normal(mu::Real, std::Real, lb::Real, ub::Real)
Samples a `Float64` value from a normal distribution.
"""
const trunc_normal = TruncatedNormal()

(d::TruncatedNormal)(mu, std, lb, ub) = Gen.random(d, mu, std, lb, ub)

Gen.random(::TruncatedNormal, mu::Real, std::Real, lb::Real, ub::Real) =
    rand(truncated(Distributions.Normal(mu, std), lb, ub))

function Gen.logpdf(::TruncatedNormal, x::Real, mu::Real, std::Real, lb::Real, ub::Real)
    d = truncated(Distributions.Normal(mu, std), lb, ub)
    untrunc_lpdf = Distributions.logpdf(d.untruncated, x)
    if d.tp > 0
        untrunc_lpdf - d.logtp
    elseif cdf(d.untruncated, lb) ≈ 1.0
        untrunc_lpdf - Distributions.logccdf(d.untruncated, lb)
    elseif cdf(d.untruncated, ub) ≈ 0.0
        untrunc_lpdf - Distributions.logcdf(d.untruncated, ub)
    end
end

Gen.is_discrete(::TruncatedNormal) = false

## Other utility functions ##

function segment_find(idx, seglens)
    count = 0
    for (i, l) in enumerate(seglens)
        count += l
        if count >= idx return (i, idx-count+l) end
    end
    return nothing
end


text:
module SupportFunc
export updaterec2, getrisk2, pidcontrol;

"""

    pidcontrol(crthresh, pK, dK, iK, conchist, Dt, sig, exp_adv, expnum_advper)

Proportional-Derivative-Integral (PID) control of concentration or risk using a gaussian weight.

Arguments:

- `crthresh`: concentration or risk threshold used as control reference.
- `pK`: proportional coefficient.
- `dK`: derivative coefficient.
- `iK`: integral coefficient.
- `conchist`: contaminant concentration record (all years until then)
- `Dt`: integration time increments.
- `sig`: gaussian weight standard deviation.
- `exp_adv`: expert influence (0), or media influence (1).
- `expnum_advper` : number of experts (if `exp_adv = 0`), or advertising intervals (if `exp_adv = 1`).

Returns:

- `belief` : boolean belief of experts or advertisers.
- `exp_ad` : number of experts (if `exp_adv = 0`), or advertising intervals (if `exp_adv = 1`).

"""

function pidcontrol(crthresh, pK, dK, iK, conchist, Dt, sig, exp_adv, expnum_advper)

	e = conchist[end] - crthresh;
	if(pK+dK+iK != 0)
		if(length(conchist) > 2)
			deriv = ( (conchist[end] - conchist[end-1])/Dt + (conchist[end-1] - conchist[end-2])/Dt )/2;
			integ = sum(conchist .- crthresh) * Dt;
		else
			deriv = 0;
			integ = 0;
		end

		K = pK * e - dK * deriv + iK * integ;

		if(exp_adv == 0)
			exp_ad = round(Int, expnum_advper*(1 - exp(-( K^2 / (2*sig^2)))));
		else
			exp_ad = round(Int, 1/(expnum_advper*(1 - exp(-( K^2 / (2*sig^2))))));
		end

	else
		exp_ad = expnum_advper;
	end

	if(e > 0)
		belief = true;
	else
		belief = false;
	end

	return belief, exp_ad

end

###########################################################

# This function does not wait for people to die in order to assign risk for their usage
# Instead it assumes that people will use water their whole life based on the usage of the last year.

"""

    updaterec2(soc, cc, period, agelimit)

Updates the `soc` array containing social and contamination record.
This method of storing records assumes that individuals will be exposed with the same concentration as that of this year for the rest of their life (for calculating `cont-accumulation`).

Arguments:

- `soc`: `[belief, usage, cont-accumulation, age]` of every individual of the past year
- `cc`: contaminant concentration
- `period` : coupling time increment (years)
- `agelimit` : age limit of each individual

Returns:

- `soc`

"""
function updaterec2(soc, cc, period, agelimit)

soc[:,4] = soc[:,4] .+ period; 				# Adding years to their age.
soc[:,3] = soc[:,3] + !map(Bool, soc[:,2]) .* cc;		# Assuming their usage is the same for all following years.
background_concentration = 5;

# for all those who die, a new person is born.
for i = 1:length(soc[:,1])
	if(soc[i,4] >= agelimit)

		# For each person who dies, one is born, clean, age 0
		soc[i,3:4] = 0;

		# inherits neighbor's usage and beliefs.
		kk = rand()*2-1;
		kk = round(Int, kk/abs(kk));
			if(i == 1)			# Making sure I'm not at the edge
				kk = 1;
			elseif(i == length(soc[:,1]))
				kk = -1;
			end
		soc[i,1:2] = soc[i+kk,1:2];
	end
end

return soc

end

###########################################################

"""

    getrisk2(soc, cc, agelimit, irbw, ef, mfi, cpf)

Computes risk of exposure to cancer due to water contamination for each individual. please refer to this [link](http://pubs.acs.org/doi/abs/10.1021/es400316c) for more information.

Arguments:

- `soc` : `[belief, usage, cont-accumulation, age]` of every individual of the past year
- `cc` : contaminant concentration
- `agelimit` : age limit of each individual
- `irbw` : Ingestion rate of water per unit body weight (L kg/day)
- `ef` : Standard exposure frequency (day/year)
- `mfi`	: Metabolized fraction of contaminant
- `cpf` :  Cancer potency factor (Kg day/mg)

Returns:

- `risk` : Vector of risk (per thousand) for each individual.

"""

function getrisk2(soc, cc, agelimit, irbw, ef, mfi, cpf)
	age = soc[:, 4];
	accumulated_conc = soc[:, 3];

	estimated_usage = (agelimit .- age) .* cc + accumulated_conc; 	# Assume they will use same amount as this year their whole life.
	addin = estimated_usage .* 0.001 .* irbw .* ef ./ agelimit;
	risk = abs(1 - exp(-cpf .* float(addin) .* mfi));

return risk
end

#######################################################

"""

    updaterec(liverec, soc, cc, period, simyears, agelimit)

[deprecated] Updates the `soc` and `liverec` arrays containing social and contamination record `[belief, usage, contaminant-accumulation, time, age, simulation years]`
This method of storing records is used for risk assessment of people who already died, calculating risk given their lifetime exposure.

Arguments:

- `liverec` : live record of every individual in the past year
- `soc`: `[belief, usage, cont-accumulation, age]` of every individual of the past year
- `cc`: contaminant concentration
- `period` : coupling time increment (years)
- `simyears` : number of years individual has lived in the simulation
- `agelimit` : age limit of each individual

Returns:

- `liverec`
- `soc`

"""

function updaterec(liverec, soc, cc, period, simyears, agelimit)

soc[:,3] = soc[:,3] + !bool(soc[:,2]) .* cc;
soc[:,4] = soc[:,4] .+ period;
background_concentration = 5;

# for all those who die, a new person is born.
for i = 1:length(soc[:,1])
	if(soc[i,4] >= agelimit)
		at = agelimit;
		if(agelimit > simyears)
			at = simyears;
		end
		if(soc[i,3] > background_concentration)
			liverec = [liverec, [(soc[i,:])', at , simyears]'];
		end

		# For each person who dies, one is born, clean, age 0
		soc[i,3:4] = 0;
		kk = rand()*2-1;
		kk = kk/abs(kk);
			if(i == 1)		# Making sure I'm not at the edge
				kk = 1;
			elseif(i == length(soc[:,1]))
				kk = -1;
			end
		soc[i,1:2] = soc[i+kk,1:2];	# inherits neighbor's usage and beliefs.
	end
end

return liverec, soc

end

#################################################################

"""

    getrisk(liverec, irbw, ef, mfi, cpf)

[deprecated] Computes risk of exposure to cancer due to water contamination for each individual.

Arguments:

- `liverec` : `[belief, usage, concentration, cont-accumulation, simulation years, age]` of every individual of the past year
- `cc` : contaminant concentration
- `agelimit` : age limit of each individual
- `irbw` : Ingestion rate of water per unit body weight (L kg/day)
- `ef` : Standard exposure frequency (day/year)
- `mfi`	: Metabolized fraction of contaminant
- `cpf` :  Cancer potency factor (Kg day/mg)

Returns:

- `risk` : Vector of risk (per thousand) for each individual.

"""

function getrisk(liverec, irbw, ef, mfi, cpf)
	addin = liverec[:,3] .* 0.001 .* irbw .* ef ./ liverec[:,5];
	risk = abs(1 - exp(-cpf .* float(addin) .* mfi));

return risk
end

###############################################################
###############################################################
###############################################################

end

max_stars_repo_name: josephbakarji/socio-hydro-abm

max_issues_repo_name: josephbakarji/socio-hydro-abm

max_forks_repo_name: josephbakarji/socio-hydro-abm


text:
module SMARTBackendServices

import Dates
import HTTP
import HealthBase
import JSON3
import JSONWebTokens
import Random
import TimeZones
import URIs

const get_fhir_access_token = HealthBase.get_fhir_access_token

export BackendServicesConfig
export get_fhir_access_token
export backend_services

include("types.jl")

include("backend_services.jl")
include("jwt.jl")
include("timestamps.jl")

end # module


text:
include("master.jl")
include("child.jl")

#TODO do this correctly later
module LocalJobs
end


text:
"""
	struct EulerUpdater <: AbstractFieldUpdater

	EulerUpdater(; method =  :PECE, lam = 0.05)

A field updater implementing predictor-corrector Euler time stepping with static step size `lam`.
Available options for `method` include:

* PEC  - Predict-Evaluate-Correct scheme, with 1 density evaluation per step
* PECE - Predict-Evaluate-Correct-Evaluate scheme, with 2 density evaluations per step
"""
mutable struct EulerUpdater <: AbstractFieldUpdater
	method :: Symbol
	lam    :: Float64
	fields :: Dict{Int,FieldGrid{Float64}} # Used as pre-allocated storage for half-steps
	system :: Option{FieldSystem}

	function EulerUpdater(; method::Symbol = :PECE, lam::Real = 0.05)
		if !(method in (:PEC, :PECE))
			error("Invalid `method` parameter for EulerUpdater. Options include: (:PEC, :PECE)")
		end
		return new(method, lam, Dict(), nothing)
	end
end

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

Base.show(io::IO, updater::EulerUpdater) = @printf(io, "EulerUpdater(method = %s, lam = %.3g)", updater.method, updater.lam)

function setup!(updater::EulerUpdater, sys::FieldSystem)
	updater.system = sys
	copydict!(updater.fields, sys.fields)
	return nothing
end

function step!(updater::EulerUpdater)
	@assert !isnothing(updater.system)
	sys = updater.system

    # f(w_t): dH/dphi = dH_exc/dphi + eta - w
	# Time-stepping scheme:
	#	Predictor: wp_t+1 = w_t + dt * f(w_t)
	#	Corrector: w_t+1  = w_t + 0.5 * dt * (f(w_t) + f(wp_t+1))

	# 1) Predict
	for (mid, omega) in sys.fields
		@. omega += updater.lam * sys.residuals[mid]
		@. updater.fields[mid] = omega + 0.5 * updater.lam * sys.residuals[mid]
	end

	# 2 Evaluate
	residuals!(sys)

	# 3) Correct
	for (mid, temp) in updater.fields
		@. temp += 0.5 * updater.lam * sys.residuals[mid]
	end
	copydict!(sys.fields, updater.fields)

	# 4) Evaluate (PECE mode)
	if updater.method == :PECE
		residuals!(sys)
	end

	return nothing
end


text:
"""
`module Polys`

The exported symbols are
* `NBody`: an N-body function wrapped into a JuLIP calculator
* `blpoly_basis` : generate a basis of N-body functions
* `bapoly_basis` : generate a basis of N-body functions

## Usage

## Notation

* `N` : N-body
* `M` : number of edges, i.e., M = N (N-1)/2
* `K` : length of the tuples defining polynomials, K = M+1
"""
module Polys

import StaticPolynomials, NBodyIPs

using StaticArrays

using LinearAlgebra: dot

using NBodyIPs: NBodyFunction,
                bodyorder,
                _decode_dict,
                invariants,
                edges2bo,
                bo2edges,
                ninvariants,
                BASIS

import Base:              length,
                          Dict,
                          ==,
                          hash
import JuLIP:             cutoff
import JuLIP.Potentials:  @pot
import NBodyIPs:          fast,
                          degree,
                          combinebasis,
                          descriptor,
                          evaluate_many!,
                          evaluate_many_d!,
                          evaluate_I,
                          evaluate_I_ed,
                          basisname,
                          wtdegree

const Tup{M} = NTuple{M, Int}
const VecTup{M} = Vector{NTuple{M, Int}}

export NBPoly,
       StNBPoly,
       info   # TODO: move info to NBodyIPs and export from there? elsewhere?

# ==================================================================
#           Polynomials of Invariants
# ==================================================================

"""
`struct NBPoly`  (N-Body Polynomial, slow implementation)

A `struct` storing the information for a (pure) N-body potential, i.e.,
containing *only* terms of a specific body-order. Several `NBPoly`s can be
combined into an interatomic potential via `NBodyIP`.

### Fields

* `t::Vector{NTuple{M,TI}}` : list of M-tuples containing basis function information
e.g., if M = 7, α = t[n] is a 7-vector then this corresponds to the basis function
```
I2[α[7]] * ∏_{j=1..6} I1[j]^α[j]
```
where `I1, I2` are the 4-body invariants.

* `c::Vector{T}`: vector of coefficients for the basis functions

* `D`: a descriptor (cf `NBodyIPs.NBodyDescriptor`)
"""
struct NBPoly{N, M, T, TD} <: NBodyFunction{N, TD}
   t::VecTup{M}               # tuples M = #edges + 1
   c::Vector{T}               # coefficients
   D::TD                      # Descriptor
   valN::Val{N}               # encodes that this is an N-body function

   NBPoly(t::VecTup{M}, c::Vector{T}, D::TD, valN::Val{N}) where {N, M, T, TD} = (
      N <= 1 ? error("""NBPoly must have body-order 2 or larger;
                        use `NBodyIPs.OneBody{T}` for 1-body.""")
             : new{N, M, T, TD}(t, c, D, valN))
end

@pot NBPoly

==(V1::NBPoly, V2::NBPoly) = ( (V1.t == V2.t) && (V1.c == V2.c) && (V1.D == V2.D) )

descriptor(V::NBPoly) = V.D

basisname(::NBPoly) = "NBPoly"

hash(::BASIS, V::NBPoly) =
   hash( (hash(NBPoly), hash(bodyorder(V)), hash(BASIS(), V.D)) )

# standard constructor (N can be inferred)
NBPoly(t::VecTup{K}, c, D) where {K} = NBPoly(t, c, D, Val(edges2bo(K-1)))

# NBPoly made from a single basis function rather than a collection
NBPoly(t::Tup, c, D) = NBPoly([t], [c], D)

# collect multiple basis functions represented as NBPoly's into a single NBPoly
# (for performance reasons)
# TODO: this is not general enough!
NBPoly(B::Vector{TB}, c, D) where {TB <: NBPoly} =
      NBPoly([b.t[1] for b in B], c .* [b.c[1] for b in B], D)

# 1-body term (on-site energy)
NBPoly(c::Float64) = NBPoly([Tup{0}()], [c], nothing, Val(1))

# number of basis functions which this term is made from
length(V::NBPoly) = length(V.t)

cutoff(V::NBPoly) = cutoff(V.D)

function match_dictionary(V::NBPoly, V1::NBPoly)
   if V.D != V1.D
      if V.D.s != V1.D.s
         @warn("matching two non-matching dictionaries!")
      end
   end
   return NBPoly(V.t, V.c, V1.D, V.valN)
end

# combinebasis(basis::AbstractVector{TV}, coeffs) where {TV <: NBPoly} =
#       NBPoly(basis, coeffs, basis[1].D, basis[1].polytype)

function combinebasis(basis::AbstractVector{TV}, coeffs) where {TV <: NBPoly}
   # assume all basis functions have compatible descriptor
   # as well as compatible polytype, this should have been checked
   # already via the `hash(::BASIS,...)`

   # collect all tuples and coefficients into a long list
   tt = Vector{eltype(basis[1].t)}()
   cc = Vector{eltype(basis[1].c)}()
   for (b, c) in zip(basis, coeffs)
      append!(tt, b.t)
      append!(cc, c*b.c)
   end

   # sort `tt` and remember the ordering
   p = sortperm(tt)
   tt, cc = tt[p], cc[p]

   # compress into unique tuples
   t = [tt[1]]
   c = [cc[1]]
   for n = 2:length(tt)
      if tt[n] == t[end]
         c[end] += cc[n]
      else
         push!(t, tt[n])
         push!(c, cc[n])
      end
   end
   return NBPoly(t, c, basis[1].D) # , basis[1].polytype)
end

function degree(V::NBPoly)
   if length(V) == 1
      return NBodyIPs.PolyBasis.tdegree(descriptor(V), V.t[1])
   end
   error("`degree` is only defined for `NBPoly` basis functions, length == 1")
end

function wtdegree(V::NBPoly, rtF)
    if length(V) == 1
       return NBodyIPs.PolyBasis.wtdegree(descriptor(V), V.t[1], rtF)
    end
    error("`wtdegree` is only defined for `NBPoly` basis functions, length == 1")
end

function info(B::Vector{T}; indent = 2) where T <: NBPoly
   ind = repeat(" ", indent)
   println(ind * "body-order = $(bodyorder(B[1]))")
   println(ind * "    length = $(length(B))")
   if bodyorder(B[1]) > 1
      println(ind * " transform = $(B[1].D.s[1])")
      println(ind * "    cutoff = $(B[1].D.s[2])")
   end
end

include("fast_monomials.jl")

# ---------------  evaluate the n-body terms ------------------

# a tuple α = (α1, …, α6, α7) means the following:
# with f[0] = 1, f[1] = I7, …, f[5] = I11 we construct the basis set
#   f[α7] * g(I1, …, I6)
# this means, that gen_tuples must generate 7-tuples instead of 6-tuples
# with the first 6 entries restricted by degree and the 7th tuple must
# be in the range 0, …, 5

function evaluate_I(V::NBPoly, II)
   I1, I2 = II
   E = zero(eltype(I1))
   for (α, c) in zip(V.t, V.c)
      E += c * I2[1+α[end]] * monomial(α, I1)
   end
   return E
end

function evaluate_I_ed(V::NBPoly, II)
   I1, I2, dI1, dI2 = II
   E = zero(eltype(I1))
   dM = zero(typeof(I1))
   dE = zero(typeof(I1))
   #
   for (α, c) in zip(V.t, V.c)
      m, m_d = monomial_d(α, I1)
      E += c * I2[1+α[end]] * m        # just the value of the function itself
      dM += (c * I2[1+α[end]]) * m_d   # the I2 * ∇m term without the chain rule
      dE += (c * m) * dI2[1+α[end]]    # the ∇I2 * m term
   end
   # chain rule
   for i = 1:length(dI1)   # dI1' * dM
      dE += dM[i] * dI1[i]
   end
   return E, dE
end

# -------------- Infrastructure to read/write NBPoly  --------

Dict(V::NBPoly{N}) where {N} = Dict( "__id__" => "NBPoly",
                                      "t" => V.t,
                                      "c" => V.c,
                                      "D" => Dict(V.D),
                                      "N" => N )

NBPoly(D::Dict) = NBPoly([ tuple(ti...) for ti in D["t"] ],
                           Vector{Float64}(D["c"]),
                           _decode_dict(D["D"]),
                           Val(D["N"]))

Base.convert(::Val{:NBPoly}, D::Dict) = NBPoly(D)

# ==================================================================
#    StNBPoly
# ==================================================================

"""
`struct StNBPoly`  (N-Body Bond-length Polynomial)

fast evaluation of the outer polynomial using `StaticPolynomials`
"""
struct StNBPoly{N, TD, TP} <: NBodyFunction{N, TD}
   D::TD       # Descriptor
   P::TP       # a static polynomial
   valN::Val{N}
end

@pot StNBPoly

descriptor(V::StNBPoly) = V.D

function StNBPoly(V::NBPoly{N}) where {N}
   nI1, nI2 = ninvariants(V.D, N)
   nI = nI1 + nI2
   nmonomials = length(V.c)
   # generate the exponents for the StaticPolynomial
   exps = zeros(Int, nI, nmonomials)
   for (i, α) in enumerate(V.t)  # i = 1:nmonomials
      for (j, a) in enumerate(α[1:end-1])   #  ∏ I1[j]^α[j]
         exps[j, i] = a    # I1[j]^α[j]
      end
      exps[nI1+1+α[end], i] = 1   # I2[α[end]] * (...)
   end
   # generate the static polynomial
   return StNBPoly(V.D, StaticPolynomials.Polynomial(V.c, exps), V.valN)
end

cutoff(V::StNBPoly) = cutoff(V.D)

fast(Vn::StNBPoly)  = Vn
fast(Vn::NBPoly) =  StNBPoly(Vn)

evaluate_I(V::StNBPoly, II) =
      StaticPolynomials.evaluate(V.P, vcat(II...))

function evaluate_I_ed(V::StNBPoly, II)
   V, dV_dI = StaticPolynomials.evaluate_and_gradient(V.P, vcat(II[1], II[2]))

   # if length(dV_dI) != length(II[3]) + length(II[4])
   #    @show size.(II)
   #    @show size(dV_dI)
   #    @show size(vcat(II[1], II[2]))
   # end

   # TODO: check a few variants how to compute dV
   #       and test performance
   II34 = vcat(II[3], II[4])
   dV = II34[1]*dV_dI[1]
   for i = 2:length(II34)
      dV += II34[i] * dV_dI[i]
   end

   return V, dV # sum(i * dV_di for (i, dVdi) in zip( .* dV_dI)  # (dI' * dV_dI)
end

end # module

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
using ProximalOperators
using ProximalAlgorithms
using AbstractOperators
using LinearAlgebra
using SparseArrays
using FFTW
using Random
using Test
using DelimitedFiles

@testset "ProximalAlgorithms" begin

@testset "Utilities" begin
    include("test_conjugate.jl")
end

@testset "Algorithms" begin
    include("test_template.jl")
    include("test_lasso_small.jl")
    include("test_lasso_harder.jl")
    include("test_lasso_small_split_x.jl")
    include("test_lasso_small_split_f.jl")
    include("test_l1logreg_small.jl")
    include("test_afba.jl")
    include("test_afba_LP.jl")
end

end


text:
module FefferyUtilsComponents
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.0.5"

include("jl/''_fefferycaptcha.jl")
include("jl/''_fefferypasteimage.jl")
include("jl/''_fefferysyntaxhighlighter.jl")
include("jl/''_fefferywatermark.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "feffery_utils_components",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "feffery_utils_components.min.js",
    external_url = "https://unpkg.com/feffery_utils_components@0.0.5/feffery_utils_components/feffery_utils_components.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "feffery_utils_components.min.js.map",
    external_url = "https://unpkg.com/feffery_utils_components@0.0.5/feffery_utils_components/feffery_utils_components.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
)
            ]
        )

    )
end
end

max_stars_repo_name: CNFeffery/feffery-utils-components

max_issues_repo_name: CNFeffery/feffery-utils-components

max_forks_repo_name: CNFeffery/feffery-utils-components


text:
function hex_escape(c::Char)::String
    tchar = first(transcode(UInt8, string(c)))
    string("\\x", string(tchar; pad = sizeof(tchar)<<1, base = 16))
end

hex_escape(s::String)::String = join([hex_escape(c) for c in s])

hex_escape(s) = hex_escape(string(s))


text:
include("epamodule.jl")
em = Main.epamodule
rede_path = "/home/augusto/Documents/IC-2020/optimized-calibration/networks/c-town/C-Town.inp"

em.ENopen(rede_path)
em.ENopenH()

link = "P2"
rugosidade = 0.001
em.ENsetlinkvalue(em.ENgetlinkindex(link), em.EN_ROUGHNESS, rugosidade)
em.ENsolveH()
println("Teste tubulações:")
println(em.ENgetlinkvalue(em.ENgetlinkindex(link), em.EN_ROUGHNESS))
println(rugosidade)

node = "J414"
em.ENsolveH()
println("Teste nós:")
println(em.ENgetnodevalue(em.ENgetnodeindex(node), em.EN_PRESSURE))

em.ENcloseH()
em.ENclose()

max_stars_repo_name: augusnunes/epanet-julia

max_issues_repo_name: augusnunes/epanet-julia

max_forks_repo_name: augusnunes/epanet-julia


text:
module TransitionMixtureTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: SPTransitionMixtureOutNCCPX, SPTransitionMixtureIn1CNCPX, SPTransitionMixtureZCCNPX, SVBTransitionMixtureOutNCCDX, SVBTransitionMixtureIn1CNCDX, SVBTransitionMixtureZCCNDX,  SVBTransitionMixtureA, MTransitionMixtureCCCDX, MTransitionMixtureCCCNX

#-------------
# Update rules
#-------------

@testset "SPTransitionMixtureOutNCCPX" begin
    @test SPTransitionMixtureOutNCCPX <: SumProductRule{TransitionMixture}
    @test outboundType(SPTransitionMixtureOutNCCPX) == Message{Categorical}
    @test isApplicable(SPTransitionMixtureOutNCCPX, [Nothing, Message{Categorical}, Message{Categorical}, Message{PointMass}, Message{PointMass}])

    @test ruleSPTransitionMixtureOutNCCPX(nothing, Message(Univariate, Categorical, p=[0.1, 0.9]), Message(Univariate, Categorical, p=[1.0, 0.0]), Message(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), Message(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(Univariate, Categorical, p=[0.19000000000000003, 0.81])
end

@testset "SPTransitionMixtureIn1CNCPX" begin
    @test SPTransitionMixtureIn1CNCPX <: SumProductRule{TransitionMixture}
    @test outboundType(SPTransitionMixtureIn1CNCPX) == Message{Categorical}
    @test isApplicable(SPTransitionMixtureIn1CNCPX, [Message{Categorical}, Nothing, Message{Categorical}, Message{PointMass}, Message{PointMass}])

    @test ruleSPTransitionMixtureIn1CNCPX(Message(Univariate, Categorical, p=[0.2, 0.8]), nothing, Message(Univariate, Categorical, p=[1.0, 0.0]), Message(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), Message(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(Univariate, Categorical, p=[0.5211267605633803, 0.4788732394366197])
end

@testset "SPTransitionMixtureZCCNPX" begin
    @test SPTransitionMixtureZCCNPX <: SumProductRule{TransitionMixture}
    @test outboundType(SPTransitionMixtureZCCNPX) == Message{Categorical}
    @test isApplicable(SPTransitionMixtureZCCNPX, [Message{Categorical}, Message{Categorical}, Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPTransitionMixtureZCCNPX(Message(Univariate, Categorical, p=[0.2, 0.8]), Message(Univariate, Categorical, p=[1.0, 0.0]), nothing, Message(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), Message(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(Univariate, Categorical, p=[0.5441176470588237, 0.4558823529411764])
end

@testset "SVBTransitionMixtureOutNCCDX" begin
    @test SVBTransitionMixtureOutNCCDX <: StructuredVariationalRule{TransitionMixture}
    @test outboundType(SVBTransitionMixtureOutNCCDX) == Message{Categorical}
    @test isApplicable(SVBTransitionMixtureOutNCCDX, [Nothing, Message{Categorical}, Message{Categorical}, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleSVBTransitionMixtureOutNCCDX(nothing, Message(Univariate, Categorical, p=[0.1, 0.9]), Message(Univariate, Categorical, p=[0.3, 0.7]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(Univariate, Categorical, p=[0.33000000000033997, 0.6699999999996599])
end

@testset "SVBTransitionMixtureIn1CNCDX" begin
    @test SVBTransitionMixtureIn1CNCDX <: StructuredVariationalRule{TransitionMixture}
    @test outboundType(SVBTransitionMixtureIn1CNCDX) == Message{Categorical}
    @test isApplicable(SVBTransitionMixtureIn1CNCDX, [Message{Categorical}, Nothing, Message{Categorical}, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleSVBTransitionMixtureIn1CNCDX(Message(Univariate, Categorical, p=[0.2, 0.8]), nothing, Message(Univariate, Categorical, p=[0.3, 0.7]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(Univariate, Categorical, p=[0.5239616613418148, 0.4760383386581853])
end

@testset "SVBTransitionMixtureZCCNDX" begin
    @test SVBTransitionMixtureZCCNDX <: StructuredVariationalRule{TransitionMixture}
    @test outboundType(SVBTransitionMixtureZCCNDX) == Message{Categorical}
    @test isApplicable(SVBTransitionMixtureZCCNDX, [Message{Categorical}, Message{Categorical}, Nothing, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleSVBTransitionMixtureZCCNDX(Message(Univariate, Categorical, p=[0.2, 0.8]), Message(Univariate, Categorical, p=[0.1, 0.9]), nothing, ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(Univariate, Categorical, p=[0.5479233226837061, 0.4520766773162939])
end

@testset "SVBTransitionMixtureA" begin
    @test SVBTransitionMixtureA <: StructuredVariationalRule{TransitionMixture}
    @test outboundType(SVBTransitionMixtureA) == Message{Dirichlet}
    @test isApplicable(SVBTransitionMixtureA, [ProbabilityDistribution, Nothing, ProbabilityDistribution])

    @test ruleSVBTransitionMixtureA(ProbabilityDistribution(Multivariate, Contingency, p=[0.3*[0.1 0.6; 0.2 0.1], 0.7*[0.1 0.6; 0.2 0.1]]), nothing, ProbabilityDistribution(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == Message(MatrixVariate, Dirichlet, a=[1.03 1.18; 1.06 1.03])
end

@testset "MTransitionMixtureCCCDX" begin
    @test MTransitionMixtureCCCDX <: MarginalRule{TransitionMixture}
    @test isApplicable(MTransitionMixtureCCCDX, [Message{Categorical}, Message{Categorical}, Message{Categorical}, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleMTransitionMixtureCCCDX(Message(Univariate, Categorical, p=[0.2, 0.8]), Message(Univariate, Categorical, p=[0.1, 0.9]), Message(Univariate, Categorical, p=[0.3, 0.7]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == ProbabilityDistribution(Multivariate, Contingency, p=[[0.0009966777408637877 0.017940199335548173; 0.035880398671096346 0.28704318936877077], [0.00697674418604651 0.08372093023255814; 0.06511627906976744 0.5023255813953489]])
end

@testset "MTransitionMixtureCCCNX" begin
    @test MTransitionMixtureCCCNX <: MarginalRule{TransitionMixture}
    @test isApplicable(MTransitionMixtureCCCNX, [Message{Categorical}, Message{Categorical}, Message{Categorical}, Nothing, Nothing])

    @test ruleMTransitionMixtureCCCNX(Message(Univariate, Categorical, p=[0.2, 0.8]), Message(Univariate, Categorical, p=[0.1, 0.9]), Message(Univariate, Categorical, p=[0.3, 0.7]), Message(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), Message(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == ProbabilityDistribution(Multivariate, Contingency, p=[[0.0009966777408637875 0.017940199335548173; 0.035880398671096346 0.28704318936877077], [0.0069767441860465115 0.08372093023255814; 0.06511627906976744 0.5023255813953489]])
end

@testset "averageEnergy" begin
    @test averageEnergy(TransitionMixture, ProbabilityDistribution(Multivariate, Contingency, p=[0.3*[0.1 0.6; 0.2 0.1], 0.7*[0.1 0.6; 0.2 0.1]]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.1 0.2; 0.9 0.8]), ProbabilityDistribution(MatrixVariate, PointMass, m=[0.3 0.4; 0.7 0.6])) == 0.9266048040118577
end

end #module


text:
# kruskal_wallis.jl
# Kruskal-Wallis rank sum test
#
# Copyright (C) 2014   Christoph Sawade
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export KruskalWallisTest

struct KruskalWallisTest <: HypothesisTest
    n_i::Vector{Int}         # number of observations in each group
    df::Int                  # degrees of freedom
    R_i::Vector{Float64}     # rank sums
    H::Float64               # test statistic: chi-square statistic
    tie_adjustment::Float64  # adjustment for ties
end

"""
    KruskalWallisTest(groups::AbstractVector{<:Real}...)

Perform Kruskal-Wallis rank sum test of the null hypothesis that the `groups`
``\\mathcal{G}`` come from the same distribution against the alternative hypothesis that
that at least one group stochastically dominates one other group.

The Kruskal-Wallis test is an extension of the Mann-Whitney U test to more than two groups.

The p-value is computed using a ``χ^2`` approximation to the distribution of the test
statistic ``H_c=\\frac{H}{C}``:
```math
    \\begin{align*}
    H & = \\frac{12}{n(n+1)} \\sum_{g ∈ \\mathcal{G}} \\frac{R_g^2}{n_g} - 3(n+1)\\\\
    C & = 1-\\frac{1}{n^3-n}\\sum_{t ∈ \\mathcal{T}} (t^3-t),
    \\end{align*}
```
where ``\\mathcal{T}`` is the set of the counts of tied values at each tied position, ``n``
is the total number of observations across all groups, and ``n_g`` and ``R_g`` are the number of
observations and the rank sum in group ``g``, respectively. See references for further
details.

Implements: [`pvalue`](@ref)

# References

  * Meyer, J.P, Seaman, M.A., Expanded tables of critical values for the Kruskal-Wallis
    H statistic. Paper presented at the annual meeting of the American Educational Research
    Association, San Francisco, April 2006.

# External links

  * [Kruskal-Wallis test on Wikipedia
    ](https://en.wikipedia.org/wiki/Kruskal-Wallis_one-way_analysis_of_variance)
"""
function KruskalWallisTest(groups::AbstractVector{T}...) where T<:Real
    (H, R_i, tieadj, n_i) = kwstats(groups...)
    if length(groups)<=3 && any(n_i .< 6)
        @warn("This test is only asymptotically correct and might be inaccurate for the given group size")
    end
    df = length(groups) - 1
    KruskalWallisTest(n_i, df, R_i, H, tieadj)
end

testname(::KruskalWallisTest) = "Kruskal-Wallis rank sum test (chi-square approximation)"
population_param_of_interest(x::KruskalWallisTest) = ("Location parameters", "all equal", NaN) # parameter of interest: name, value under h0, point estimate
default_tail(test::KruskalWallisTest) = :right

function show_params(io::IO, x::KruskalWallisTest, ident)
    println(io, ident, "number of observation in each group: ", x.n_i)
    println(io, ident, "χ²-statistic:                        ", x.H)
    println(io, ident, "rank sums:                           ", x.R_i)
    println(io, ident, "degrees of freedom:                  ", x.df)
    println(io, ident, "adjustment for ties:                 ", x.tie_adjustment)
end

pvalue(x::KruskalWallisTest) = pvalue(Chisq(x.df), x.H; tail=:right)

## helper

# Get H, rank sums, and tie adjustment for Kruskal-Wallis test
function kwstats(groups::AbstractVector{T}...) where T<:Real
    n_i = [length(g) for g in groups]
    n = sum(n_i)

    # get ranks and adjustment for ties
    (ranks, tieadj) = tiedrank_adj([groups...;])
    C = 1-tieadj/(n^3 - n)

    # compute rank sums
    R_i = Vector{Float64}(undef, length(groups))
    n_end = 0
    for i in 1:length(groups)
        R_i[i] = sum(ranks[n_end+1:n_end+n_i[i]])
        n_end += n_i[i]
    end

    # compute test statistic and correct for ties
    H = 12 * sum(R_i .^ 2 ./ n_i) / (n * (n + 1)) - 3 * (n + 1)
    H /= C

    (H, R_i, C, n_i)
end


text:
using LazySets: _isapprox

for N in [Float64, Float32, Rational{Int}]
    # random polygons
    rand(HPolygon)
    rand(HPolygonOpt)
    rand(VPolygon)

    # Empty polygon
    p = HPolygon{N}()

    # Test that constraints are automatically sorted
    c1 = LinearConstraint(N[2, 2], N(12))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(0))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p, c3)
    addconstraint!(p, c1)
    addconstraint!(p, c4)
    addconstraint!(p, c2)
    @test p.constraints[1] == c1
    @test p.constraints[2] == c2
    @test p.constraints[3] == c3
    @test p.constraints[4] == c4

    # redundancy of constraints
    h1 = HalfSpace([N(1), N(1)], N(1))
    h2 = HalfSpace([N(0), N(1)], N(1)) # redundant
    h3 = HalfSpace([N(-1), N(1)], N(1))
    h4 = HalfSpace([N(0), N(-1)], N(0))
    @test LazySets.isredundant(h2, h1, h3)
    @test !LazySets.isredundant(h3, h2, h4)
    @test LazySets.isredundant(h2, h2, h3) && LazySets.isredundant(h3, h2, h3)
    p2 = HPolygon([h1, h2, h3, h4]; sort_constraints=false) # sorted in the right order
    @test length(p2.constraints) == 4
    remove_redundant_constraints!(p2)
    @test length(p2.constraints) == 3
    h1 = HalfSpace([N(1), N(0)], N(1))
    h2 = HalfSpace([N(0), N(1)], N(1))
    h3 = HalfSpace([N(-1), N(0)], N(1))
    h4 = HalfSpace([N(0), N(-1)], N(1))
    p2 = HPolygon{N}()
    addconstraint!(p2, h1)
    addconstraint!(p2, h2)
    addconstraint!(p2, h3)
    addconstraint!(p2, h4)
    addconstraint!(p2, h1)
    @test ispermutation(p2.constraints, [h1, h2, h3, h4])
    h5 = HalfSpace([N(1), N(1)], N(0))
    addconstraint!(p2, h5)
    @test ispermutation(p2.constraints, [h3, h4, h5])

    # constructor from simple H-rep
    A = N[2 0; 1 3]
    b = N[-1, 1]
    @test constraints_list(HPolygon(A, b)) ==
          constraints_list(HPolygonOpt(A, b)) ==
          [HalfSpace(N[2, 0], N(-1)), HalfSpace(N[1, 3], N(1))]

    # conversion to optimized polygon
    po = convert(HPolygonOpt, p)
    # conversion back
    convert(HPolygon, po)
    # conversion from HPolytope
    polytope = HPolytope{N}()
    addconstraint!(polytope, c1)
    addconstraint!(polytope, c2)
    addconstraint!(polytope, c3)
    addconstraint!(polytope, c4)
    convert(HPolygon, polytope)
    convert(HPolygonOpt, polytope)
    # conversion to HPolytope
    HPolytope(constraints_list(p))
    HPolytope(constraints_list(po))

    # conversion from other set type
    H = Hyperrectangle(low=N[-1, -1], high=N[1, 1])
    HPolygon(constraints_list(H))
    HPolygonOpt(constraints_list(H))

    # support vector of polygon with no constraints
    @test_throws AssertionError σ(N[0], HPolygon{N}())
    @test_throws AssertionError σ(N[0], HPolygonOpt{N}())

    # boundedness
    @test isbounded(p) && isbounded(po)
    @test !isbounded(HPolygon{N}(), false) &&
          !isbounded(HPolygonOpt{N}(), false)
    @test_throws AssertionError HPolygon(LinearConstraint{N, Vector{N}}[];
                                        check_boundedness=true)
    @test_throws AssertionError HPolygonOpt(LinearConstraint{N, Vector{N}}[];
                                           check_boundedness=true)

    # isempty
    @test !isempty(p)
    @test !isempty(po)

    # subset
    b1 = Ball1(N[0, 0], N(1))
    b2 = Ball1(N[1, 1], N(4))
    l1 = LinearMap(N[1 0; 0 1], b1)
    l2 = LinearMap(N[1 0; 0 1], b2)
    subset, point = ⊆(p, b1, true)
    @test !subset && point ∈ p && point ∉ b1
    subset, point = ⊆(p, l1, true)
    @test !subset && point ∈ p && point ∉ l1
    subset, point = ⊆(p, b2, true)
    @test subset && p ⊆ b2 && point == N[]
    subset, point = ⊆(p, l2, true)
    @test subset && p ⊆ l2 && point == N[]
    @test p ⊆ p

    # HPolygon/HPolygonOpt tests
    for (hp, t_hp) in [(p, HPolygon), (po, HPolygonOpt)]
        # Test Dimension
        @test dim(hp) == 2

        # test support vector, with linear and binary search
        d = N[1, 0]
        @test σ(d, hp) == N[4, 2]
        @test σ(d, hp, linear_search=true) == σ(d, hp, linear_search=false)
        d = N[0, 1]
        @test σ(d, hp) == N[2, 4]
        @test σ(d, hp, linear_search=true) == σ(d, hp, linear_search=false)
        d = N[0, -1]
        @test σ(d, hp) == N[0, 0]
        @test σ(d, hp, linear_search=true) == σ(d, hp, linear_search=false)
        d = N[-1, 0]
        @test σ(d, hp) == N[-1, 1]
        @test σ(d, hp, linear_search=true) == σ(d, hp, linear_search=false)
        d = N[1, -1]
        @test σ(d, hp) == N[4, 2]
        @test σ(d, hp, linear_search=true) == σ(d, hp, linear_search=false)

        # membership
        @test N[0, 0] ∈ hp
        @test N[4, 2] ∈ hp
        @test N[2, 4] ∈ hp
        @test N[-1, 1] ∈ hp
        @test N[2, 3] ∈ hp
        @test N[1, 1] ∈ hp
        @test N[3, 2] ∈ hp
        @test N[5 / 4, 7 / 4] ∈ hp
        @test N[4, 1] ∉ hp
        @test N[5, 2] ∉ hp
        @test N[3, 4] ∉ hp
        @test N[-1, 5] ∉ hp

        # an_element function
        @test an_element(hp) ∈ hp
        hp_shallow = HPolygon{N}()
        @test_throws AssertionError an_element(hp_shallow)
        addconstraint!(hp_shallow, c1)
        @test_throws AssertionError an_element(hp_shallow)

        # isuniversal
        answer, w = isuniversal(hp, true)
        @test !isuniversal(hp) && !answer && w ∉ hp

        # hrep conversion
        @test tohrep(hp) == hp

        # constraints list getter
        @test constraints_list(hp) == hp.constraints

        # test constraints list of a VPolygon
        vp = tovrep(hp)
        @test ispermutation(constraints_list(vp), hp.constraints)

        # translation
        @test translate(hp, N[1, 2]) == t_hp(
            [HalfSpace(N[2, 2], N(18)), HalfSpace(N[-3, 3], N(9)),
             HalfSpace(N[-1, -1], N(-3)), HalfSpace(N[2, -4], N(-6))])

        # test for concrete minkowski sum
        A = [N[4, 0], N[6, 2], N[4, 4]]
        B = [N[-2, -2], N[2, 0], N[2, 2], N[-2, 4]]

        P = VPolygon(A)
        Q = VPolygon(B)
        PQ = minkowski_sum(P, Q)
        @test is_cyclic_permutation(PQ.vertices, [N[2, -2], N[6, 0], N[8, 2],
                                                N[8, 4], N[6, 6], N[2, 8]])

        # test for different starting points in vertices_list of minkowski sum
        P = VPolygon([N[4, 0], N[6, 2], N[4, 4]])
        P2 = VPolygon([N[4, 4], N[4, 0], N[6, 2]])
        Q = VPolygon([N[-2, -2], N[2, 0], N[2, 2], N[-2, 4]])
        @test is_cyclic_permutation(vertices_list(minkowski_sum(P, Q)), vertices_list(minkowski_sum(P2, Q)))

        # test for corner case of parallel edges in minkowski sum
        C = [N[10, 5], N[10, 10], N[8, 10], N[5, 8], N[5, 5]]
        D = [N[-1, -2], N[1, -2], N[-1, 2]]

        R = VPolygon(C)
        S = VPolygon(D)
        RS = minkowski_sum(R, S)
        @test is_cyclic_permutation(RS.vertices, [N[4, 3], N[11, 3], N[11, 8],
                                                N[9,12], N[7, 12], N[4, 10]])
    end

    # concrete intersection of H-rep
    p2 = HPolygon{N}()
    c1 = LinearConstraint(N[2, 2], N(14))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(-2))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p2, c3)
    addconstraint!(p2, c1)
    addconstraint!(p2, c4)
    addconstraint!(p2, c2)
    p3 = intersection(p, p2)
    @test length(constraints_list(p3)) == 4

    # check that tighter constraints are used in intersection (#883)
    h1 = HalfSpace([N(1), N(0)], N(3))
    h2 = HalfSpace([N(-1), N(0)], N(-3))
    h3 = HalfSpace([N(0), N(1)], N(7))
    h4 = HalfSpace([N(0), N(-1)], N(-5))
    h5 = HalfSpace([N(0), N(1)], N(6))
    h6 = HalfSpace([N(0), N(-1)], N(-4))
    p1 = HPolygon([h1, h3, h2, h4])
    p2 = HPolygon([h1, h5, h2, h6])
    c = intersection(p1, p2).constraints
    @test c == [h1, h5, h2, h4]
    h1 = HalfSpace([N(1), N(1)], N(1))
    h2 = HalfSpace([N(-1), N(1)], N(1))
    h3 = HalfSpace([N(0), N(-1)], N(0))
    p1 = HPolygon([h1, h2, h3])
    h4 = HalfSpace([N(0), N(1)], N(1))
    h5 = HalfSpace([N(-1), N(-1)], N(0))
    h6 = HalfSpace([N(1), N(-1)], N(0))
    p2 = HPolygon([h4, h5, h6])
    c = intersection(p1, p2).constraints
    @test c == [h1, h2, h5, h6]

    # check that empty polygon (infeasible constraints) has no vertices (#918)
    P = HPolygon([HalfSpace(N[1, 1], N(0)), HalfSpace(N[-1, 0], N(-1)),
                  HalfSpace(N[0, -1], N(-1))])
    @test vertices_list(P) == Vector{Vector{N}}()
    P = HPolygon([HalfSpace(N[1, 0], N(0)), HalfSpace(N[-1, 0], N(-1))])
    @test vertices_list(P) == Vector{Vector{N}}()
    @test isempty(P)

    # empty intersection results in empty set
    p3 = tohrep(VPolygon([N[0, 0]]))
    p4 = tohrep(VPolygon([N[1, 1]]))
    @test intersection(p3, p4) isa EmptySet{N}

    # concrete linear map
    # in 2D and for an invertible map we get an HPolygon (see #631 and #1093)
    HP = convert(HPolygon, BallInf(N[0, 0], N(1)))
    @test linear_map(N[1 0; 0 2], HP) isa HPolygon{N}
    # in higher dimensions we get an HPolytope (#2168)
    if test_suite_polyhedra
        p4 = convert(HPolygon, BallInf(zeros(N, 2), N(1)))
        A = ones(N, 4, 2)
        @test linear_map(A, p4) isa HPolytope
    end

    # vertices_list removes duplicates by default (#1405)
    p3 = HPolygon([HalfSpace(N[1, 0], N(0)), HalfSpace(N[0, 1], N(0)),
                   HalfSpace(N[-1, 0], N(0)), HalfSpace(N[0, -1], N(0))])
    @test vertices_list(p3) == [N[0, 0]]

    # Test VRepresentation
    vp = tovrep(p)
    @test N[2, 4] ∈ vertices_list(vp)
    @test N[-1, 1] ∈ vertices_list(vp)
    @test N[0, 0] ∈ vertices_list(vp)
    @test N[4, 2] ∈ vertices_list(vp)
    @test tovrep(vp) == vp

    # test convex hull of a set of points using the default algorithm
    v1 = N[1//10, 3//10]
    v2 = N[1//5, 1//10]
    v3 = N[2//5, 3//5]
    v4 = N[9//10, 1//5]
    v5 = N[3//10, 7//25]
    points = [v4, v3, v2, v1, v5]
    vp = VPolygon(points) # by default, a convex hull is run
    @test ispermutation(vertices_list(vp), [v1, v2, v3, v4])

    vp = VPolygon(points, apply_convex_hull=false) # we can turn it off
    @test ispermutation(vertices_list(vp), [v1, v2, v3, v4, v5])

    # test for pre-sorted points
    v3_new = N[2//5, 3//10]
    points = [v1, v2, v3_new, v3, v4]
    vp = VPolygon(points, algorithm="monotone_chain_sorted")
    @test ispermutation(vertices_list(vp), [v1, v2, v3, v4])

    # test that #83 is fixed
    p = VPolygon([N[2, 3]])
    @test N[2, 3] ∈ p
    @test N[3, 2] ∉ p
    v = N[-1, -17//5]
    p = VPolygon([N[2, 3], v])
    @test v ∈ p

    # an_element function
    p = VPolygon([N[2, 3]])
    @test an_element(p) ∈ p

    # membership
    point = N[0, 1]
    @test point ∉ VPolygon(Vector{Vector{N}}(), apply_convex_hull=false)
    @test point ∉ VPolygon([N[0, 2]])
    @test point ∈ VPolygon([N[0, 0], N[0, 2]])
    @test point ∉ VPolygon([N[1, 0], N[1, 2]])

    # translation
    vp2 = VPolygon([N[0, 0], N[1, 0], N[0, 1]])
    @test translate(vp2, N[1, 2]) == VPolygon([N[1, 2], N[2, 2], N[1, 3]])

    # subset
    p1 = VPolygon([N[0, 0], N[2, 0]])
    p2 = VPolygon([N[1, 0]])
    b = BallInf(N[2, 0], N(1))
    @test p2 ⊆ p1 && ⊆(p2, p1, true)[1]
    subset, witness = ⊆(p1, b, true)
    @test p1 ⊈ b && !subset && witness ∈ p1 && witness ∉ b
    subset, witness = ⊆(p2, b, true)
    @test p2 ⊆ b && subset

    v1 = N[1, 0]
    v2 = N[1, 1]
    v3 = N[0, 1]
    v4 = N[-1, 1]
    v5 = N[-1, 0]
    v6 = N[-1, -1]
    v7 = N[0, -1]
    v8 = N[1, -1]
    v = [v1, v2, v3, v4, v5, v6, v7, v8]

    for (i, vi) in enumerate(v)
        for j in i:8
                @test vi <= v[j]
        end
    end

    # hrep conversion
    v1 = N[9//10, 1//5]
    v2 = N[2//5, 3//5]
    v3 = N[1//5, 1//10]
    v4 = N[1//10, 3//10]
    v5 = N[7//10, 2//5]
    v6 = N[4//5, 3//10]
    v7 = N[3//10, 11//20]
    v8 = N[1//5, 9//20]
    points5 = [v1, v2, v3, v4, v5, v6, v7, v8]
    for i in [0, 1, 2, 4, 8]
        points = i == 0 ? Vector{Vector{N}}() : points5[1:i]
        vp = VPolygon(points, apply_convex_hull=i > 0)
        h1 = tohrep(vp)
        @test convert(HPolygon, vp) == h1
        if i == 0
            @test isempty(h1)
            continue
        elseif i == 1
            @test v1 ∈ h1
        elseif i == 2
            c = h1.constraints[1]
            @test c.a ≈ N[2//5, 1//2] && c.b ≈ N(23//50)
            c = h1.constraints[3]
            @test c.a ≈ N[-2//5, -1//2] && c.b ≈ N(-23//50)
        elseif i == 4
            @test length(h1.constraints) == 4
            c = h1.constraints[1]
            @test c.a ≈ N[2//5, 1//2] && c.b ≈ N(23//50)
            c = h1.constraints[2]
            @test c.a ≈ N[-3//10, 3//10] && c.b ≈ N(3//50)
            c = h1.constraints[3]
            @test c.a ≈ N[-1//5, -1//10] && c.b ≈ N(-1//20)
            c = h1.constraints[4]
            @test c.a ≈ N[1//10, -7//10] && c.b ≈ N(-1//20)
        end

        # check that constraints are sorted correctly
        h2 = HPolygon{N}()
        for c in h1.constraints
            addconstraint!(h2, c)
        end
        @test h1.constraints == h2.constraints
    end

    # empty VPolygon: conversion to hrep
    @test tohrep(VPolygon{N}()) isa EmptySet{N}

    # test VPolygon constructor given the matrix of vertices
    m = N[4 0; 6 2; 4 4]'
    P = VPolygon(m)
    @test is_cyclic_permutation(vertices_list(P), [N[4, 0], N[6, 2], N[4, 4]])
end

function same_constraints(v::Vector{<:LinearConstraint{N}})::Bool where N<:Real
    c1 = v[1]
    for k = 2:length(v)
        c2 = v[2]
        if c1.a != c2.a || c1.b != c2.b
            return false
        end
    end
    return true
end

for N in [Float64, Float32]
    v1 = N[1//10, 3//10]
    v2 = N[1//5, 1//10]
    v3 = N[2//5, 3//10]
    v4 = N[2//5, 3//5]
    v5 = N[9//10, 1//5]

    # test support vector of a VPolygon
    points = [v1, v2, v3, v4, v5]
    vp = VPolygon(points, algorithm="monotone_chain_sorted")
    d = N[1, 0]
    @test σ(d, vp) == points[5]
    d = N[0, 1]
    @test σ(d, vp) == points[4]
    d = N[-1, 0]
    @test σ(d, vp) == points[1]
    d = N[0, -1]
    @test σ(d, vp) == points[2]
    dirs = [N[1, 0], N[1, 1//2], N[1//2, 1//2], N[1//2, 1], N[0, 1],
            N[-1//2, 1], N[-1//2, 1//2], N[-1, 1//2], N[-1, 0], N[1, -1//2],
            N[1//2, -1//2], N[1//2, -1], N[0, -1], N[-1//2, -1],
            N[-1//2, -1//2], N[-1, -1//2]]
    B = Ball2(zeros(N, 2), N(1))
    P = HPolygon([HalfSpace(di, ρ(di, B)) for di in dirs])
    vlistP = vertices_list(P)
    V = VPolygon(vlistP)
    all(x -> _isapprox(σ(x, V), x), vlistP)

    # test adding constraints, with linear and binary search
    p1 = HPolygon{N}()
    p2 = HPolygon{N}()
    po1 = HPolygonOpt{N}()
    po2 = HPolygonOpt{N}()
    n = 10
    for i in 1:n
        constraint = LinearConstraint(rand(N, 2), rand(N))
        addconstraint!(p1, constraint, linear_search=true)
        addconstraint!(p2, constraint, linear_search=i<=2)
        addconstraint!(po1, constraint, linear_search=true)
        addconstraint!(po2, constraint, linear_search=i<=2)
    end
    n = length(p1.constraints)
    @test n == length(p2.constraints) == length(po2.constraints) ==
          length(po2.constraints)
    for i in 1:n
        @test same_constraints([p1.constraints[i], p2.constraints[i],
                                po1.constraints[i], po2.constraints[i]])
    end

    # check redundancy removal
    p2 = HPolygon([
        HalfSpace(N[-1.29817, 1.04012], N(6.07731)),
        HalfSpace(N[-1.29348, -0.0920708], N(1.89515)),
        HalfSpace(N[-1.42909, -0.347449], N(1.50577)),
        HalfSpace(N[0.349446, -0.130477], N(-0.575904)),
        HalfSpace(N[1.50108, -0.339291], N(-2.18331)),
        HalfSpace(N[2.17022, -0.130831], N(-2.14411))])
    addconstraint!(p2, p2.constraints[2])
    @test length(p2.constraints) == 6
end

for N in [Float64]
    # check boundedness after conversion
    H = Hyperrectangle(low=N[-1, -1], high=N[1, 1])
    HPolygon(constraints_list(H); check_boundedness=true)
    HPolygonOpt(constraints_list(H); check_boundedness=true)

    p = HPolygon{N}()
    c1 = LinearConstraint(N[2, 2], N(12))
    c2 = LinearConstraint(N[-3, 3], N(6))
    c3 = LinearConstraint(N[-1, -1], N(0))
    c4 = LinearConstraint(N[2, -4], N(0))
    addconstraint!(p, c3)
    addconstraint!(p, c1)
    addconstraint!(p, c4)
    addconstraint!(p, c2)
    po = convert(HPolygonOpt, p)

    # test boundedness
    @test isbounded(p, false) && isbounded(po, false)

    # is intersection empty
    p3 = convert(HPolygon, Hyperrectangle(low=N[-1, -1], high=N[1, 1]))
    I1 = Interval(N(9//10), N(11//10))
    I2 = Interval(N(1//5), N(3//10))
    I3 = Interval(N(4), N(5))
    @test !is_intersection_empty(I1 × I2 , p3)
    @test is_intersection_empty(I1 × I3 , p3)

    # redundancy with almost-parallel constraints (numeric issues, #2102)
    Hs = [HalfSpace([0.05125918978413134, 0.0], 0.956012965730266),
          HalfSpace([-0.0363879863894044, 0.3036279569180349], -0.370213880061176),
          HalfSpace([0.03529692609786872, -0.2945239519959644], 0.36249837124872736)]
    H = HalfSpace([-0.012707966980287463, 3.809113859067846e-15], -0.2348665397215645)
    addconstraint!(Hs, H)
    @test length(Hs) == 4

    # redundancy with almost-parallel constraints and order issues (#2169)
    Hs = [
        HalfSpace([1.3877787807814457e-17, 0.0], -7.27918967693285e-18)
        HalfSpace([3.642919299551295e-17, 2.7755575615628914e-17], 3.194076233118193e-17)
        HalfSpace([-6.938893903907228e-17, 0.03196481037863655], 0.1022675780516848)
        HalfSpace([-0.022382796523655775, -3.469446951953614e-18], 0.05268494728581621)
        HalfSpace([-8.673617379884035e-19, -0.01331598581298936], 0.020336411343083463)
        HalfSpace([0.01620615792275367, -2.949029909160572e-17], 0.03845370605895632)
       ]
    P = HPolygon(copy(Hs))
    @test ispermutation(P.constraints, Hs[3:6])
end

# default Float64 constructors
@test HPolygon() isa HPolygon{Float64, Vector{Float64}}
@test HPolygonOpt() isa HPolygonOpt{Float64, Vector{Float64}}
@test VPolygon() isa VPolygon{Float64}


text:
import LinearAlgebra.LAPACK: gels!, gesv!, getrs!, getri!

@testset "test magma_gels! $T by gels $interface" for T in [Float32, Float64, ComplexF32, ComplexF64], interface in ["CPU", "GPU"]

    # randomly generate a 2 by 2 matrix for testing
    size = 2
    A = rand(T, size, size)
    B = rand(T, size, size)

    dA = CuArray(A)
    dB = CuArray(B)

    A_copy = copy(A)
    B_copy = copy(B)

    (A_test, B_test) = interface == "GPU" ? (dA, dB) : (A_copy, B_copy)

    right_answer = gels!('N', A, B)

    magma_init()
    # result = interface == "GPU" ? magma_gels!(MAGMA.MagmaNoTrans, dA, dB) : magma_gels!(MAGMA.MagmaNoTrans, A_copy, B_copy)
    result = magma_gels!(MAGMA.MagmaNoTrans, A_test, B_test)
    magma_finalize()

    # if S is approximately equal to s, we defined then it's alright
    for i in 1:length(result)
        @test Array(result[i]) ≈ Array(right_answer[i])
    end

end

@testset "test magma_gesv! $T by gesv $interface" for T in [Float32, Float64, ComplexF32, ComplexF64], interface in ["CPU", "GPU"]

    # randomly generate a 2 by 2 matrix for testing
    size_test = 2
    A = rand(T, size_test, size_test)
    B = rand(T, size_test, size_test)

    dA = CuArray(A)
    dB = CuArray(B)

    A_copy = copy(A)
    B_copy = copy(B)

    (A_test, B_test) = interface == "GPU" ? (dA, dB) : (A_copy, B_copy)

    right_answer = gesv!(A, B)

    magma_init()
    result = magma_gesv!(A_test, B_test)
    magma_finalize()

    # if S is approximately equal to s, we defined then it's alright
    for i in 1:length(result)-1
        @test Array(result[i]) ≈ Array(right_answer[i])
    end

end

@testset "getrf/getri" begin
    @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for interface in (CuMatrix, Matrix)
            A = rand(elty,2,2)
            B = (copy(A))

            iA = inv(A)
            A, ipiv = LAPACK.getrf!(A)
            # println("A = ", A)
            # println("ipiv = ", ipiv)
            # println("inverse of A = ", iA)
            A = LAPACK.getri!(A, ipiv)

            iB = inv(B)
            B, ipivB = LAPACK.getrf!(B)

            # println("B = ", B)
            # println("ipiv = ", ipiv)
            # println("inverse of B = ", iB)
            magma_init()
            B = interface(B)
            B = magma_getri!(B, ipivB)
            magma_finalize()

            @test Array(A) ≈ Array(B)
        end
    end
end

@testset "getrs" begin
    @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
        @testset for interface in (CuMatrix, Matrix)
            A = rand(elty,2,2)
            B = rand(elty,2,2)
            trans_char = 'N'

            Atest = copy(A)
            Btest = copy(B)
            trans = MAGMA.MagmaNoTrans

            A, ipiv = LAPACK.getrf!(A)
            # println("ipiv is ", ipiv)
            # println("A = ", A)
            # println("ipiv = ", ipiv)
            # println("inverse of A = ", iA)
            B = getrs!(trans_char, A, ipiv, B)

            # println("B = ", B)
            # println("ipiv = ", ipiv)
            # println("inverse of B = ", iB)
            Atest, ipivtest = LAPACK.getrf!(Atest)
            # println("ipiv of test matrix is ", ipivtest)
            magma_init()
            Btest = magma_getrs!(trans, Atest, ipivtest, Btest)
            magma_finalize()
            # println("A ≈ Atest: ", A≈Atest)
            # println("ipiv", ipiv ≈ ipivtest)
            @test Array(B) ≈ Array(Btest)
        end
    end
end

@testset "getrf" begin
    @testset for elty in MAGMA.magmaTypeTuple
        @testset for interface in (Array, CuArray)
            # println(magmaTypeTuple)
            A = rand(elty, 2, 2)
            B = copy(A)

            magma_init()
            A, ipiv = magma_getrf!(A)
            magma_finalize()
            B, ipivB,infoB= LAPACK.getrf!(B)

            @test Array(A) ≈ Array(B)
            # @test ipiv ≈ ipivB
        end
    end
end


text:
################################
##Miscellaneous useful functions
################################

"
Returns whether to accept a Metropolis-Hastings move with acceptance probability `min(1, exp(logα))`.
"
function accMH(logα::Float64)
    logα + rand(Distributions.Exponential()) > 0.0
end

"
Returns effective sample size for weights `w`, which do not need to be normalised.
"
function ess(w::Array{Float64, 1})
    if maximum(w)==0.0
        return 0.0
    end
    sum(w)^2 / sum(w.^2)
end

"
Inserts `x` into a sorted vector `v`. Sorting order is denoted by `rev`.
"
insertsorted!(v::Vector, x; rev=false) = insert!(v, searchsortedfirst(v,x,rev=rev), x)


text:
using ModelingToolkitStandardLibrary.Electrical, ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkitStandardLibrary.Electrical: _step, _square_wave, _triangular_wave,
                                 _cos_wave, _damped_sine_wave, _ramp
# using Plots

@parameters t
@named ground = Ground()
R = 5
@named resistor = Resistor(R=R)

@info "Testing the voltage sources..."
@testset "voltage sources" begin
    @named resistor = Resistor(R=R)
    @named voltage_sensor = VoltageSensor()
    offset = 1
    freq = 50
    @named sinesource = SineVoltage(offset=offset, amplitude=1.0, frequency=freq, starttime=.5, phase=0.0)
    @named constsource = ConstantVoltage(V=1.0)
    @named stepsource = StepVoltage(height=10, offset=1.0, starttime=1.0)
    sources = [constsource, sinesource, stepsource]
    for source in sources
        rc_eqs = [
            connect(voltage_sensor.p, source.p, resistor.p)
            connect(voltage_sensor.n, source.n, resistor.n, ground.g)
             ]
        @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, voltage_sensor, ground])
        sys = structural_simplify(rc_model)
        u0 = [
            resistor.p.i => 10.0
            constsource.v => 1.0
        ]
        prob = ODEProblem(sys, u0, (0, 2.0))
        sol = solve(prob, Rosenbrock23())

        if source == constsource
            @test sol[voltage_sensor.v][1:20000:end] ≈ ones(length(sol.t[1:20000:end])) atol=1e-3
        elseif source == stepsource
            @test sol[voltage_sensor.v][1:20000:end] ≈ [(t < 1) ? offset : offset + height for t in sol.t[1:20000:end]] atol=1e-3
        else
            @test sol[voltage_sensor.v][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6
        end
    end
end

@info "Testing the sensors..."
@testset "sensors" begin
    offset = 1
    freq = 50
    @named source = SineVoltage(offset=offset, amplitude=1.0, frequency=freq, starttime=.5, phase=0.0)

    @testset "current" begin
        @named current_sensor = CurrentSensor()

        rc_eqs = [
            connect(source.p, resistor.p)
            connect(resistor.n, current_sensor.p)
            connect(current_sensor.n, source.n, ground.g)
            ]
            @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, current_sensor, ground])
            sys = structural_simplify(rc_model)
            u0 = [
                resistor.p.i => 10.0
                ]
                prob = ODEProblem(sys, u0, (0, 2.0))
                sol = solve(prob, Rosenbrock23())

                @test sol[current_sensor.i][1:20000:end] ≈ [(t < 0.5) ? offset/R : (offset + sin(2π*freq*(t .- 0.5)))/R for t in sol.t[1:20000:end]] atol=1e-6
            end

    @testset "potential" begin
        @named potential_sensor1 = PotentialSensor()
        @named potential_sensor2 = PotentialSensor()

        rc_eqs = [
                connect(source.p, resistor.p, potential_sensor1.p)
                connect(resistor.n, potential_sensor2.p, source.n, ground.g)
                ]
        @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, potential_sensor1, potential_sensor2, ground])
        sys = structural_simplify(rc_model)
        u0 = [
            resistor.p.i => 10.0
        ]
        prob = ODEProblem(sys, u0, (0, 2.0))
        sol = solve(prob, Rosenbrock23())
        @test sol[potential_sensor1.phi][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6
        @test iszero(sol[potential_sensor2.phi][1:20000:end])
    end

    @testset "voltage" begin
        @named voltage_sensor = VoltageSensor()

        rc_eqs = [
                connect(source.p, resistor.p, voltage_sensor.p)
                connect(voltage_sensor.n, source.n, resistor.n, ground.g)
                ]
        @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, voltage_sensor, ground])
        sys = structural_simplify(rc_model)
        u0 = [
            resistor.p.i => 10.0
        ]
        prob = ODEProblem(sys, u0, (0, 2.0))
        sol = solve(prob, Rosenbrock23())
        @test sol[voltage_sensor.v][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6
    end

    @testset "power" begin
        @named power_sensor = PowerSensor()

        rc_eqs = [
                connect(source.p, resistor.p, power_sensor.pv)
                connect(power_sensor.nv, resistor.n, power_sensor.nc)
                connect(power_sensor.pc, source.n, ground.g)
                ]
        @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, power_sensor, ground])
        sys = structural_simplify(rc_model)
        u0 = [
            resistor.p.i => 10.0
        ]
        prob = ODEProblem(sys, u0, (0, 2.0))
        sol = solve(prob, Rosenbrock23())
        @test sol[power_sensor.power][1:20000:end] ≈ [(t < 0.5) ? offset^2/R : (offset + sin(2π*freq*(t .- 0.5)))^2 / R for t in sol.t[1:20000:end]] atol=1e-6
    end

    @testset "multi" begin
        @named multi_sensor = MultiSensor()

        rc_eqs = [
                connect(source.p, resistor.p, multi_sensor.pv)
                connect(multi_sensor.nv, resistor.n, multi_sensor.nc)
                connect(multi_sensor.pc, source.n, ground.g)
                ]
        @named rc_model = ODESystem(rc_eqs, t, systems = [resistor, source, multi_sensor, ground])
        sys = structural_simplify(rc_model)
        u0 = [
            resistor.p.i => 10.0
        ]
        prob = ODEProblem(sys, u0, (0, 2.0))
        sol = solve(prob, Rosenbrock23())
        @test sol[multi_sensor.i][1:20000:end] ≈ [(t < 0.5) ? offset/R : (offset + sin(2π*freq*(t .- 0.5)))/R for t in sol.t[1:20000:end]] atol=1e-6
        @test sol[multi_sensor.v][1:20000:end] ≈ [(t < 0.5) ? offset : offset + sin(2π*freq*(t .- 0.5)) for t in sol.t[1:20000:end]] atol=1e-6
    end
end

@info "Testing the inductor..."
@testset "Inductor" begin
    freq, offset = 5000, 1
    @named sinesource = SineVoltage(offset=offset, starttime=0.5, amplitude=1.0, frequency=5, phase=0.0)
    @named l1 = Inductor()
    @named l2 = Inductor()
    @named voltage_sensor = VoltageSensor()
    l_eqs = [
        connect(voltage_sensor.p, sinesource.p, l1.p)
        connect(l1.n, l2.p)
        connect(voltage_sensor.n, sinesource.n, l2.n, ground.g)
            ]
    @named l_model = ODESystem(l_eqs, t, systems = [l1, l2, sinesource, voltage_sensor, ground])
    sys = structural_simplify(l_model)
    equations(sys)
    u0 = [
        l1.p.i => 10.0
        sinesource.v => 1.0
        l2.v => 0.0
    ]
    prob = ODEProblem(sys, u0, (0, 10.0))
    sol = solve(prob, Rosenbrock23())

    @test sol[l1.v] + sol[l2.v] ≈ sol[voltage_sensor.v]
end

@info "Constructing an integrator circuit..."
# This tests Capacitor, IdealOpAmp
@testset "Integrator" begin
    @named ground = Ground()
    @named res1 = Resistor()
    @named c1 = Capacitor()
    @named opamp = IdealOpAmp()
    @named square = SquareVoltage()

    in_eqs = [
        connect(square.p, res1.p)
        connect(res1.n, c1.p, opamp.p1)
        connect(opamp.n2, c1.n)
        connect(opamp.n1, ground.g, opamp.p2, square.n)
    ]
    @named in_model = ODESystem(in_eqs, t, systems = [res1,  opamp, square, c1, ground])
    sys = structural_simplify(in_model)
    u0 = [
        c1.v => 1
        res1.v => 1
    ]
    prob = ODEProblem(sys, u0, (0, 10.0))
    sol = solve(prob, Rosenbrock23())
    @test sol[opamp.v2] == sol[c1.v] # Not a great one however. Rely on the plot

    # plt = plot(sol)
    # savefig(plt, "integrator")
end

@info "Testing voltage function generators..."
@testset "Voltage function generators" begin
    st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001

    @named ground = Ground()
    @named res = Resistor()
    @named voltage_sensor = VoltageSensor()
    @named vstep = StepVoltage(starttime=st, offset=o, height=h)
    @named vsquare = SquareVoltage(offset=o, starttime=st, amplitude=A, frequency=f)
    @named vtri = TriangularVoltage(offset=o, starttime=st, amplitude=A, frequency=f)
    # @named vsawtooth = SawToothVoltage(amplitude=A, starttime=st, frequency=f, offset=o)
    @named vcosine = CosineVoltage(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ)
    @named vdamped_sine = DampedSineVoltage(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ, damping_coef=d)
    @named vramp = RampVoltage(offset=o, starttime=st, endtime=et, height=h)

    vsources = [vtri, vsquare, vstep, vcosine, vdamped_sine, vramp]
    waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _square_wave.(x, δ, f, A, st),
                                o .+  _step.(x, δ, h, st),
                                # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ),
                                o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d),
                                o .+ _ramp.(x, δ, st, et, h)], i)
    for i in 1:length(vsources)
        vsource = vsources[i]
        @info Symbolics.getname(vsource)
        eqs = [
            connect(vsource.p, voltage_sensor.p, res.p)
            connect(ground.g, voltage_sensor.n, vsource.n, res.n)
        ]
        @named vmodel = ODESystem(eqs, t, systems = [voltage_sensor, res, vsource, ground])
        vsys = structural_simplify(vmodel)

        u0 = [
            vsource.v => 1
            res.v => 1
        ]

        prob = ODEProblem(vsys, u0, (0, 10.0))
        vsol = solve(prob, dt=0.1, Rosenbrock23())

        @test vsol[vsource.v][1150:end] ≈ waveforms(i, vsol.t)[1150:end] atol=1e-1
        # For visual inspection
        # plt = plot(sol)
        # savefig(plt, "test_current_$(Symbolics.getname(source))")
    end
end

@info "Testing the Current generators..."
@testset "Current function generators" begin
    st, o, h, f, A, et, ϕ, d, δ = 0.7, 1.25, 3, 2, 2.5, 2.5, π/4, 0.1, 0.0001

    @named ground = Ground()
    @named res = Resistor()
    @named current_sensor = CurrentSensor()
    @named istep = StepCurrent(starttime=st, offset=o, height=h)
    @named isquare = SquareCurrent(offset=o, starttime=st, amplitude=A, frequency=f)
    @named itri = TriangularCurrent(offset=o, starttime=st, amplitude=A, frequency=f)
    # @named isawtooth = SawToothCurrent(amplitude=A, starttime=st, frequency=f, offset=o)
    @named icosine = CosineCurrent(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ)
    @named idamped_sine = DampedSineCurrent(offset=o, amplitude=A, frequency=f, starttime=st, phase=ϕ, damping_coef=d)
    @named iramp = RampCurrent(offset=o, starttime=st, endtime=et, height=h)

    isources = [itri, isquare, istep, icosine, idamped_sine, iramp]
    waveforms(i, x) = getindex([o .+ (x .> st) .* _triangular_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _square_wave.(x, δ, f, A, st),
                                o .+  _step.(x, δ, h, st),
                                # o .+ (x .> st). * _sawtooth_wave.(x, δ, f, A, st),
                                o .+ (x .> st) .* _cos_wave.(x, f, A, st, ϕ),
                                o .+ (x .> st) .* _damped_sine_wave.(x, f, A, st, ϕ, d),
                                o .+ _ramp.(x, δ, st, et, h)], i)

    for i in 1:length(isources)
        isource = isources[i]
        eqs = [
            connect(isource.p, current_sensor.n)
            connect(current_sensor.p, res.n)
            connect(isource.n, res.p)
        ]
        @named model = ODESystem(eqs, t, systems = [current_sensor, isource, res])
        isys = alias_elimination(model)

        u0 = [
            isource.i => 1
            res.v => 1
        ]
        prob = ODEProblem(isys, u0, (0, 2.0))
        sol = solve(prob, Rosenbrock23())

        @test sol[isource.i][1150:end] ≈ waveforms(i, sol.t)[1150:end] atol=1e-1
        # For visual inspection
        # plt = plot(sol)
        # savefig(plt, "test_current_$(Symbolics.getname(source))")
    end
end


text:
using ModernGL

"""
	mutable struct GLVertexBuffer

Submitting vertex data for rendering requires creating a stream of vertices, and then telling OpenGL how to interpret that stream. `GLVertexBuffer` provides methods for uploading vertex data (position, normal vector, color, etc.) to the device for non-immediate-mode rendering.
A Vertex Buffer Object is the common term for a Buffer Object when it is used as a source for vertex array data.
Two inner constructor methods, with and without a `Vector{Float32}`.
"""
mutable struct GLVertexBuffer

	id::Int32
	vector::Vector{Float32}

	# constructor
	function GLVertexBuffer()
		ret=new(-1,[])
		finalizer(releaseGpuResources, ret)
		return ret
	end

	# constructor
	function GLVertexBuffer(vector::Vector{Float32})
		ret=new(-1,vector)
		finalizer(releaseGpuResources, ret)
		return ret
	end

end

"""

# Example

Release an array of `GLVertexBuffer`, a Gpu resource. It is a list of indices that will select which vertices to use and in which order.
"""
function releaseGpuResources(buffer::GLVertexBuffer)
	global __release_gpu_resources__
	if buffer.id>=0
		id=buffer.id
		buffer.id=-1
		glDeleteLater(function()  glDeleteBuffers(1,[id]) end)
	end
end

"""
	enableAttribute(location::Int32,buffer::GLVertexBuffer,num_components::Int64)

If enabled, the values in the generic vertex attribute array will be accessed and used for rendering when calls are made to vertex array commands.
"""
function enableAttribute(location::Int32,buffer::GLVertexBuffer,num_components::Int64)
	if length(buffer.vector)==00 || location<0 return end
	if buffer.id<0 buffer.id=glGenBuffer() end
	# Enable blending
	glEnable(GL_BLEND)
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA)

	glLineWidth()
	glPointSize()
	glBindBuffer(GL_ARRAY_BUFFER, buffer.id)
	glBufferData(GL_ARRAY_BUFFER, sizeof(buffer.vector), buffer.vector, GL_STATIC_DRAW)
	glVertexAttribPointer(location,num_components,GL_FLOAT,false,0,C_NULL)
	glEnableVertexAttribArray(location)
	glBindBuffer(GL_ARRAY_BUFFER, 0)
end

"""
	disableAttribute(location::Int32,buffer::GLVertexBuffer)

By default, all client-side capabilities are disabled, including all generic vertex attribute arrays.
"""
function disableAttribute(location::Int32,buffer::GLVertexBuffer)
	if length(buffer.vector)==00 || location<0 return end
	glDisableVertexAttribArray(location)
end