text:
using Random
using DelimitedFiles
using Statistics

# Tested under Julia 1.6.2

distance(emb, i, j) =
    sqrt(sum((emb[i+1, 2:end] - emb[j+1, 2:end]).^ 2))

# the code assumes
# - nodes are labelled from 0
# - embedding is sorted on node id; node id can be ignored and is a first column
# - files are space delimited
function reshuffle_edges(edge_file::AbstractString,
                         community_file::AbstractString,
                         embedding_file::AbstractString,
                         out_edge_file::AbstractString,
                         shufflepct::Real;
                         seed::Int=1234,
                         check_correctness::Bool=true)
    edges_raw = readdlm(edge_file, Int)
    communities = readdlm(community_file, Int)
    embedding = readdlm(embedding_file, Float64)

    @assert size(edges_raw, 2) == 2
    vmin, vmax = extrema(edges_raw)
    @assert vmin >= 0
    @assert vmax < length(communities)
    @assert size(communities, 2) == 1
    @assert size(embedding,1) == length(communities)

    comm_edges = Dict{Int, Set{Tuple{Int,Int}}}()
    background = Set{Tuple{Int,Int}}()

    for (a, b) in eachrow(edges_raw)
        com = communities[a+1]
        if com == communities[b+1]
            comm_set = get!(comm_edges, com, Set{Tuple{Int, Int}}())
            @assert !((a, b) in comm_set)
            push!(comm_set, (a, b))
        else
            @assert !((a, b) in background)
            push!(background, (a, b))
        end
    end
    @assert sum(length.(values(comm_edges))) + length(background) == size(edges_raw, 1)

    Random.seed!(seed)
    for comm_set in values(comm_edges)
        if check_correctness
            avg_dist = mean(distance(embedding, a, b) for (a,b) in comm_set)
            @info "old avg distance: $avg_dist"
        end
        len = length(comm_set)
        must_rewire = len * shufflepct
        did_rewire = 0
        while did_rewire < must_rewire
            (a1, a2) = rand(comm_set)
            (b1, b2) = rand(comm_set)
            (a1, a2) == (b1, b2) && continue
            (a1, b2) == (b1, a2) && continue
            (a1, b2) in comm_set && continue
            (b1, a2) in comm_set && continue
            @assert a1 >= 0
            @assert a2 >= 0
            @assert b1 >= 0
            @assert b2 >= 0
            old_sum = distance(embedding, a1, a2) + distance(embedding, b1, b2)
            new_sum = distance(embedding, a1, b2) + distance(embedding, b1, a2)
            old_sum >= new_sum && continue
            pop!(comm_set, (a1, a2))
            pop!(comm_set, (b1, b2))
            push!(comm_set, (a1, b2))
            push!(comm_set, (b1, a2))
            did_rewire += 1
            @assert length(comm_set) == len
        end
        if check_correctness
            avg_dist = mean(distance(embedding, a, b) for (a,b) in comm_set)
            @info "new avg distance: $avg_dist"
            println()
        end
    end

    @assert sum(length.(values(comm_edges))) + length(background) == size(edges_raw, 1)
    open(out_edge_file, "w") do io
        for (a, b) in background
            println(io, a, " ", b)
        end
        for comm_set in values(comm_edges)
            for (a, b) in comm_set
                println(io, a, " ", b)
            end
        end
    end

    check_correctness || return

    edges_raw2 = readdlm(out_edge_file, Int)
    if iszero(shufflepct)
        @assert sort!(collect(eachrow(edges_raw))) ==
                sort!(collect(eachrow(edges_raw2)))
    else
        @assert sort!(edges_raw[:, 1]) == sort!(edges_raw2[:, 1])
        @assert sort!(edges_raw[:, 2]) == sort!(edges_raw2[:, 2])
    end
end

reshuffle_edges(ARGS[1],
                ARGS[2],
                ARGS[3],
                ARGS[4], parse(Float64,ARGS[5]))


text:
# included in runtests.jl
# The first test set is an close adaption of /stdlib/Logging/test/runtests.jl
# The second test set checks additional functionality
using Test
using WebSockets
import WebSockets: Logging
import Logging: min_enabled_level, shouldlog, handle_message, with_logger
import Base.CoreLogging: Info,
                         Debug,
                         Warn,
                         Error,
                         LogLevel

@noinline func1() = backtrace()

@testset "WebSocketLogger" begin
    # First pass log limiting
    @test min_enabled_level(WebSocketLogger(devnull, Debug)) == Debug
    @test min_enabled_level(WebSocketLogger(devnull, Error)) == Error

    # Second pass log limiting
    logger = WebSocketLogger(devnull)
    @test shouldlog(logger, Info, Base, :group, :asdf) === true
    handle_message(logger, Info, "msg", Base, :group, :asdf, "somefile", 1, maxlog=2)
    @test shouldlog(logger, Info, Base, :group, :asdf) === true
    handle_message(logger, Info, "msg", Base, :group, :asdf, "somefile", 1, maxlog=2)
    @test shouldlog(logger, Info, Base, :group, :asdf) === false

    # Check that maxlog works without an explicit ID (#28786)
    buf = IOBuffer()
    io = IOContext(buf, :displaysize=>(30,80), :color=>false)
    logger = WebSocketLogger(io)
    # This covers an issue (#28786) in stdlib/Julia v 0.7, where log_record_id is not defined
    # in the expanded macro. Fixed Nov 2018, labelled for backporting to v1.0.0
    if VERSION >= v"1.0.0"
        with_logger(logger) do
            for i in 1:2
                @info "test" maxlog=1
            end
        end
        @test String(take!(buf)) ==
        """
        [ Info: test
        """
        with_logger(logger) do
            for i in 1:2
                @info "test" maxlog=0
            end
        end
        @test String(take!(buf)) == ""
    end
    @testset "Default metadata formatting" begin
        @test Logging.default_metafmt(Debug, Base, :g, :i, expanduser("~/somefile.jl"), 42) ==
            (:blue,      "Debug:",   "@ Base ~/somefile.jl:42")
        @test Logging.default_metafmt(Info,  Main, :g, :i, "a.jl", 1) ==
            (:cyan,      "Info:",    "")
        @test Logging.default_metafmt(Warn,  Main, :g, :i, "b.jl", 2) ==
            (:yellow,    "Warning:", "@ Main b.jl:2")
        @test Logging.default_metafmt(Error, Main, :g, :i, "", 0) ==
            (:light_red, "Error:",   "@ Main :0")
        # formatting of nothing
        @test Logging.default_metafmt(Warn,  nothing, :g, :i, "b.jl", 2) ==
            (:yellow,    "Warning:", "@ b.jl:2")
        @test Logging.default_metafmt(Warn,  Main, :g, :i, nothing, 2) ==
            (:yellow,    "Warning:", "@ Main")
        @test Logging.default_metafmt(Warn,  Main, :g, :i, "b.jl", nothing) ==
            (:yellow,    "Warning:", "@ Main b.jl")
        @test Logging.default_metafmt(Warn,  nothing, :g, :i, nothing, 2) ==
            (:yellow,    "Warning:", "")
        @test Logging.default_metafmt(Warn,  Main, :g, :i, "b.jl", 2:5) ==
            (:yellow,    "Warning:", "@ Main b.jl:2-5")
    end

    function dummy_metafmt(level, _module, group, id, file, line)
        :cyan,"PREFIX","SUFFIX"
    end

    # Log formatting
    function genmsg(message; level=Info, _module=Main,
                    file="some/path.jl", line=101, color=false, width=75,
                    meta_formatter=dummy_metafmt, show_limited=true,
                    right_justify=0, kws...)
        buf = IOBuffer()
        io = IOContext(buf, :displaysize=>(30,width), :color=>color)
        logger = WebSocketLogger(io, Debug,
                               meta_formatter=meta_formatter,
                               show_limited=show_limited,
                               right_justify=right_justify)
        Logging.handle_message(logger, level, message, _module, :a_group, :an_id,
                               file, line; kws...)
        String(take!(buf))
    end

    # Basic tests for the default setup
    @test genmsg("msg", level=Info, meta_formatter=Logging.default_metafmt) ==
    """
    [ Info: msg
    """
    @test genmsg("line1\nline2", level=Warn, _module=Base,
                 file="other.jl", line=42, meta_formatter=Logging.default_metafmt) ==
    """
    ┌ Warning: line1
    │ line2
    └ @ Base other.jl:42
    """
    # Full metadata formatting
    @test genmsg("msg", level=Debug,
                 meta_formatter=(level, _module, group, id, file, line)->
                                (:white,"Foo!", "$level $_module $group $id $file $line")) ==
    """
    ┌ Foo! msg
    └ Debug Main a_group an_id some/path.jl 101
    """

    @testset "Prefix and suffix layout" begin
        @test genmsg("") ==
        replace("""
        ┌ PREFIX EOL
        └ SUFFIX
        """, "EOL"=>"")
        @test genmsg("msg") ==
        """
        ┌ PREFIX msg
        └ SUFFIX
        """
        # Behavior with empty prefix / suffix
        @test genmsg("msg", meta_formatter=(args...)->(:white, "PREFIX", "")) ==
        """
        [ PREFIX msg
        """
        @test genmsg("msg", meta_formatter=(args...)->(:white, "", "SUFFIX")) ==
        """
        ┌ msg
        └ SUFFIX
        """
    end

    @testset "Metadata suffix, right justification" begin
        @test genmsg("xxx", width=20, right_justify=200) ==
        """
        [ PREFIX xxx  SUFFIX
        """
        @test genmsg("xxx\nxxx", width=20, right_justify=200) ==
        """
        ┌ PREFIX xxx
        └ xxx         SUFFIX
        """
        # When adding the suffix would overflow the display width, add it on
        # the next line:
        @test genmsg("xxxx", width=20, right_justify=200) ==
        """
        ┌ PREFIX xxxx
        └             SUFFIX
        """
        # Same for multiline messages
        @test genmsg("""xxx
                        xxxxxxxxxx""", width=20, right_justify=200) ==
        """
        ┌ PREFIX xxx
        └ xxxxxxxxxx  SUFFIX
        """
        @test genmsg("""xxx
                        xxxxxxxxxxx""", width=20, right_justify=200) ==
        """
        ┌ PREFIX xxx
        │ xxxxxxxxxxx
        └             SUFFIX
        """
        # min(right_justify,width) is used
        @test genmsg("xxx", width=200, right_justify=20) ==
        """
        [ PREFIX xxx  SUFFIX
        """
        @test genmsg("xxxx", width=200, right_justify=20) ==
        """
        ┌ PREFIX xxxx
        └             SUFFIX
        """
    end

    # Keywords
    @test genmsg("msg", a=1, b="asdf") ==
    """
    ┌ PREFIX msg
    │   a = 1
    │   b = "asdf"
    └ SUFFIX
    """
    # Exceptions shown with showerror
    @test genmsg("msg", exception=DivideError()) ==
    """
    ┌ PREFIX msg
    │   exception = DivideError: integer division error
    └ SUFFIX
    """

    # Attaching backtraces
    bt = func1()
    @test startswith(genmsg("msg", exception=(DivideError(),bt)),
    """
    ┌ PREFIX msg
    │   exception =
    │    DivideError: integer division error
    │    Stacktrace:
    │     $(VERSION < v"1.6" ? "" : " ")[1] func1()""")

    @testset "Limiting large data structures" begin
        @test genmsg("msg", a=fill(1.00001, 100,100), b=fill(2.00002, 10,10)) ==
        """
        ┌ PREFIX msg
        │   a =
        │    100×100 $(Matrix{Float64}):
        │     1.00001  1.00001  1.00001  1.00001  …  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001     1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001     1.00001  1.00001  1.00001
        │     ⋮                                   ⋱                    $(VERSION < v"1.1" ? "       " : "")
        │     1.00001  1.00001  1.00001  1.00001     1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001     1.00001  1.00001  1.00001
        │   b =
        │    10×10 $(Matrix{Float64}):
        │     2.00002  2.00002  2.00002  2.00002  …  2.00002  2.00002  2.00002
        │     2.00002  2.00002  2.00002  2.00002     2.00002  2.00002  2.00002
        │     2.00002  2.00002  2.00002  2.00002     2.00002  2.00002  2.00002
        │     ⋮                                   ⋱                    $(VERSION < v"1.1" ? "       " : "")
        │     2.00002  2.00002  2.00002  2.00002     2.00002  2.00002  2.00002
        │     2.00002  2.00002  2.00002  2.00002     2.00002  2.00002  2.00002
        └ SUFFIX
        """
        # Limiting the amount which is printed
        @test genmsg("msg", a=fill(1.00001, 10,10), show_limited=false) ==
        """
        ┌ PREFIX msg
        │   a =
        │    10×10 $(Matrix{Float64}):
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        │     1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001  1.00001
        └ SUFFIX
        """
    end

    # Basic colorization test
    @test genmsg("line1\nline2", color=true) ==
    """
    \e[36m\e[1m┌ \e[22m\e[39m\e[36m\e[1mPREFIX \e[22m\e[39mline1
    \e[36m\e[1m│ \e[22m\e[39mline2
    \e[36m\e[1m└ \e[22m\e[39m\e[90mSUFFIX\e[39m
    """

end

@testset "WebSocketLogger_special" begin
    # Log formatting
    function spec_msg(message; level=Info,
                    color = false,
                    shouldlog = WebSockets.shouldlog_default)
        buf = IOBuffer()
        io = IOContext(buf, :color => color)
        logger = WebSocketLogger(io, Debug,
                               shouldlog = shouldlog)
        with_logger(logger) do
                if level == Info
                    @info message
                elseif level == Wslog
                    @wslog message
                elseif level == Warn
                    @warn message
                elseif level == Debug
                    @debug message
                else
                    @error message
                end
            end
        String(take!(buf))
    end

    # Normal info
    @test spec_msg("---") == "[ Info: ---\n"

    # Debug, own format
    msg = spec_msg("---", level = Debug)
    @test startswith(msg, "┌ Debug: ---\n└ ")
    @test '@' in msg

    # Loglevel Wslog, own format
    # This covers an issue related to #28786 in stdlib/Julia v 0.7, where log_record_id is not defined
    # in the expanded macro. Why it is triggered here, but not in the LogLevel Info test is unknown.
    if VERSION >= v"1.0.0"
        msg = spec_msg("---", level = Wslog)
        @test startswith(msg, "[ Wslog ")
        @test endswith(msg, ": ---\n")
    end
    # Blocking all log levels except Wslog.
    spec_shouldlog(logger, level, _module, group, id) = level == Wslog
    # Issue 28786
    if VERSION >= v"1.0.0"
        msg = spec_msg("---", level = Wslog, shouldlog = spec_shouldlog)
        @test startswith(msg, "[ Wslog ")
        @test endswith(msg, ": ---\n")
    end
    @test spec_msg("---", shouldlog = spec_shouldlog) == ""
    @test spec_msg("---", level = Debug, shouldlog = spec_shouldlog) == ""
    @test spec_msg("---", level = Warn, shouldlog = spec_shouldlog) == ""
    @test spec_msg("---", level = Error, shouldlog = spec_shouldlog) == ""

    # Check two lines output for all except log level Info and Wslog
    @test startswith(spec_msg("---", level = Debug), "┌ Debug: ---\n└ ")
    @test startswith(spec_msg("---", level = Warn), "┌ Warn: ---\n└ ")
    @test startswith(spec_msg("---", level = Error), "┌ Error: ---\n└ ")

    # Check colors in first argument show methods are preserved
    ws = WebSocket(IOBuffer(), true)
    @test spec_msg(ws) == "[ Info: WebSocket{GenericIOBuffer}(server, CONNECTED)\n"
    @test spec_msg(ws, color= true) == "\e[36m\e[1m[ \e[22m\e[39m\e[36m\e[1mInfo: \e[22m\e[39mWebSocket{GenericIOBuffer}(server, \e[32mCONNECTED\e[39m)\n"

    # Check tuples in first argument, resembling println(arg1, arg2)
    @test spec_msg(("Now testing ", ws), color=true) == "\e[36m\e[1m[ \e[22m\e[39m\e[36m\e[1mInfo: \e[22m\e[39mNow testing WebSocket{GenericIOBuffer}(server, \e[32mCONNECTED\e[39m)\n"

    # Check tuples in first argument, resembling println(arg1, arg2)
    @test spec_msg(("Now testing ", ws), color=true) == "\e[36m\e[1m[ \e[22m\e[39m\e[36m\e[1mInfo: \e[22m\e[39mNow testing WebSocket{GenericIOBuffer}(server, \e[32mCONNECTED\e[39m)\n"

    # Check default shouldlog filters away WebSockets. HTTP.Servers, i.e. HTTP.Servers
    logger = WebSocketLogger()
    @test shouldlog(logger, Info, WebSockets.HTTP.Servers, :group, :asdf) == false
    @test shouldlog(logger, Info, Main, :group, :asdf) == true

    # Check that error handling works with @wslog
    if VERSION >= v"1.0.0"
        # This covers an issue (#28786) in stdlib/Julia v 0.7, where log_record_id is not defined
        # in the expanded macro. Fixed Nov 2018, labelled for backporting to v1.0.0
        buf = IOBuffer()
        io = IOContext(buf, :displaysize=>(30,80), :color=>true)
        logger = WebSocketLogger(io)
        with_logger(logger) do
               @wslog sqrt(-2)
               end
       @test length(String(take!(buf))) > 1900
       """
       [ Info: test
       """
   end
end
nothing

ext: jl


text:
# Get status of the registration
function status(r::HTTP.Request)
    r.method == "GET" || return json(405; error="Method not allowed")

    # Only allow authenticated users
    state = getcookie(r, "state")
    if isempty(state) || !haskey(USERS, state)
        return json(400; error="Missing or invalid state cookie")
    end

    id = getquery(r, "id")
    haskey(REGISTRATIONS, id) || return json(404; state="unknown")

    # 10 second arbitrary timeout.
    for _ in 1:10
        regstate = REGISTRATIONS[id]
        if regstate.state !== :pending
            delete!(REGISTRATIONS, id)
            return json(; state=regstate.state, message=regstate.msg)
        end
        sleep(1)
    end

    return json(; state="pending", message="Registration in progress...")
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct AdminEmojiApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminEmojiAdd(_api::AdminEmojiApi, name::String, url::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.emoji.add", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_param(_ctx.form, "url", url)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Add an emoji.
Param: name::String (required)
Param: url::String (required)
Return: Dict{String, Any}
"""
function adminEmojiAdd(_api::AdminEmojiApi, name::String, url::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiAdd(_api, name, url; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminEmojiAdd(_api::AdminEmojiApi, response_stream::Channel, name::String, url::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiAdd(_api, name, url; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminEmojiAddAlias(_api::AdminEmojiApi, name::String, alias_for::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.emoji.addAlias", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_param(_ctx.form, "alias_for", alias_for)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Add an emoji alias.
Param: name::String (required)
Param: alias_for::String (required)
Return: Dict{String, Any}
"""
function adminEmojiAddAlias(_api::AdminEmojiApi, name::String, alias_for::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiAddAlias(_api, name, alias_for; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminEmojiAddAlias(_api::AdminEmojiApi, response_stream::Channel, name::String, alias_for::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiAddAlias(_api, name, alias_for; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminEmojiList(_api::AdminEmojiApi; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.emoji.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

List emoji for an Enterprise Grid organization.
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function adminEmojiList(_api::AdminEmojiApi; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiList(_api; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminEmojiList(_api::AdminEmojiApi, response_stream::Channel; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiList(_api; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminEmojiRemove(_api::AdminEmojiApi, name::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.emoji.remove", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Remove an emoji across an Enterprise Grid organization
Param: name::String (required)
Return: Dict{String, Any}
"""
function adminEmojiRemove(_api::AdminEmojiApi, name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiRemove(_api, name; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminEmojiRemove(_api::AdminEmojiApi, response_stream::Channel, name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiRemove(_api, name; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminEmojiRename(_api::AdminEmojiApi, name::String, new_name::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.emoji.rename", ["slackAuth"])
    Swagger.set_param(_ctx.form, "name", name)  # type String
    Swagger.set_param(_ctx.form, "new_name", new_name)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Rename an emoji.
Param: name::String (required)
Param: new_name::String (required)
Return: Dict{String, Any}
"""
function adminEmojiRename(_api::AdminEmojiApi, name::String, new_name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiRename(_api, name, new_name; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminEmojiRename(_api::AdminEmojiApi, response_stream::Channel, name::String, new_name::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminEmojiRename(_api, name, new_name; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminEmojiAdd, adminEmojiAddAlias, adminEmojiList, adminEmojiRemove, adminEmojiRename


text:
#!/usr/bin/env julia

# Given some model parameters, plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "--config"
    help = "a YAML configuration file"
    default = "config.yaml"
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using DiskJockey.constants
using DiskJockey.image
using DiskJockey.model
using HDF5

import PyPlot.plt
using LaTeXStrings
import Images

# used to plot masked values
using PyCall
@pyimport numpy.ma as ma
import NaNMath

species = config["species"]
transition = config["transition"]
lam0 = lam0s[species * transition]
model = config["model"]
params = config["parameters"]

# Choose a diverging colorscheme for in front of/ behind the plane
cmap = plt[:get_cmap]("coolwarm")
cmap_img = plt[:get_cmap]("Blues")

# Use the same configuration as plot_chmaps, but we'll use different scaling.
function plot_tausurf(img::image.TausurfImg; fname = "tausurf.png")

    plotcbar = false

    data = img.data ./ AU # Convert from cm to AU

    vmin = abs(NaNMath.minimum(data)) # [AU]
    vmax = NaNMath.maximum(data) # [AU]
    vvmax = maximum([vmin, vmax])

    norm = PyPlot.matplotlib[:colors][:Normalize](-vvmax, vvmax)
    # norm = PyPlot.matplotlib[:colors][:SymLogNorm](linthresh=1,vmin=-10, vmax=10)

    (im_ny, im_nx) = size(img.data)[1:2] # y and x dimensions of the image

    xxx = ((Float64[i for i = 0:im_nx - 1] + 0.5) - im_nx / 2.) * img.pixsize_x ./ AU
    yyy = ((Float64[i for i = 0:im_ny - 1] + 0.5) - im_ny / 2.) * img.pixsize_y ./ AU

    ext = (xxx[1], xxx[end], yyy[1], yyy[end])

    # Figure out how many plots we'll have.
    ncols = 8
    nrows = ceil(Int, nlam / ncols)

    xx = 1.5 * 9
    dx = 1.5
    yy = (nrows + 1) * 1.5
    dy = 1.5

    fig, ax = plt[:subplots](nrows = nrows, ncols = ncols, figsize = (xx, yy))

    for row = 1:nrows
        for col = 1:ncols
            iframe = col + (row - 1) * ncols

            if col != 1 || row != nrows
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$x$ (AU)")
                ax[row, col][:set_ylabel](L"$y$ (AU)")
            end

            if iframe > nlam
                continue # just pass through for empty channels

            else
                # No flip is needed for the tausurf because it was never put into a SkyImage. It can be plotted as is.
                frame = data[:,:,iframe]

                # Stuff the actual frame into a numpy MaskedArray object, where the masked values are the NaNs.
                mframe = pycall(ma.array, Any, frame, mask = isnan.(frame))

                # If the entire image is NaN, then just skip plotting it.
                if all(isnan, frame)
                    continue
                else
                    im = ax[row, col][:imshow](mframe, extent = ext, interpolation = "none", origin = "lower", cmap = cmap, norm = norm)

                    if !plotcbar
                        # Plot the colorbar
                        cax = fig[:add_axes]([(xx - 0.35 * dx) / xx, (yy - 1.5 * dy) / yy, (0.1 * dx) / xx, dy / yy])
                        cbar = fig[:colorbar](mappable = im, cax = cax)

                        cbar[:ax][:tick_params](labelsize = 6)
                        fig[:text](0.99, (yy - 1.7 * dy) / yy, "AU", size = 8, ha = "right")
                        plotcbar = true
                    end

                end

                ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords = "axes fraction", size = 8)
            end

        end
    end

    fig[:subplots_adjust](hspace = 0.00, wspace = 0.00, top = (yy - 0.5 * dy) / yy, bottom = (0.5 * dy) / yy, left = (0.5 * dx) / xx, right = (xx - 0.5 * dy) / xx)

    plt[:savefig](fname, dpi = 600)
end

# Use the same configuration as plot_chmaps, but we'll use different scaling.
function plot_taupos(pos::image.TausurfPos; fname = "taupos.png")

    data_x = pos.data_x ./ AU # Convert from cm to AU
    data_y = pos.data_y ./ AU # Convert from cm to AU
    data_z = pos.data_z ./ AU # Convert from cm to AU

    # Figure out how many plots we'll have.
    ncols = 8
    nrows = ceil(Int, nlam / ncols)

    xx = 1.5 * 9
    dx = 1.5
    yy = (nrows + 1) * 1.5
    dy = 1.5

    fig, ax = plt[:subplots](nrows = nrows, ncols = ncols, figsize = (xx, yy), sharex = true, sharey = true)

    for row = 1:nrows
        for col = 1:ncols
            iframe = col + (row - 1) * ncols

            if col != 1 || row != nrows
                ax[row, col][:xaxis][:set_ticklabels]([])
                ax[row, col][:yaxis][:set_ticklabels]([])
            else
                ax[row, col][:set_xlabel](L"$xx$ (AU)")
                ax[row, col][:set_ylabel](L"$yy$ (AU)")
            end

            if iframe > nlam
                continue # skip the plot

            else

                ax[row, col][:plot](data_x[:,:,iframe], data_y[:,:,iframe], "k.", ms = 2)
                ax[row, col][:annotate](@sprintf("%.1f", vels[iframe]), (0.1, 0.8), xycoords = "axes fraction", size = 8)
            end

        end
    end

    fig[:subplots_adjust](hspace = 0.00, wspace = 0.00, top = (yy - 0.5 * dy) / yy, bottom = (0.5 * dy) / yy, left = (0.5 * dx) / xx, right = (xx - 0.5 * dy) / xx)

    plt[:savefig](fname, dpi = 600)
end

"Plot the (x,y), (y,z), and (x,z) views through the disk, along with the (projected) direction of the observer. Assumes that P.A. of the synthesized disk is 0."
function plot_taupos_tryptych(pos::image.TausurfPos, img::image.SkyImage, incl::Real; fname = "taupos_tryptych.png")
    # Figure out how many plots we'll have.

    data_x = pos.data_x ./ AU # Convert from cm to AU
    data_y = pos.data_y ./ AU # Convert from cm to AU
    data_z = pos.data_z ./ AU # Convert from cm to AU

    ncols = 4
    nrows = nlam

    vvmax = maximum(abs, img.data)
    norm = PyPlot.matplotlib[:colors][:Normalize](0, vvmax)

    # Image needs to be flipped along RA dimension
    ext = (img.ra[end], img.ra[1], img.dec[1], img.dec[end])

    dx = 2.5
    xx = dx * (ncols + 1)
    dy = 2.5
    yy = (nrows + 1) * dy

    fig, ax = plt[:subplots](nrows = nrows, ncols = ncols, figsize = (xx, yy))

    for row = 1:nrows

        ax[row,1][:set_xlabel](L"$\Delta \alpha$ ('')")
        ax[row,1][:set_ylabel](L"$\Delta \delta$ ('')")
        ax[row,2][:set_xlabel](L"$x$ (AU)")
        ax[row,2][:set_ylabel](L"$y$ (AU)")
        ax[row,3][:set_xlabel](L"$y$ (AU)")
        ax[row,3][:set_ylabel](L"$z$ (AU)")
        ax[row,4][:set_xlabel](L"$x$ (AU)")
        ax[row,4][:set_ylabel](L"$z$ (AU)")

        # Plot the SkyImage
        frame = flipdim(img.data[:,:,row], 2)
        im = ax[row, 1][:imshow](frame, extent = ext, interpolation = "none", origin = "lower", cmap = cmap_img, norm = norm)

        # Plot and set the bounds to have equal aspect area.
        ax[row, 2][:plot](data_x[:,:,row], data_y[:,:,row], "k.", ms = 2)
        ax[row, 2][:set_aspect]("equal", "datalim")

        ax[row, 3][:plot](data_y[:,:,row], data_z[:,:,row], "k.", ms = 2)
        ax[row, 3][:set_aspect]("equal", "datalim")

        ax[row, 4][:plot](data_x[:,:,row], data_z[:,:,row], "k.", ms = 2)
        ax[row, 4][:set_aspect]("equal", "datalim")

        ax[row, 1][:annotate](@sprintf("%.1f", vels[row]), (0.1, 0.8), xycoords = "axes fraction", size = 8)
    end

    fig[:subplots_adjust](hspace = 0.35, wspace = 0.35, top = (yy - 0.5 * dy) / yy, bottom = (0.5 * dy) / yy, left = (0.5 * dx) / xx, right = (xx - 0.5 * dy) / xx)

    plt[:savefig](fname, dpi = 120)
end

# Read the image
im = taureadImg()
# Do the velocity conversion here for plot labels
global nlam = length(im.lams)
# convert wavelengths to velocities
global vels = c_kms * (im.lams .- lam0) / lam0

pos = taureadPos()

# Make sure we are always plotting channels from blueshift to redshift.
if vels[2] < vels[1]
  println("Flipping channel maps to plot in order of increasing velocity.")
  im.data = im.data[:,:,end:-1:1]

  pos.data_x = pos.data_x[:,:,end:-1:1]
  pos.data_y = pos.data_y[:,:,end:-1:1]
  pos.data_z = pos.data_z[:,:,end:-1:1]

  vels = vels[end:-1:1]
end

pars = convert_dict(config["parameters"], config["model"])

# Assume that image has already been generated by synthesize_model.jl
img = imread()
skimg = imToSky(img, pars.dpc)

# plot_tausurf(im)
plot_taupos_tryptych(pos, skimg, params["incl"])


text:
@testset "CANDECOMP" begin
    Random.seed!(12345)

    r = 2 # core tensor dimensions
    T = rand_kruskal3(r, (10, 20, 30), false)
    initial_guess = ntuple(k -> randn(size(T, k), r), ndims(T))

    @testset "ALS (Alternating least squares)" begin
        @time factors = candecomp(T, r, initial_guess, tol=1e-6, compute_error=true, method=:ALS)
        @test length(factors.factors) == ndims(T)
        @test size.(factors.factors) == ntuple(i -> (size(T, i), r), ndims(T))
        @test rank(factors) == r
        @test rel_residue(factors) < 1e-5
    end

    @testset "SGSD (Simultaneous generalized Schur decomposition)" begin
        @time factors = candecomp(T, r, initial_guess, compute_error=true, method=:SGSD)
        @test length(factors.factors) == ndims(T)
        @test size.(factors.factors) == ntuple(i -> (size(T, i), r), ndims(T))
        @test rank(factors) == r
        @test length(factors.lambdas) == r
        @test_broken rel_residue(factors) < 1e-5
    end

end


text:
"""
    get_next_type(T)
Return a data with a Integer Datatype twice the size of `T`
"""
function get_next_type(T::DataType)
    (T== UInt8) && return  (UInt16,8)
    (T== UInt16) && return (UInt32,16)
    (T== UInt32) && return (UInt64,32)
    (T== UInt64) && return (UInt128,64)
end

"""
    kernel_1(ij)
Create Graph from edge list `ij` based on kernel_1 benchmark(http://graph500.org/?page_id=12#sec-4)
Return `LightGraphs.SimpleGraph` structure
"""
function kernel_1(
  ij::Matrix{T}
  ) where T<:Integer

  # Creating Graph
  nv = maximum(ij) # Number of vertices
  g = Graph{T}(nv) # Creating a graph with `nv` number of vertices
  T_upper, sftamt = get_next_type(T)
  v = Vector{T_upper}()
  vertex_degree = zeros(T_upper, nv)  # Edgefactor could not exceed 2^(SCALE)...Increase SCALE
  sizehint!(v, size(ij,2)*2)
  ne = 0

  # Construct an integer with src as MSB and dst as LSB
  for i in 1:size(ij,2)
    if(ij[1,i] != ij[2,i])
        push!(v,  (T_upper(ij[1,i]) << sftamt) + ij[2,i])
        push!(v,  (T_upper(ij[2,i]) << sftamt) + ij[1,i])
        vertex_degree[T_upper(ij[1,i])]+=1
        vertex_degree[T_upper(ij[2,i])]+=1
    end
  end

 sort!(v)

 for i in one(T):nv
      sizehint!(g.fadjlist[i], vertex_degree[i])
 end

 # Create adjacency list and check for duplicate vertices
 last_added = zero(Int128)
 for edge in v
     if(edge != last_added)
         push!(g.fadjlist[T(edge >> sftamt)], T(edge & typemax(T)))
         last_added = edge
         g.ne += 1
     end
 end

  g.ne /= 2
  return g
end


text:
# You ask the submarine to determine the best route out of the deep-sea cave,
# but it only replies:
#
# Syntax error in navigation subsystem on line: all of them
#
# All of them?! The damage is worse than you thought. You bring up a copy of
# the navigation subsystem (your puzzle input).
#
# The navigation subsystem syntax is made of several lines containing chunks.
# There are one or more chunks on each line, and chunks contain zero or more
# other chunks. Adjacent chunks are not separated by any delimiter; if one
# chunk stops, the next chunk (if any) can immediately start. Every chunk must
# open and close with one of four legal pairs of matching characters:
#
# If a chunk opens with (, it must close with ).
# If a chunk opens with [, it must close with ].
# If a chunk opens with {, it must close with }.
# If a chunk opens with <, it must close with >.
#
# So, () is a legal chunk that contains no other chunks, as is []. More complex
# but valid chunks include ([]), {()()()}, <([{}])>, [<>({}){}[([])<>]], and
# even (((((((((()))))))))).
#
# Some lines are incomplete, but others are corrupted. Find and discard the
# corrupted lines first.
#
# A corrupted line is one where a chunk closes with the wrong character - that
# is, where the characters it opens and closes with do not form one of the four
# legal pairs listed above.
#
# Examples of corrupted chunks include (], {()()()>, (((()))}, and
# <([]){()}[{}]). Such a chunk can appear anywhere within a line, and its
# presence causes the whole line to be considered corrupted.
#
# For example, consider the following navigation subsystem:
#
# [({(<(())[]>[[{[]{<()<>>
# [(()[<>])]({[<{<<[]>>(
# {([(<{}[<>[]}>{[]{[(<()>
# (((({<>}<{<{<>}{[]{[]{}
# [[<[([]))<([[{}[[()]]]
# [{[{({}]{}}([{[{{{}}([]
# {<[[]]>}<{[{[{[]{()[[[]
# [<(<(<(<{}))><([]([]()
# <{([([[(<>()){}]>(<<{{
# <{([{{}}[<[[[<>{}]]]>[]]
#
# Some of the lines aren't corrupted, just incomplete; you can ignore these
# lines for now. The remaining five lines are corrupted:
#
# {([(<{}[<>[]}>{[]{[(<()> - Expected ], but found } instead.
# [[<[([]))<([[{}[[()]]] - Expected ], but found ) instead.
# [{[{({}]{}}([{[{{{}}([] - Expected ), but found ] instead.
# [<(<(<(<{}))><([]([]() - Expected >, but found ) instead.
# <{([([[(<>()){}]>(<<{{ - Expected ], but found > instead.
#
# Stop at the first incorrect closing character on each corrupted line.
#
# Did you know that syntax checkers actually have contests to see who can get
# the high score for syntax errors in a file? It's true! To calculate the
# syntax error score for a line, take the first illegal character on the line
# and look it up in the following table:
#
# ): 3 points.
# ]: 57 points.
# }: 1197 points.
# >: 25137 points.
#
# In the above example, an illegal ) was found twice (2*3 = 6 points), an
# illegal ] was found once (57 points), an illegal } was found once (1197
# points), and an illegal > was found once (25137 points). So, the total syntax
# error score for this file is 6+57+1197+25137 = 26397 points!
#
# Find the first illegal character in each corrupted line of the navigation
# subsystem. What is the total syntax error score for those errors?

const scores = Dict(')' => 3, ']' => 57, '}' => 1197, '>' => 25137)
const match = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
stack = Array{Char}(undef, 100)
cnt = foldl(eachline("input.txt"); init=0) do a,l
  p = 1
  for c in l
    if haskey(match, c)
      stack[p] = match[c]
      p += 1
    elseif p > 1 && stack[p-1] == c
      p -= 1
    else
      return a + scores[c]
    end
  end
  return a
end

println("Result: ", cnt)


text:
"""
    `function dftrace{T}(sono::Array{T, 2}, syllables::AbstractArray{SyllableBox, 1})`

Given a sonogram sono and a list of boxes around syllables that it contains,
return the dominant frequency as a function of time in each syllable

Parameters:
-----------
 - `sono` The sonogram
 - `syllables` The list of boxes around syllables

Return
------
 - An array with one member for each element of `syllables`. Each member is an array containing the index (along dimension 1) in `sono` of the dominant frequency  at each time slice in that particular syllable

    `function dftrace{T}(sono::Array{T, 2}, syllables::AbstractArray{Syllable, 1})`

A convenience method that allows syllables to be a an array of `Syllable`s rather than an array
of boxes around syllables

**Note**: see Birdsong.aliased_types for documentation on the `Syllable` and `SyllableBox` types
"""
function dftrace{T, T2}(sono::Array{T, 2}, syllables::AbstractArray{SyllableBox{T2}, 1})
    traces = similar(syllables, Array{Int, 1})

    for (i, syll) in enumerate(syllables)
	tmin = syll.tmin
	tmax = syll.tmax
	fmin = syll.hmin
	fmax = syll.hmax

	traces[i] = Array{Int, 1}((tmax-tmin)+1)

	for (j, jnd) in enumerate(tmin:tmax)
	    traces[i][j] = fmin+findmax(sono[fmin:fmax, jnd])[2]-1
	end
    end

    return traces
end

function dftrace{T, T2}(sono::Array{T, 2}, syllables::AbstractArray{Syllable{T2}, 1})
    return dftrace(sono, syllable_boxes(syllables))
end

export dftrace


text:
using Test
using DEBStandardModel, Plots  #NB: for test to pass, GR backend required
gr()

#DEBStandardModel.rundemo("../src/Msplendens_AmPfittedparams.csv", (0.0, 3.0))

@test typeof(DEBStandardModel.makeparsdict("src/Msplendens_AmPfittedparams.csv")) == Dict{SubString{String}, Real}
#@test typeof(DEBStandardModel.makeparsdict("../src/Msplendens_AmPfittedparams.csv")) == Dict
@test  typeof(DEBStandardModel.rundemo("src/Msplendens_AmPfittedparams.csv", (0.0, 3.0)))  == Plots.Plot{Plots.GRBackend}

max_stars_repo_name: HenriLaurie/DEBStandardModel

max_issues_repo_name: HenriLaurie/DEBStandardModel

max_forks_repo_name: HenriLaurie/DEBStandardModel


text:
using CSV # This is a pacakge we use for loading CSV Files.
using DataFrames
@sk_import model_selection:train_test_split
@sk_import neural_network:MLPClassifier
@sk_import metrics:confusion_matrix
@sk_import metrics:classification_report
@sk_import model_selection:cross_val_score
@sk_import metrics:plot_confusion_matrix
pathtodata = joinpath("julia-scripts", "model-zoo", "covid_cleaned.csv")
data = CSV.read(pathtodata, DataFrame)

X = convert(Array, data[!, Not(:covid_res)])
y = convert(Array, data[!, :covid_res])

X_train, X_test, y_train, y_test = train_test_split(X, y; test_size=0.33, random_state=42)

mlp_6layer = MLPClassifier(; hidden_layer_sizes=(30, 50, 60, 10, 10, 10))

fit!(mlp_6layer, X_train, y_train)

y_pred = predict(mlp_6layer, X_train)
print(classification_report(y_train, y_pred))
cross_val_score(
    MLPClassifier(; hidden_layer_sizes=(30, 50, 60, 10, 10, 10)), X_train, y_train
)
using PyPlot
plot_confusion_matrix(mlp_6layer, X_train, y_train)
p1 = gcf()

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
using Mimi

@defcomp grosseconomy begin
    regions = Index() # a regional index is defined here

    YGROSS = Variable(index = [time, regions]) # Gross output
    K = Variable(index = [time, regions]) # Capital
    l = Parameter(index = [time, regions]) # Labor
    tfp = Parameter(index = [time, regions]) # TFP
    s = Parameter(index = [time, regions]) # Savings rate
    depk = Parameter(index = [regions]) # Depreciation rate on capital
    k0 = Parameter(index = [regions]) # Initial level of capital
    share = Parameter() # Capital share

    function run_timestep(p, v, d, t)
        # the regional dimension is defined in d and parameters and variables indexed

        # define an equation for K
        for r in d.regions
            if is_first(t)
                v.K[t,r] = p.k0[r]
            else
                v.K[t,r] = (1 - p.depk[r])^5 * v.K[t-1, r] + v.YGROSS[t-1, r] * p.s[t-1, r] * 5
            end
        end

        # Define an equation for YGROSS
        for r in d.regions
            v.YGROSS[t,r] = p.tfp[t,r] * v.K[t,r]^p.share * p.l[t,r]^(1-p.share)
        end
    end
end

max_stars_repo_name: tlorans/SustainableFinance

max_issues_repo_name: tlorans/SustainableFinance

max_forks_repo_name: tlorans/SustainableFinance


text:
## Esta función da la distancia de seguridad que debe de haber entre dos vehículos
function distancia_seguridad(i, Carretera1, N, delta, alfa::Float32 = 0.75f0)

    pvfo = Vehiculo_Adelante(i, Carretera1, N)
    d_ifo = (pvfo < N ? (pvfo-i-1) : delta*5 ) # 5 = vmax[2]
    dfo = d_ifo + int8( floor( (1.-alfa) * Carretera1[pvfo].velocidad + 0.5) )
    alfa = 0; pvfo = 0; d_ifo = 0
    return dfo
end

## Esta función hace que un vehículo cambie de carril
function Cambiar_Carril!(v0::Vehiculo, Cuenta0, C1::Carretera1D)

    i = v0.posicion
    if C1.carretera[i].tipo == -1

        pvbo = Vehiculo_Atras(i, C1.carretera, C1.N)
        ############################ cambio ###################################
        d_bo = (pvbo < (C1.N) ? (i-pvbo-1) : 5 ) # 5 = vmax[2] (i-pvbo-longitud)
        vb = C1.carretera[pvbo].velocidad
        if vb <= d_bo

            Cambiar_Vehiculo!( C1.carretera[i], v0.velocidad, v0.posicion, v0.tipo, 1, v0.num)
            Celda_Vacia!( v0 )
            Cuenta0 -= 1
            C1.cuenta += 1
        end
    end
end

## Esta función hace los cambios de carril desde un carril izquierdo hacia uno derecho
function CambioIzq_Der!(Ck, Ckmenos1, k)

    for v in Ck.carretera[end:-1:9]
        if (v.tipo == 1 || v.tipo == 2) && v.cambio == 0

            delta = v.tipo == 2 ? int8(3) : int8(1)
            l = v.velocidad
            dfo = distancia_seguridad(v.posicion, Ckmenos1.carretera, Ckmenos1.N, delta)
            dfs = distancia_seguridad(v.posicion, Ck.carretera, Ck.N, delta)

            if v.tipo == 2 && dfo >= l*delta && (dfs  >= delta*l || l > dfs)
                Cambiar_Carril!(v, Ck.cuenta, Ckmenos1)
            end

            if v.tipo == 1 && dfo >= l*delta #(dfs  >= delta*l || l > dfs)
                Cambiar_Carril!(v, Ck.cuenta, Ckmenos1)
            end
            l = dfo = dfs = delta = 0
        end
    end
end

## Esta función hace los cambios de carril desde un carril derecho hacia uno izquierdo
function CambioDer_Izq!(Ck, Ckmas1, k)

    for v in Ck.carretera[end:-1:9]
        #if (v.tipo == 2 || (v.tipo == 1 && k < 2)) && v.cambio == 0
        if v.tipo == 2 && v.cambio == 0
            l = v.velocidad
            dfs = distancia_seguridad(v.posicion, Ck.carretera, Ck.N, int8(1))

            if dfs < l
                if Ckmas1.carretera[v.posicion].tipo == -1
                    dfo = distancia_seguridad(v.posicion, Ckmas1.carretera, Ckmas1.N, int8(1))
                    if dfo >= l
                        Cambiar_Carril!(Ck.carretera[v.posicion], Ck.cuenta, Ckmas1)
                    end
                end
            end
            dfs = l = dfo = 0
        end
    end
end

## Esta función simula una rampa de entrada (tipo == 1) o de salida (tipo == 0)
function Rampa!( tipo, x0, lramp, pin, p, Ck, num, vmax::Array{Int8, 1} = Int8[3, 5])
    #vmax = [3, 5]

    # Rampa de salida, tipo == 0
    if tipo == 0
    ############################ cambio ###################################
        x = Pos_der(Ck.carretera[x0:x0+lramp])  # if longitud > 1 Ck.carretera[xi] && ... && Ck.carretera[xi+londitud] esten vacia
        if x <= x0+lramp && rand() < pin
            if Ck.cuenta > 0
                Ck.cuenta -= 1
            end
            Celda_Vacia!(Ck.carretera[x])
        end

    # Rampa de entrada, tipo == 1
    elseif tipo == 1
    ############################ cambio ###################################
        x = Pos_siguiente(Ck.carretera[x0:x0+lramp])
        if x <= x0+lramp && rand() < pin
          Ck.cuenta += 1
          v_nueva = (rand() <= p ? vmax[1] : vmax[2])
          tipo_nuevo = (v_nueva == vmax[1] ? 1 : 2)
          Cambiar_Vehiculo!( Ck.carretera[x], v_nueva, x, tipo_nuevo, 1, num)
        end
        v_nueva = 0; tipo_nuevo = 0
    end
    x = 0
end


text:
"""
Sampling from DPP and k-DPP.

Methods:
--------
    rand(dpp::DeterminantalPointProcess, N::Int)
    randmcmc(dpp::DeterminantalPointProcess, N::Int)
    rand(dpp::DeterminantalPointProcess, N::Int, k::Int)
    randmcmc(dpp::DeterminantalPointProcess, N::Int, k::Int)

References:
-----------
    [1] Kulesza, A., and B. Taskar. Determinantal point processes for machine
        learning. arXiv preprint arXiv:1207.6083, 2012.
    [2] Kang, B. Fast determinantal point process sampling with application to
        clustering. NIPS, 2013.
"""

"""Sample a mask for an elementary DPP.
"""
function _sample_mask(Λ::SharedArray{Float64},
                      M::SharedMatrix{Bool},
                      i::Int, seed::Int)
    rng = MersenneTwister(seed)

    for j in 1:length(Λ)
        M[j, i] = (rand(rng) < Λ[j] / (Λ[j] + 1))
    end
end

"""Sample a mask for an elementary k-DPP.
"""
function _sample_k_mask(Λ::SharedArray{Float64},
                        M::SharedMatrix{Bool},
                        E::SharedMatrix{Float64},
                        k::Int, i::Int, seed::Int)
    rng = MersenneTwister(seed)

    j = length(Λ)
    remaining = k

    # iteratively sample a k-mask
    while remaining > 0
        # compute marginal of j given that we choose remaining values from 1:j
        if j == remaining
            marg = 1
        else
            marg = Λ[j] * E[remaining, j] / E[remaining + 1, j + 1];
        end

        # sample marginal
        if rand(rng) <= marg
            M[j, i] = true
            remaining -= 1
        end
        j -= 1
      end
end

"""Exact sampling from an elementary DPP. The algorithm based on [1].
"""
function _sample_from_elementary(V::SharedMatrix,
                                 M::SharedMatrix{Bool},
                                 i::Int, seed::Int)
    rng = MersenneTwister(seed)

    # select the elementary DPP
    V_mask = M[:, i]

    # edge case: empty sample
    if !any(V_mask)
        return Int[]
    end

    # select the kernel of the elementary DPP
    L = V[:, V_mask]

    Y = Int[]
    mask = ones(Bool, size(L, 2))
    prob = Array{Float64}(undef,size(L, 1))

    for i in 1:size(L, 2)
        # compute probabilities
        fill!(prob, 0)
        for c in 1:size(L, 2)
            !mask[c] && continue
            for r in 1:size(L, 1)
                prob[r] += L[r, c].^2
            end
        end
        prob ./= sum(prob)

        # sample a point in the original space
        h = findfirst(rand(rng) .<= cumsum(prob))
        push!(Y, h)

        # select and mask-out an element
        j = get_first_nz_idx(L[h, :], mask)
        mask[j] = false

        if any(mask)
            # Subtract scaled Lj from other columns so that their
            # projections on e_s[i] turns into 0. This operation
            # preserves the rank of L_{-j}.
            for c in 1:size(L, 2)
                !mask[c] && continue
                for r in 1:size(L, 1)
                    L[r, c] -= L[r, j] * L[h, c] / L[h, j]
                end
            end

            # Gram-Schmidt orthogonalization
            L[:, mask] = Matrix(qr(L[:, mask]).Q)
        end
    end

    sort(Y)
end

"""Exact sampling from a DPP [1].
"""
function Base.rand(dpp::DeterminantalPointProcess, N::Int)
    Λ = SharedArray{Float64}(dpp.Lfact.values)
    V = SharedMatrix{Float64}(dpp.Lfact.vectors)
    M = SharedMatrix{Bool}(zeros(Bool, dpp.size, N))

    # step I: sample masks for elementary DPPs
    pmap((i, seed) -> _sample_mask(Λ, M, i, seed),
         1:N, abs.(rand(dpp.rng, Int, N)))

    # step II: iteratively sample from a mixture of elementary DPPs
    pmap((i, seed) -> _sample_from_elementary(V, M, i, seed),
         1:N, abs.(rand(dpp.rng, Int, N)))
end

"""Exact sampling from a k-DPP [1].
"""
function Base.rand(dpp::DeterminantalPointProcess, N::Int, k::Int)
    Λ = SharedArray{Float64}(dpp.Lfact.values)
    V = SharedMatrix{Float64}(dpp.Lfact.vectors)
    M = SharedMatrix{Bool}(zeros(Bool, dpp.size, N))

    # compute elementary symmetric polynomials
    E = SharedMatrix{Float64}(elem_symm_poly(dpp.Lfact.values, k))

    # step I: sample masks for elementary DPPs
    pmap((i, seed) -> _sample_k_mask(Λ, M, E, k, i, seed),
         1:N, abs.(rand(dpp.rng, Int, N)))

    # step II: iteratively sample from a mixture of elementary DPPs
    pmap((i, seed) -> _sample_from_elementary(V, M, i, seed),
         1:N, abs.(rand(dpp.rng, Int, N)))
end

"""Perform one MCMC accept-reject transition for DPP.
"""
function _do_mcmc_step!(dpp::DeterminantalPointProcess, state::MCMCState)
    # propose an element to swap
    u = rand(dpp.rng, 1:dpp.size)
    insert = !state[1][u]

    # attempt to make a transition
    if insert
        p = _comp_accept_prob(dpp, state, u, insert)
        rand(dpp.rng) < p && _update_mcmc_state!(dpp, state, u, insert)
    else  # delete
        new_state = _update_mcmc_state(dpp, state, u, insert)
        p = _comp_accept_prob(dpp, new_state, u, insert)
        if rand(dpp.rng) < p
            state[1] .= new_state[1]
            state[2] .= new_state[2]
        end
    end
end

"""Perform one MCMC accept-reject transition for k-DPP.
"""
function _do_mcmc_k_step!(dpp::DeterminantalPointProcess, state::MCMCState)
    z, L_z_inv = state

    # propose the elements to swap
    u, v = rand(findall(z)), rand(findall(z.==true))

    # copy the state and delete the u element
    new_state = _update_mcmc_state(dpp, state, u, false)

    # attempt to make a transition
    p = _comp_accept_prob(dpp, new_state, u, v)
    if rand(dpp.rng) < p
        # insert the v element into the new state
        _update_mcmc_state!(dpp, new_state, v, true)
        state[1] .= new_state[1]
        state[2] .= new_state[2]
    end
end

function _comp_accept_prob(dpp::DeterminantalPointProcess, state::MCMCState,
                           u::Int, insert::Bool)
    """Compute accept probability to insert / delete u from the state.
    """
    z, L_z_inv = state

    d_u = dpp.L[u, u]
    if any(z)
        b_u = dpp.L[z, u]
        d_u -= dot(b_u, L_z_inv[z, z] * b_u)
    end

    insert ? min(1.0, d_u) : min(1.0, 1.0 / d_u)
end

function _comp_accept_prob(dpp::DeterminantalPointProcess, state::MCMCState,
                           u::Int, v::Int)
    """Compute accept probability to swap u and v.
    """
    z, L_z_inv = state

    d_u, d_v = dpp.L[u, u], dpp.L[v, v]
    if any(z)
        b_u, b_v = dpp.L[z, u], dpp.L[z, v]
        d_u -= dot(b_u, L_z_inv[z, z] * b_u)
        d_v -= dot(b_v, L_z_inv[z, z] * b_v)
    end

    min(1.0, d_v / d_u)
end

function _update_mcmc_state!(dpp::DeterminantalPointProcess, state::MCMCState,
                             u::Int, insert::Bool)
    """Update the state after u is inserted / deleted.
    """
    z, L_z_inv = state

    if insert
        d_u = dpp.L[u, u]
        if any(z)
            b_u = dpp.L[z, u]
            x_u = L_z_inv[z, z] * b_u
            d_u -= dot(b_u, x_u)

            L_z_inv[z, z] += (x_u * x_u') / d_u
            L_z_inv[z, u] = L_z_inv[u, z] = - x_u / d_u
        end

        L_z_inv[u, u] = 1. / d_u
        z[u] = true
    else  # delete
        z[u] = false

        e = L_z_inv[z, u]
        f = L_z_inv[u, u]

        L_z_inv[z, z] -= (e * e') / f
    end
end

function _update_mcmc_state(dpp::DeterminantalPointProcess, state::MCMCState,
                            u::Int, insert::Bool)
    """Update the state after u is inserted / deleted.
    """
    new_state = deepcopy(state)
    _update_mcmc_state!(dpp, new_state, u, insert)
    new_state
end

function randmcmc(dpp::DeterminantalPointProcess, N::Int;
                  init_state=nothing,
                  return_final_state::Bool=false,
                  mix_eps::Float64=1e-1,
                  mixing_steps::Int=ceil(Int, dpp.size*log(dpp.size/mix_eps)),
                  steps_between_samples::Int=mixing_steps)
    """MCMC sampling from a DPP [2].

    TODO: Add support for running MCMC in parallel, similar as rand.
          Make sure parallelization produces unbiased and consistent samples.
    """
    # initialize the Markov chain
    state = init_state
    if state == nothing
        L_z_inv = Array{Float64}(undef,size(dpp.L))
        z = bitrand(dpp.rng, dpp.size)  # TODO: improve initialization (?)
        if any(z)
            L_z_inv[z, z] = pinv(dpp.L[z, z])
        end
        state = (z, L_z_inv)
    end

    # sanity check
    @assert typeof(state) == MCMCState

    # mix the Markov chain
    for t in 1:mixing_steps
        _do_mcmc_step!(dpp, state)
    end

    Y = []
    for i in 1:N
        push!(Y, findall(state[1]))
        for t in 1:steps_between_samples
            _do_mcmc_step!(dpp, state)
        end
    end

    return_final_state ? (Y, state) : Y
end

function randmcmc(dpp::DeterminantalPointProcess, N::Int, k::Int;
                  init_state=nothing,
                  return_final_state::Bool=false,
                  mix_eps::Float64=1e-1,
                  mixing_steps::Int=ceil(Int, k*log(k / mix_eps)),
                  steps_between_samples::Int=mixing_steps)
    """MCMC sampling from a k-DPP [2].

    TODO: Add support for running MCMC in parallel, similar as rand.
          Make sure parallelization produces unbiased and consistent samples.
    """
    # initialize the Markov chain
    state = init_state
    if state == nothing
        L_z_inv = Array{Float64}(undef,size(dpp.L))
        z = falses(dpp.size)  # TODO: improve initialization (?)
        z[1:k] .= true
        if any(z)
            L_z_inv[z, z] = pinv(dpp.L[z, z])
        end
        state = (z, L_z_inv)
    end

    # sanity check
    @assert typeof(state) == MCMCState
    @assert sum(state[1]) == k

    # mix the Markov chain
    for t in 1:mixing_steps
        _do_mcmc_k_step!(dpp, state)
    end

    Y = []
    for i in 1:N
        push!(Y, findall(state[1]))
        for t in 1:steps_between_samples
            _do_mcmc_k_step!(dpp, state)
        end
    end

    return_final_state ? (Y, state) : Y
end


text:
@testset "clean" begin
    remove(join(TEST_TMP_DIR, "test-playground"), recursive=true)
    @test !exists(join(TEST_TMP_DIR, "test-playground"))
    @test islink(join(TEST_CONFIG.bin, "julia-bin"))
    @test islink(join(TEST_CONFIG.share, "myproject"))

    clean(TEST_CONFIG)
    @test !exists(join(TEST_TMP_DIR, "test-playground"))
    @test islink(join(TEST_CONFIG.bin, "julia-bin"))
    @test !islink(join(TEST_CONFIG.share, "myproject"))
end

@testset "remove" begin
    rm(TEST_CONFIG, dir=join(TEST_DIR, ".playground"))
    @test !exists(join(TEST_DIR, ".playground"))

    rm(TEST_CONFIG, name="julia-nightly-dir")
    @test !exists(join(TEST_CONFIG.bin, "julia-nightly-dir"))

    remove(join(TEST_TMP_DIR, "test-playground2"), recursive=true)
    rm(TEST_CONFIG, name="otherproject")
    @test !islink(join(TEST_CONFIG.share, "otherproject"))
    @test !exists(join(TEST_TMP_DIR, "test-playground2"))
end


text:
# https://sudonull.com/post/7909-Julia-Scripts-and-parsing-command-line-arguments
# import Pkg; Pkg.add("PyPlot")
using DelimitedFiles, SharedArrays, Random, LinearAlgebra, Clustering, Base
using DataFrames, CSV
using Statistics, ArgParse

nmats = 100

function max_eigv(gene) #returns largest eigenvalue of a randomized matrix generated from data.
    gran = zeros(size(gene))
    for i = 1:size(gene,2)
        gran[:,i] = gene[shuffle(1:end),i]
    end
 #    @time mr = corr_mat(gran)
    # println(mr[[end-1],end])
    mr = cor(transpose(gran))
    # println(mr[[end-1],end])
    return eigvals(mr)
end

parsed_args = Dict([("c", true), ("expfile", "Data/clean/Basal-clean.tsv")])

##============ 0.5 - Input
#! gname s1 s2 ... gstart chromname
df = CSV.read(parsed_args["expfile"], delim='\t')

df = df[df.chromname .== "17", :]
gene = convert(Array{Float64,2},df[:, 2:end-2])

evs = SharedArray{Float64, 2}((nmats,size(gene)[1]))
@time for i in 1:nmats
    @time evs[i,:] = max_eigv(gene) # return array of ngenes size
end
# All randomized eigenvalues
new_evs = vec(evs)
new_evs = filter(x->x>0.0000001, new_evs) #Drop Zeros

# Original eigenvalues
mr = cor(transpose(gene))
mr_e = eigvals(mr)
mr_e = filter(x->x>0.0000001, mr_e) #Drop Zeros

using PyPlot

#Figura 2a
plt.clf()
# color = '0.75'
# color = '#eeefff'
plt.hist(new_evs,density=true,bins=100,label="Randomized",color=[:C1])
# grid("on")
xlabel("Eigenvalues")
ylabel("Probability density")
legend()
ax = gca()
ax.spines["top"].set_visible(false)
ax.spines["right"].set_visible(false)
plot(size=(1200,800))
savefig("Plots/Basal-ch17-eigendist.png")

#Figura 2b
plt.clf()
plt.hist(mr_e,density=true,bins=100,label="Empirical data")
plt.hist(new_evs,density=true,bins=100,alpha=0.7,label="Randomized")
# grid("on")
xlabel("Eigenvalues")
ylabel("Probability density")
legend()
ax = gca()
ax.spines["top"].set_visible(false)
ax.spines["right"].set_visible(false)
plot(size=(1200,800))
savefig("Plots/Basal-ch17-eigenRandEmpirical.png")

max_stars_repo_name: josemaz/gene-matrices

max_issues_repo_name: josemaz/gene-matrices

max_forks_repo_name: josemaz/gene-matrices


text:
include("SDOF_Model.jl")

using Plots, Plots.PlotMeasures, LaTeXStrings

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

# ## Setup

p = sdof()

prob = IVP(p)
q = InitialValueProblem_quad(p)

x0 = initial_state(prob)

T = period(p)

tmax = 4.0
dom = 0:0.001:tmax

# some plotting options
COLOR_initial = :green
COLOR_transformed = :yellow
COLOR_convexhull = :blue
COLOR_Ω0 = :orange1
COLOR_Ω0_box = :lightblue
COLOR_analytic = :magenta
COLOR_numeric = :red
COLOR_zono = :chartreuse;

#-

# ## Singleton initial conditions

# ### Discretization

model = StepIntersect(setops=:concrete)

α = 0.05
δ = α*T

probn = normalize(prob)
probd = discretize(probn, δ, model)

Φ = state_matrix(probd)
Ω0 = initial_state(probd);

function plot_singleton_discretization()
    fig = plot(legend=:topleft, xlab=L"u(t)", ylab=L"v(t)",
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([0.95, 0.96, 0.97, 0.98, 0.99, 1.0], [L"0.95", L"0.96", L"0.97", L"0.98", L"0.99", L"1.0"]),
               ytick=([0.0, -1.0, -2.0, -3.0, -4.0], [L"0.0", L"-1.0", L"-2.0", L"-3.0", L"-4.0"]),
               xlims=(0.94, 1.0),
               ylims=(-4.0, 1.0),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, size=(900, 600))
    plot!(fig, box_approximation(Ω0), lab=L"\textrm{Box}(\Omega_0)", c=COLOR_Ω0_box)
    plot!(fig, Ω0, lab=L"\Omega_0", c=COLOR_Ω0, alpha=0.6)

    ## analytic solution of the ODE
    tdom1 = 0:0.0001:α*T

    asol_x = analytic_solution(p)
    asol_v = analytic_derivative(p)

    ## refine the domain of interest
    plot!(fig, asol_x.(tdom1), asol_v.(tdom1), seriestype=:path, c=COLOR_analytic, lw=1.5, lab=L"\textrm{Analytic} = \Omega_0^e")

    ## initial state: singleton
    plot!(fig, x0, c=COLOR_initial, markersize=12, alpha=0.5, markershape=:circle, lab=L"\textrm{Initial } \mathcal{X}_0")

    ## transformed state: linear map of the singleton
    plot!(fig, Φ * x0, c=COLOR_transformed, markersize=12, alpha=0.5, markershape=:circle, lab=L"\textrm{Transformation }\Phi \mathcal{X}_0")

    ## initial state: singleton
    plot!(fig, CH(x0, Φ * x0), c=COLOR_convexhull, markersize=0, linestyle=:solid, lw=2.0, markershape=:circle, lab=L"\textrm{Convex hull CH}(\mathcal{X}_0, \Phi \mathcal{X}_0)")

    ## numerical solution with smaller step size
    algn = Newmark(Δt=α*T/5, α=0.25, δ=0.5)
    soln = solve(q, algn, NSTEPS=round(Int, tmax / (α*T)));
    plot!(fig, soln.U[1:10], soln.U′[1:10], markershape=:utriangle, markersize=8, color=COLOR_numeric, lab="")
    plot!(fig, soln.U[11:11], soln.U′[11:11], markershape=:utriangle, markersize=8, color=COLOR_numeric, lab=L"\textrm{Numeric}")

    return fig
end

fig = plot_singleton_discretization()
savefig(fig, joinpath(TARGET_FOLDER, "sdof_singleton_discretization.pdf"))

# ### Flowpipe

model = StepIntersect(setops=:concrete)
α = 0.05

alg = VREP(δ=α*T, approx_model=model)
solvrep = solve(prob, tspan=(0.0, tmax), alg=alg);

alg = GLGM06(δ=α*T, approx_model=model)
solglg = solve(prob, tspan=(0.0, tmax), alg=alg);

alg = BOX(δ=α*T, approx_model=model)
solbox = solve(prob, tspan=(0.0, tmax), alg=alg);

function plot_singleton_flowpipe()
    fig = plot(legend=:topleft, xlab=L"u(t)", ylab=L"v(t)",
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([-1.0, -0.5, 0.0, 0.5, 1.0], [L"-1.0", L"-0.5", L"0.0", L"0.5", L"1.0"]),
               ytick=([0.0, -3.0, -6.0, -9.0, -12.0], [L"0.0", L"-3.0", L"-6.0", L"-9.0", L"-12.0"]),
               xlims=(-1.1, 1.2),
               ylims=(-14, 1.4),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=8mm, size=(900, 600))

    NMAX = 9

    ## propagation of Ω0 without overapproximation
    plot!(fig, solbox[1:NMAX], vars=(1, 2), c=COLOR_Ω0_box, lab=L"\textrm{Flowpipe (hyperrectangles)}")

    ## propagation of Ω0 using zonotopes
    plot!(fig, solglg[1:NMAX], vars=(1, 2), c=COLOR_zono, lw=1.0, alpha=.6, lab=L"\textrm{Flowpipe (zonotopes)}")

    ## propagation of Ω0 withot overapproximation
    plot!(fig, solvrep[1:NMAX], vars=(1, 2), c=COLOR_Ω0, alpha=0.8, lw=1.0, lab=L"\textrm{Flowpipe (polygons)}")

    ## analytic solution of the ODE
    tdom1 = 0:0.0001:NMAX*α*T
    asol_x = analytic_solution(p)
    asol_v = analytic_derivative(p)
    ## refine the domain of interest
    plot!(fig, asol_x.(tdom1), asol_v.(tdom1), seriestype=:path, c=COLOR_analytic, lw=2.0, lab=L"\textrm{Analytic} = F^e")

    ## numerical solution with smaller step size
    algn = Newmark(Δt=α*T/5, α=0.25, δ=0.5)
    soln = solve(q, algn, NSTEPS=round(Int, tmax / (α*T)));
    plot!(fig, soln.U[1:1:NMAX*5], soln.U′[1:1:NMAX*5], markersize=3.0, markershape=:utriangle, color=COLOR_numeric, lab="")
    plot!(fig, soln.U[NMAX*10+1:NMAX*10+1], soln.U′[NMAX*10+1:NMAX*10+1], markershape=:utriangle, color=COLOR_numeric, lab=L"\textrm{Numeric}")
    ##  \Delta t = \delta / 5

    fig
end

fig = plot_singleton_flowpipe()
savefig(fig, joinpath(TARGET_FOLDER, "sdof_singleton_flowpipe.pdf"))

# ## Distributed initial conditions

# ### Discretization

p = sdof()

prob = IVP(p)
q = InitialValueProblem_quad(p)

X0 = Hyperrectangle(low=[0.9, -0.1], high=[1.1, 0.1])

probset = IVP(prob.s, X0);

T = period(p)

tmax = 4.0
dom = 0:0.001:tmax

model = StepIntersect(setops=:concrete)

α = 0.05
δ = α*T

probset_d = discretize(normalize(probset), δ, model)

Φ = state_matrix(probset_d)
Ω0 = initial_state(probset_d);

#-

import Random

Random.seed!(3333)

function plot_distributed_discretization()
    fig = plot(legend=:left, xlab=L"u(t)", ylab=L"v(t)",
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([0.85, 0.90, 0.95, 1.00, 1.05, 1.10], [L"0.85", L"0.90", L"0.95", L"1.00", L"1.05", L"1.10"]),
               ytick=([0.0, -1.0, -2.0, -3.0, -4.0], [L"0.0", L"-1.0", L"-2.0", L"-3.0", L"-4.0"]),
               xlims=(0.84, 1.12),
               ylims=(-4.5, 0.5),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, size=(900, 600))

    plot!(fig, box_approximation(Ω0), lab=L"\textrm{Box}(\Omega_0)", c=COLOR_Ω0_box)
    plot!(fig, Ω0, lab=L"\Omega_0", c=COLOR_Ω0, alpha=0.6)

    ## convex hull of the initial state and the transformed state
    plot!(fig, CH(X0, Φ * X0), c=COLOR_convexhull, lab=L"\textrm{Convex hull CH}(\mathcal{X}_0, \Phi \mathcal{X}_0)")

    ## initial state: singleton
    plot!(fig, X0, c=COLOR_initial, alpha=1., lab=L"\textrm{Initial } \mathcal{X}_0")

    ## transformed state: linear map of the singleton
    plot!(fig, Φ * X0, c=COLOR_transformed, alpha=1., lab=L"\textrm{Transformation }\Phi \mathcal{X}_0")

    ## -------- Random sample of analytic solutions ------
    tdom1 = 0:0.0001:α*T

    nsamples = 50
    X0sampled = sample(X0, nsamples);
    ω = p.ω

    s = X0sampled[1]
    asol_x = analytic_solution(p, x0=s[1], v0=s[2])
    asol_v = analytic_derivative(p, x0=s[1], v0=s[2])
    plot!(fig, asol_x.(tdom1), asol_v.(tdom1), seriestype=:path, lw=1.5, c=COLOR_analytic, lab=L"\textrm{Analytic (random)}")

    for k in 2:nsamples
        s = X0sampled[k]
        asol_x = analytic_solution(p, x0=s[1], v0=s[2])
        asol_v = analytic_derivative(p, x0=s[1], v0=s[2])
        plot!(fig, asol_x.(tdom1), asol_v.(tdom1), seriestype=:path, c=COLOR_analytic, lw=1.5, lab="")
    end

    return fig
end

fig = plot_distributed_discretization()
savefig(fig, joinpath(TARGET_FOLDER, "sdof_distributed_discretization.pdf"))

# ### Flowpipe

model = StepIntersect(setops=:concrete)
α = 0.05

alg = VREP(δ=α*T, approx_model=model)
solvrepset = solve(probset, tspan=(0.0, tmax), alg=alg);

alg = GLGM06(δ=α*T, approx_model=model)
solglgset = solve(probset, tspan=(0.0, tmax), alg=alg);

alg = BOX(δ=α*T, approx_model=model)
solboxset = solve(probset, tspan=(0.0, tmax), alg=alg);

import Random

Random.seed!(3333)

function plot_distributed_flowpipe()
    fig = plot(legend=:top, xlab=L"u(t)", ylab=L"v(t)",
               legendfontsize=15,
               tickfont=font(20, "Times"),
               guidefontsize=20,
               xguidefont=font(30, "Times"),
               yguidefont=font(30, "Times"),
               xtick=([-1.0, -0.5, 0.0, 0.5, 1.0], [L"-1.0", L"-0.5", L"0.0", L"0.5", L"1.0"]),
               ytick=([0.0, -3.0, -6.0, -9.0, -12.0], [L"0.0", L"-3.0", L"-6.0", L"-9.0", L"-12.0"]),
               xlims=(-1.15, 1.2),
               ylims=(-15, 1.1),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm, size=(900, 600))
    NMAX = 9

    ## propagation of Ω0 without overapproximation
    plot!(fig, solboxset[1:NMAX], vars=(1, 2), c=COLOR_Ω0_box, lab=L"\textrm{Flowpipe (hyperrectangles)}")

    ## propagation of Ω0 using zonotopes
    plot!(fig, solglgset[1:NMAX], vars=(1, 2), c=COLOR_zono, lw=1.0, alpha=.6, lab=L"\textrm{Flowpipe (zonotopes)}")

    ## propagation of Ω0 withot overapproximation
    plot!(fig, solvrepset[1:NMAX], vars=(1, 2), c=COLOR_Ω0, alpha=0.8, lw=1.0, lab=L"\textrm{Flowpipe (polygons)}")

    ## -------- Random sample of analytic solutions ------
    tdom1 = 0:0.0001:NMAX*α*T

    nsamples = 50
    X0sampled = sample(X0, nsamples);
    ω = p.ω

    s = X0sampled[1]
    asol_x = analytic_solution(p, x0=s[1], v0=s[2])
    asol_v = analytic_derivative(p, x0=s[1], v0=s[2])
    plot!(fig, asol_x.(tdom1), asol_v.(tdom1), seriestype=:path, lw=1.0, c=COLOR_analytic, lab=L"\textrm{Analytic (random)}")

    for k in 2:nsamples
        s = X0sampled[k]
        asol_x = analytic_solution(p, x0=s[1], v0=s[2])
        asol_v = analytic_derivative(p, x0=s[1], v0=s[2])
        plot!(fig, asol_x.(tdom1), asol_v.(tdom1), seriestype=:path, c=COLOR_analytic, lw=1.0, lab="")
    end

    return fig
end

fig = plot_distributed_flowpipe()
savefig(fig, joinpath(TARGET_FOLDER, "sdof_distributed_flowpipe.pdf"))

max_stars_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_issues_repo_name: JuliaReach/SetPropagation-FEM-Examples

max_forks_repo_name: JuliaReach/SetPropagation-FEM-Examples


text:
"""
`sort_by_size(list)` returns a sorted copy of `list`
in which the elements are sorted by the `length` function
(from smallest to largest).
"""
function sort_by_size(list::Array)
    sort(list,lt=((x,y) -> length(x) < length(y)))
end

import Base.sqrt

function weave(c::Array{Int,1})::Array{Int,1}
    nc = length(c)
    if nc%2 == 0
        error("Attempting to weave $c, an even-length cycle, with itself")
    end
    half = Int((nc+1)/2)
    result = zeros(Int,nc)
    j = 1  # index into c
    jj = 1  # index into result
    while j <= nc
        result[j] = c[jj]
        j += 1
        jj = (jj+half) % nc
        if jj == 0
            jj = nc
        end
    end
    return result
end

function weave(c1::Array{Int,1}, c2::Array{Int,1})::Array{Int,1}
    n1 = length(c1)
    n2 = length(c2)
    @assert n1==n2 "Cycles $c1 and $c2 have different lengths"
    @assert n1%2==0 "Cycles must have even length to be weaved together"

    result = zeros(Int,2n1)
    for j=1:n1
        result[2j-1]=c1[j]
        result[2j]=c2[j]
    end
    return result
end
"""
`sqrt(p::Permutation)` returns a `Permutation` `q` such that
`q*q==p`. Note: There may be other square roots besides the one
returned. If `p` does not have a square root, an error is thrown.
"""
function sqrt(p::Permutation)::Permutation
    err_msg = "This permutation does not have a square root"
    n = length(p)
    clist = sort_by_size(cycles(p))
    nc = length(clist)
    result = zeros(Int,n)

    idx = 1
    while idx <= nc
        cyc = clist[idx]  # next cycle to process
        m = length(cyc)
        if m%2 == 1   # an odd cycle
            tmp = weave(cyc)
            for j=1:m-1
                result[tmp[j]] = tmp[j+1]
            end
            result[tmp[m]] = tmp[1]
            idx += 1
        else
            if idx == nc
                error(err_msg)
            end
            cyc2 = clist[idx+1]
            if length(cyc2) != length(cyc)
                error(err_msg)
            end
            tmp = weave(cyc,cyc2)
            for j=1:2m-1
                result[tmp[j]] = tmp[j+1]
            end
            result[tmp[2m]] = tmp[1]
            idx += 2
        end
    end
    return Permutation(result)
end

√(p::Permutation) = sqrt(p)