text:
# parsing tests

@test v"1" == VersionNumber(1)
@test v"2.1" == VersionNumber(2, 1)
@test v"3.2.1" == VersionNumber(3, 2, 1)

@test v"2-1" == VersionNumber(2, 0, 0, (1,), ())
@test v"3.2-1" == VersionNumber(3, 2, 0, (1,), ())
@test v"4.3.2-1" == VersionNumber(4, 3, 2, (1,), ())

@test v"1-a" == VersionNumber(1, 0, 0, ("a",), ()) == v"1a"
@test v"2.1-a" == VersionNumber(2, 1, 0, ("a",), ()) == v"2.1a"
@test v"3.2.1-a" == VersionNumber(3, 2, 1, ("a",), ()) == v"3.2.1a"

@test v"2-a1" == VersionNumber(2, 0, 0, ("a1",), ()) == v"2a1"
@test v"3.2-a1" == VersionNumber(3, 2, 0, ("a1",), ()) == v"3.2a1"
@test v"4.3.2-a1" == VersionNumber(4, 3, 2, ("a1",), ()) == v"4.3.2a1"

@test v"2-1a" == VersionNumber(2, 0, 0, ("1a",), ())
@test v"3.2-1a" == VersionNumber(3, 2, 0, ("1a",), ())
@test v"4.3.2-1a" == VersionNumber(4, 3, 2, ("1a",), ())

@test v"2-a.1" == VersionNumber(2, 0, 0, ("a", 1), ()) == v"2a.1"
@test v"3.2-a.1" == VersionNumber(3, 2, 0, ("a", 1), ()) == v"3.2a.1"
@test v"4.3.2-a.1" == VersionNumber(4, 3, 2, ("a", 1), ()) == v"4.3.2a.1"

@test v"2-1.a" == VersionNumber(2, 0, 0, (1, "a"), ())
@test v"3.2-1.a" == VersionNumber(3, 2, 0, (1, "a"), ())
@test v"4.3.2-1.a" == VersionNumber(4, 3, 2, (1, "a"), ())

@test v"2+1" == VersionNumber(2, 0, 0, (), (1,))
@test v"3.2+1" == VersionNumber(3, 2, 0, (), (1,))
@test v"4.3.2+1" == VersionNumber(4, 3, 2, (), (1,))

@test v"1+a" == VersionNumber(1, 0, 0, (), ("a",))
@test v"2.1+a" == VersionNumber(2, 1, 0, (), ("a",))
@test v"3.2.1+a" == VersionNumber(3, 2, 1, (), ("a",))

@test v"2+a1" == VersionNumber(2, 0, 0, (), ("a1",))
@test v"3.2+a1" == VersionNumber(3, 2, 0, (), ("a1",))
@test v"4.3.2+a1" == VersionNumber(4, 3, 2, (), ("a1",))

@test v"2+1a" == VersionNumber(2, 0, 0, (), ("1a",))
@test v"3.2+1a" == VersionNumber(3, 2, 0, (), ("1a",))
@test v"4.3.2+1a" == VersionNumber(4, 3, 2, (), ("1a",))

@test v"2+a.1" == VersionNumber(2, 0, 0, (), ("a", 1))
@test v"3.2+a.1" == VersionNumber(3, 2, 0, (), ("a", 1))
@test v"4.3.2+a.1" == VersionNumber(4, 3, 2, (), ("a", 1))

@test v"2+1.a" == VersionNumber(2, 0, 0, (), (1, "a"))
@test v"3.2+1.a" == VersionNumber(3, 2, 0, (), (1, "a"))
@test v"4.3.2+1.a" == VersionNumber(4, 3, 2, (), (1, "a"))

# basic comparison

@test VersionNumber(2,3,1) == VersionNumber(int8(2),uint32(3),int32(1)) == v"2.3.1"
@test v"2.3.0" < v"2.3.1" < v"2.4.8" < v"3.7.2"

# advanced comparison & manipulation

import Base: thispatch, thisminor, thismajor,
             nextpatch, nextminor, nextmajor, check_new_version

@test v"1.2.3" == thispatch(v"1.2.3-")
@test v"1.2.3" == thispatch(v"1.2.3-pre")
@test v"1.2.3" == thispatch(v"1.2.3")
@test v"1.2.3" == thispatch(v"1.2.3+post")
@test v"1.2.3" == thispatch(v"1.2.3+")

@test v"1.2" == thisminor(v"1.2.3-")
@test v"1.2" == thisminor(v"1.2.3-pre")
@test v"1.2" == thisminor(v"1.2.3")
@test v"1.2" == thisminor(v"1.2.3+post")
@test v"1.2" == thisminor(v"1.2.3+")

@test v"1" == thismajor(v"1.2.3-")
@test v"1" == thismajor(v"1.2.3-pre")
@test v"1" == thismajor(v"1.2.3")
@test v"1" == thismajor(v"1.2.3+post")
@test v"1" == thismajor(v"1.2.3+")

@test v"1.2.3" == nextpatch(v"1.2.3-")
@test v"1.2.3" == nextpatch(v"1.2.3-pre")
@test v"1.2.4" == nextpatch(v"1.2.3")
@test v"1.2.4" == nextpatch(v"1.2.3+post")
@test v"1.2.4" == nextpatch(v"1.2.3+")

@test v"1.2" == nextminor(v"1.2-")
@test v"1.2" == nextminor(v"1.2-pre")
@test v"1.3" == nextminor(v"1.2")
@test v"1.3" == nextminor(v"1.2+post")
@test v"1.3" == nextminor(v"1.2+")

@test v"1.3" == nextminor(v"1.2.3-")
@test v"1.3" == nextminor(v"1.2.3-pre")
@test v"1.3" == nextminor(v"1.2.3")
@test v"1.3" == nextminor(v"1.2.3+post")
@test v"1.3" == nextminor(v"1.2.3+")

@test v"1" == nextmajor(v"1-")
@test v"1" == nextmajor(v"1-pre")
@test v"2" == nextmajor(v"1")
@test v"2" == nextmajor(v"1+post")
@test v"2" == nextmajor(v"1+")

@test v"2" == nextmajor(v"1.2-")
@test v"2" == nextmajor(v"1.2-pre")
@test v"2" == nextmajor(v"1.2")
@test v"2" == nextmajor(v"1.2+post")
@test v"2" == nextmajor(v"1.2+")

@test v"2" == nextmajor(v"1.2.3-")
@test v"2" == nextmajor(v"1.2.3-pre")
@test v"2" == nextmajor(v"1.2.3")
@test v"2" == nextmajor(v"1.2.3+post")
@test v"2" == nextmajor(v"1.2.3+")

for major=0:3, minor=0:3, patch=0:3
	a = VersionNumber(major,minor,patch,("",),())
	b = VersionNumber(major,minor,patch,("pre",),())
	c = VersionNumber(major,minor,patch,(),())
	d = VersionNumber(major,minor,patch,(),("post",))
	e = VersionNumber(major,minor,patch,(),("",))
	@test a < b < c < d < e
	for x in [a,b,c,d,e]
		@test thispatch(x) == VersionNumber(major,minor,patch)
		@test thisminor(x) == VersionNumber(major,minor,0)
		@test thismajor(x) == VersionNumber(major,0,0)
		@test x < nextpatch(x) <= nextminor(x) <= nextmajor(x)
		@test x < thispatch(x) ? nextpatch(x) == thispatch(x) : thispatch(x) < nextpatch(x)
		@test x < thisminor(x) ? nextminor(x) == thisminor(x) : thisminor(x) < nextminor(x)
		@test x < thismajor(x) ? nextmajor(x) == thismajor(x) : thismajor(x) < nextmajor(x)
	end
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using FDM: jvp, j′vp

const _fdm = central_fdm(5, 1)

"""
    frule_test(f, (x, ẋ)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function for which the `frule` should be tested.
- `x`: input at which to evaluate `f` (should generally be set randomly).
- `ẋ`: differential w.r.t. `x` (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function frule_test(f, (x, ẋ); rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    return frule_test(f, ((x, ẋ),); rtol=rtol, atol=atol, fdm=fdm, kwargs...)
end

function frule_test(f, xẋs::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    xs, ẋs = collect(zip(xẋs...))
    Ω, dΩ_rule = ChainRules.frule(f, xs...)
    @test f(xs...) == Ω

    dΩ_ad, dΩ_fd = dΩ_rule(ẋs...), jvp(fdm, xs->f(xs...), (xs, ẋs))
    @test isapprox(dΩ_ad, dΩ_fd; rtol=rtol, atol=atol, kwargs...)
end

"""
    rrule_test(f, ȳ, (x, x̄)...; rtol=1e-9, atol=1e-9, fdm=central_fdm(5, 1), kwargs...)

# Arguments
- `f`: Function to which rule should be applied.
- `ȳ`: adjoint w.r.t. output of `f` (should generally be set randomly).
- `x`: input at which to evaluate `f` (should generally be set randomly).
- `x̄`: currently accumulated adjoint (should generally be set randomly).

All keyword arguments except for `fdm` are passed to `isapprox`.
"""
function rrule_test(f, ȳ, (x, x̄)::Tuple{Any, Any}; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    # Check correctness of evaluation.
    fx, dx = ChainRules.rrule(f, x)
    @test fx ≈ f(x)

    # Correctness testing via finite differencing.
    x̄_ad, x̄_fd = dx(ȳ), j′vp(fdm, f, ȳ, x)
    @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)

    # Assuming x̄_ad to be correct, check that other ChainRules mechanisms are correct.
    test_accumulation(x̄, dx, ȳ, x̄_ad)
    test_accumulation(Zero(), dx, ȳ, x̄_ad)
end

function _make_fdm_call(fdm, f, ȳ, xs, ignores)
    sig = Expr(:tuple)
    call = Expr(:call, f)
    newxs = Any[]
    arginds = Int[]
    i = 1
    for (x, ignore) in zip(xs, ignores)
        if ignore
            push!(call.args, x)
        else
            push!(call.args, Symbol(:x, i))
            push!(sig.args, Symbol(:x, i))
            push!(newxs, x)
            push!(arginds, i)
        end
        i += 1
    end
    fdexpr = :(j′vp($fdm, $sig -> $call, $ȳ, $(newxs...)))
    fd = eval(fdexpr)
    fd isa Tuple || (fd = (fd,))
    args = Any[nothing for _ in 1:length(xs)]
    for (dx, ind) in zip(fd, arginds)
        args[ind] = dx
    end
    return (args...,)
end

function rrule_test(f, ȳ, xx̄s::Tuple{Any, Any}...; rtol=1e-9, atol=1e-9, fdm=_fdm, kwargs...)
    # Check correctness of evaluation.
    xs, x̄s = collect(zip(xx̄s...))
    y, rules = rrule(f, xs...)
    @test f(xs...) == y

    # Correctness testing via finite differencing.
    x̄s_ad = map(rules) do rule
        rule isa DNERule ? DNE() : rule(ȳ)
    end
    x̄s_fd = _make_fdm_call(fdm, f, ȳ, xs, x̄s .== nothing)
    for (x̄_ad, x̄_fd) in zip(x̄s_ad, x̄s_fd)
        if x̄_fd === nothing
            # The way we've structured the above, this tests that the rule is a DNERule
            @test x̄_ad isa DNE
        else
            @test isapprox(x̄_ad, x̄_fd; rtol=rtol, atol=atol, kwargs...)
        end
    end

    # Assuming the above to be correct, check that other ChainRules mechanisms are correct.
    for (x̄, rule, x̄_ad) in zip(x̄s, rules, x̄s_ad)
        x̄ === nothing && continue
        test_accumulation(x̄, rule, ȳ, x̄_ad)
        test_accumulation(Zero(), rule, ȳ, x̄_ad)
    end
end

function Base.isapprox(ad::Wirtinger, fd; kwargs...)
    error("Finite differencing with Wirtinger rules not implemented")
end
function Base.isapprox(d_ad::Casted, d_fd; kwargs...)
    return all(isapprox.(extern(d_ad), d_fd; kwargs...))
end
function Base.isapprox(d_ad::DNE, d_fd; kwargs...)
    error("Tried to differentiate w.r.t. a DNE")
end
function Base.isapprox(d_ad::Thunk, d_fd; kwargs...)
    return isapprox(extern(d_ad), d_fd; kwargs...)
end

function test_accumulation(x̄, dx, ȳ, partial)
    @test all(extern(add(x̄, partial)) .≈ extern(x̄) .+ extern(partial))
    test_accumulate(x̄, dx, ȳ, partial)
    test_accumulate!(x̄, dx, ȳ, partial)
    test_store!(x̄, dx, ȳ, partial)
    return nothing
end

function test_accumulate(x̄::Zero, dx, ȳ, partial)
    @test extern(accumulate(x̄, dx, ȳ)) ≈ extern(partial)
    return nothing
end

function test_accumulate(x̄::Number, dx, ȳ, partial)
    @test extern(accumulate(x̄, dx, ȳ)) ≈ extern(x̄) + extern(partial)
    return nothing
end

function test_accumulate(x̄::AbstractArray, dx, ȳ, partial)
    x̄_old = copy(x̄)
    @test all(extern(accumulate(x̄, dx, ȳ)) .≈ (extern(x̄) .+ extern(partial)))
    @test x̄ == x̄_old
    return nothing
end

test_accumulate!(x̄::Zero, dx, ȳ, partial) = nothing

function test_accumulate!(x̄::Number, dx, ȳ, partial)
    @test accumulate!(x̄, dx, ȳ) ≈ accumulate(x̄, dx, ȳ)
    return nothing
end

function test_accumulate!(x̄::AbstractArray, dx, ȳ, partial)
    x̄_copy = copy(x̄)
    accumulate!(x̄_copy, dx, ȳ)
    @test extern(x̄_copy) ≈ (extern(x̄) .+ extern(partial))
    return nothing
end

test_store!(x̄::Zero, dx, ȳ, partial) = nothing
test_store!(x̄::Number, dx, ȳ, partial) = nothing

function test_store!(x̄::AbstractArray, dx, ȳ, partial)
    x̄_copy = copy(x̄)
    store!(x̄_copy, dx, ȳ)
    @test all(x̄_copy .≈ extern(partial))
    return nothing
end

ext: jl


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

# the 0mq clustermanager depends on package ZMQ. For testing purposes, at least
# make sure the code loads without it.
try
    using ZMQ
end

import Base: launch, manage, connect, kill

const BROKER_SUB_PORT = 8100
const BROKER_PUB_PORT = 8101

const SELF_INITIATED = 0
const REMOTE_INITIATED = 1

const PAYLOAD_MSG = "J"
const CONTROL_MSG = "Z"

const REQUEST_ACK = "R"
const ACK_MSG = "A"
const KILL_MSG = "K"

mutable struct ZMQCMan <: ClusterManager
    map_zmq_julia::Dict{Int, Tuple}
    c::Condition
    isfree::Bool
    ctx
    pub
    sub
    zid_self
    ZMQCMan() = new(Dict{Int, Tuple}(), Condition(), true)
end

const manager = ZMQCMan()

function lock_for_send()
    if manager.isfree == true
        manager.isfree = false
    else
        while manager.isfree == false
            wait(manager.c)
            if manager.isfree == true
                manager.isfree = false
                return
            end
        end
    end
end

function release_lock_for_send()
    manager.isfree = true
    notify(manager.c, all=true)
end

function init_node(zid=0)
    manager.ctx = Context(1)
    pub=Socket(manager.ctx, PUB)    # Outbound
    connect(pub, "tcp://127.0.0.1:$BROKER_SUB_PORT")

    sub=Socket(manager.ctx, SUB)    # In bound
    connect(sub, "tcp://127.0.0.1:$BROKER_PUB_PORT")
    ZMQ.set_subscribe(sub, string(zid))

    manager.pub = pub
    manager.sub = sub
    manager.zid_self = zid

    (pub, sub)
end

function send_data(zid, mtype, data)
    lock_for_send()
    ZMQ.send(manager.pub, Message(string(zid)), SNDMORE)
    ZMQ.send(manager.pub, Message(string(manager.zid_self)), SNDMORE)
    #println("Sending message of type $mtype to $zid")
    ZMQ.send(manager.pub, Message(mtype), SNDMORE)
    ZMQ.send(manager.pub, Message(data))
    release_lock_for_send()
end

function setup_connection(zid, initiated_by)
    try
        read_stream=BufferStream()
        write_stream=BufferStream()

        if initiated_by == REMOTE_INITIATED
            test_remote = false
        else
            test_remote = true
        end

        manager.map_zmq_julia[zid] = (read_stream, write_stream, test_remote)

        @schedule begin
            while true
                (r_s, w_s, do_test_remote) = manager.map_zmq_julia[zid]
                if do_test_remote
                    send_data(zid, CONTROL_MSG, REQUEST_ACK)
                    sleep(0.5)
                else
                    break
                end
            end
            (r_s, w_s, do_test_remote) = manager.map_zmq_julia[zid]

            while true
                data = readavailable(w_s)
                send_data(zid, PAYLOAD_MSG, data)
            end
        end
        (read_stream, write_stream)
    catch e
        Base.show_backtrace(STDOUT,catch_backtrace())
        println(e)
        rethrow(e)
    end
end

# BROKER
function start_broker()
    ctx=Context(1)
    xpub=Socket(ctx, XPUB)
    xsub=Socket(ctx, XSUB)

    ZMQ.bind(xsub, "tcp://127.0.0.1:$(BROKER_SUB_PORT)")
    ZMQ.bind(xpub, "tcp://127.0.0.1:$(BROKER_PUB_PORT)")

    ccall((:zmq_proxy, :libzmq), Cint,  (Ptr{Void}, Ptr{Void}, Ptr{Void}), xsub.data, xpub.data, C_NULL)
#    proxy(xsub, xpub)

    # control never comes here
    ZMQ.close(xpub)
    ZMQ.close(xsub)
    ZMQ.close(ctx)
end

function recv_data()
    try
        #println("On $(manager.zid_self) waiting to recv message")
        zid = parse(Int,String(ZMQ.recv(manager.sub)))
        assert(zid == manager.zid_self)

        from_zid = parse(Int,String(ZMQ.recv(manager.sub)))
        mtype = String(ZMQ.recv(manager.sub))

        #println("$zid received message of type $mtype from $from_zid")

        data = ZMQ.recv(manager.sub)
        if mtype == CONTROL_MSG
            cmsg = String(data)
            if cmsg == REQUEST_ACK
                #println("$from_zid REQUESTED_ACK from $zid")
                # send back a control_msg
                send_data(from_zid, CONTROL_MSG, ACK_MSG)
            elseif cmsg == ACK_MSG
                #println("$zid got ACK_MSG from $from_zid")
                (r_s, w_s, test_remote) = manager.map_zmq_julia[from_zid]
                manager.map_zmq_julia[from_zid] = (r_s, w_s, false)
            elseif cmsg == KILL_MSG
                exit(0)
            else
                error("Unknown control message : ", cmsg)
            end
            data = ""
        end

        (from_zid, data)
    catch e
        Base.show_backtrace(STDOUT,catch_backtrace())
        println(e)
        rethrow(e)
    end

end

# MASTER
function start_master(np)
    init_node()
    @schedule begin
        try
            while true
                (from_zid, data) = recv_data()

                #println("master recv data from $from_zid")

                (r_s, w_s, t_r) = manager.map_zmq_julia[from_zid]
                unsafe_write(r_s, pointer(data), length(data))
            end
        catch e
            Base.show_backtrace(STDOUT,catch_backtrace())
            println(e)
            rethrow(e)
        end
    end

    addprocs(manager; np=np)
end

function launch(manager::ZMQCMan, params::Dict, launched::Array, c::Condition)
    #println("launch $(params[:np])")
    for i in 1:params[:np]
        io, pobj = open(`$(params[:exename]) worker.jl $i $(Base.cluster_cookie())`, "r")

        wconfig = WorkerConfig()
        wconfig.userdata = Dict(:zid=>i, :io=>io)
        push!(launched, wconfig)
        notify(c)
    end
end

function connect(manager::ZMQCMan, pid::Int, config::WorkerConfig)
    #println("connect_m2w")
    if myid() == 1
        zid = get(config.userdata)[:zid]
        config.connect_at = zid # This will be useful in the worker-to-worker connection setup.

        print_worker_stdout(get(config.userdata)[:io], pid)
    else
        #println("connect_w2w")
        zid = get(config.connect_at)
        config.userdata = Dict{Symbol, Any}(:zid=>zid)
    end

    streams = setup_connection(zid, SELF_INITIATED)

    udata = get(config.userdata)
    udata[:streams] = streams

    streams
end

# WORKER
function start_worker(zid, cookie)
    #println("start_worker")
    Base.init_worker(cookie, ZMQCMan())
    init_node(zid)

    while true
        (from_zid, data) = recv_data()

        #println("worker recv data from $from_zid")

        streams = get(manager.map_zmq_julia, from_zid, nothing)
        if streams === nothing
            # First time..
            (r_s, w_s) = setup_connection(from_zid, REMOTE_INITIATED)
            Base.process_messages(r_s, w_s)
        else
            (r_s, w_s, t_r) = streams
        end

        unsafe_write(r_s, pointer(data), length(data))
    end
end

function manage(manager::ZMQCMan, id::Int, config::WorkerConfig, op)
    nothing
end

function kill(manager::ZMQCMan, pid::Int, config::WorkerConfig)
    send_data(get(config.userdata)[:zid], CONTROL_MSG, KILL_MSG)
    (r_s, w_s) = get(config.userdata)[:streams]
    close(r_s)
    close(w_s)

    # remove from our map
    delete!(manager.map_zmq_julia, get(config.userdata)[:zid])

    nothing
end

function print_worker_stdout(io, pid)
    @schedule while !eof(io)
        line = readline(io)
        println("\tFrom worker $(pid):\t$line")
    end
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: xorJane/julia

max_forks_repo_name: xorJane/julia


text:
"""
    benchmark_solving!(;
        valueSelectionArray::Union{T, Array{T, 1}},
        generator::AbstractModelGenerator,
        nbEpisodes::Int64=10,
        strategy::Type{DFSearch}=DFSearch,
        variableHeuristic=selectVariable,
        metricsFun,
        verbose::Bool=true
    ) where T <: ValueSelection

Used to benchmark the performance of some heuristics. It basically put all the LearnedHeuristics in test mode, which
means that the weights of the approximators are no more updated. Consequently, it's also possible to measure a good
aproximation of the time performances.

This function might evolve for a more precise one which will provide better analysis of the performances.
"""
function benchmark_solving(;
        valueSelectionArray::Union{T, Array{T, 1}},
        generator::AbstractModelGenerator,
        strategy::S=DFSearch(),
        variableHeuristic::AbstractVariableSelection=MinDomainVariableSelection(),
        out_solver::Bool=false,
        metricsFun=((;kwargs...) -> nothing),
        verbose::Bool=true,
        evaluator::Union{Nothing, AbstractEvaluator}
    ) where{T <: ValueSelection,S <: SearchStrategy}
    println("Benchmarking...")
    if isa(valueSelectionArray, T)
        valueSelectionArray = [valueSelectionArray]
    end

    for valueSelection in valueSelectionArray
        # give information to the learned heuristic
        if isa(valueSelection, LearnedHeuristic)
            if valueSelection.fitted_problem != typeof(generator)
                @warn "This learned heuristic was trained on a different problem type."
            end

            if valueSelection.fitted_strategy != typeof(strategy)
                @warn "This learned heuristic was trained with a different search strategy."
            end

            # make sure it is in testing mode
            testmode!(valueSelection, true)
        end
    end

    nbHeuristics = length(valueSelectionArray)
    nbInstances = evaluator.nbInstances

    if isnothing(evaluator)
        evaluator=SameInstancesEvaluator(valueSelectionArray, generator)
    end
    evaluate(evaluator, variableHeuristic, strategy)

    return evaluator.metrics
end


text:
using BinaryTraits
using BinaryTraits.Prefix: Is, Not

# See https://docs.julialang.org/en/v1/manual/interfaces

# -----------------------------------------------------------------------------
# Iteration interface
# -----------------------------------------------------------------------------
import Base: iterate
@trait Iterable
@implement Is{Iterable} by iterate(_)::Any
@implement Is{Iterable} by iterate(_, state::Any)::Any

# Example from https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1
struct Squares
    count::Int
end

Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)

# Let's assign the Squares type to Iterable
@assign Squares with Is{Iterable}
@check(Squares)
#=
julia> @check(Squares)
✅ Squares has implemented:
1. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹)::Any
2. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹, ::Any)::Any
=#

# -----------------------------------------------------------------------------
# Indexing interface
# -----------------------------------------------------------------------------
import Base: getindex, setindex!, firstindex, lastindex

@trait Indexable
@implement Is{Indexable} by getindex(_, i::Any)

@trait IndexableFromBeginning
@implement Is{IndexableFromBeginning} by firstindex(_)

@trait IndexableAtTheEnd
@implement Is{IndexableAtTheEnd} by lastindex(_)

# Make sure that `i` is untyped (i.e. `Any`) to adhere to the contract
function Base.getindex(S::Squares, i)
    1 <= i <= S.count || throw(BoundsError(S, i))
    return i*i
end

@assign Squares with Is{Indexable}
@check(Squares)
#=
julia> @check(Squares)
✅ Squares has implemented:
1. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹)::Any
2. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹, ::Any)::Any
3. BinaryTrait{Indexable}: Positive{Indexable} ⇢ getindex(🔹, ::Any)::Any
=#

# We want to have the traits for indexing from beginning and at the end
@assign Squares with Is{IndexableFromBeginning}, Is{IndexableAtTheEnd}
@check(Squares)
#=
julia> @check(Squares)
┌ Warning: Missing implementation
│   contract = BinaryTrait{IndexableFromBeginning}: Positive{IndexableFromBeginning} ⇢ firstindex(🔹)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits.jl/src/interface.jl:59
┌ Warning: Missing implementation
│   contract = BinaryTrait{IndexableAtTheEnd}: Positive{IndexableAtTheEnd} ⇢ lastindex(🔹)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits.jl/src/interface.jl:59
✅ Squares has implemented:
1. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹)::Any
2. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹, ::Any)::Any
3. BinaryTrait{Indexable}: Positive{Indexable} ⇢ getindex(🔹, ::Any)::Any
❌ Squares is missing these implementations:
1. BinaryTrait{IndexableFromBeginning}: Positive{IndexableFromBeginning} ⇢ firstindex(🔹)::Any
2. BinaryTrait{IndexableAtTheEnd}: Positive{IndexableAtTheEnd} ⇢ lastindex(🔹)::Any
=#

# Let's implement them now.
Base.firstindex(S::Squares) = 1
Base.lastindex(S::Squares) = length(S)
@check(Squares)
#=
julia> @check(Squares)
✅ Squares has implemented:
1. BinaryTrait{IndexableFromBeginning}: Positive{IndexableFromBeginning} ⇢ firstindex(🔹)::Any
2. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹)::Any
3. BinaryTrait{Iterable}: Positive{Iterable} ⇢ iterate(🔹, ::Any)::Any
4. BinaryTrait{Indexable}: Positive{Indexable} ⇢ getindex(🔹, ::Any)::Any
5. BinaryTrait{IndexableAtTheEnd}: Positive{IndexableAtTheEnd} ⇢ lastindex(🔹)::Any
=#


text:
import StaticArrays.SVector
import Colors: RGB, distinguishable_colors
import GeometryBasics: Point

export Layer, PenPlot

"""
Represents a layer of a plot.

A "layer" is a full pass of one pen. A layer has an associated color which
is used for previewing the plot in the SVG export, but the color is not used
by the plotting software.
"""
struct Layer
    name::String
    paths::MultiPath
    color::RGB
end

mutable struct LayerMaker
    num_layers::Int
    default_colors::Array{RGB,1}
    current_index::Int

    function LayerMaker(num_layers::Int)
        default_colors = distinguishable_colors(num_layers, lchoices = 0:40)
        new(num_layers, default_colors, 1)
    end
end

function makelayer(layer_maker::LayerMaker, layer::Layer)
    layer_maker.current_index += 1
    layer
end

function makelayer(layer_maker::LayerMaker, paths::MultiPath)
    name = "$(layer_maker.current_index)-layer"
    color = layer_maker.default_colors[layer_maker.current_index]
    layer_maker.current_index += 1
    Layer(name, paths, color)
end

function makelayer(layer_maker::LayerMaker, name_paths_pair::Pair{String,MultiPath})
    color = layer_maker.default_colors[layer_maker.current_index]
    layer_maker.current_index += 1
    Layer(name_paths_pair.first, name_paths_pair.second, color)
end

function makenotch(corner::Point, prev::Point, next::Point; frac=0.05)
    p1 = prev * frac + corner * (1-frac)
    p2 = corner
    p3 = next * frac + corner * (1-frac)
    Path([p1, p2, p3])
end

"""
A plot, consisting of one or more layers.

Constructed with:

    PenPlot(layers...)

Construct a pen plot consisting of the given layers.

Layers may be:
- [`Layer`](@ref) objects
- [`MultiPath`](@ref) objects
- `Pair{String, MultiPath}` objects, where the string becomes the name of
  the layer.
"""
struct PenPlot
    layers::Array{Layer}

    function PenPlot(layers...; addoutline=false)
        new_layers = []
        layer_maker = LayerMaker(length(layers) + (addoutline ? 1 : 0))

        new_layers = map(collect(layers)) do layerdata
            makelayer(layer_maker, layerdata)
        end
        plot = new(new_layers)

        if addoutline
            plx = extent(plot)

            upperright = Point(plx.lowerright[1], plx.upperleft[2])
            lowerleft = Point(plx.upperleft[1], plx.lowerright[2])

            outline = [
                makenotch(plx.upperleft, lowerleft, upperright),
                makenotch(upperright, plx.upperleft, plx.lowerright),
                makenotch(plx.lowerright, upperright, lowerleft),
                makenotch(lowerleft, plx.lowerright, plx.upperleft),
            ]
            push!(plot.layers, makelayer(layer_maker, outline))
        end

        plot
    end
end


text:
## variables (general)

function MOI.add_variable(o::Optimizer)
    allow_modification(o)
    vr = add_variable(o.mscip)
    var::Ptr{SCIP_VAR} = o.mscip.vars[vr][]
    register!(o, var, vr)
    return MOI.VariableIndex(vr.val)
end

MOI.add_variables(o::Optimizer, n) = [MOI.add_variable(o) for i=1:n]
MOI.get(o::Optimizer, ::MOI.NumberOfVariables) = length(o.mscip.vars)
MOI.get(o::Optimizer, ::MOI.ListOfVariableIndices) = [VI(k.val) for k in keys(o.mscip.vars)]
MOI.is_valid(o::Optimizer, vi::VI) = haskey(o.mscip.vars, VarRef(vi.value))

function MOI.get(o::Optimizer, ::MOI.VariableName, vi::VI)::String
    return GC.@preserve o unsafe_string(SCIPvarGetName(var(o, vi)))
end

function MOI.set(o::Optimizer, ::MOI.VariableName, vi::VI, name::String)
    @SC SCIPchgVarName(o, var(o, vi), name)
    return nothing
end

function MOI.delete(o::Optimizer, vi::VI)
    # Note: SCIP does not currently support deletion of variables that are still
    # present in some constraint. We don't want the overhead of keeping track of
    # the variable-in-constraint relation, so, to be conservative, we only allow
    # to delete a variable when there are no constraints in the model.
    if length(o.mscip.conss) > 0
        error("Can not delete variable while model contains constraints!")
    end

    allow_modification(o)
    haskey(o.binbounds, vi) && delete!(o.binbounds, vi)
    delete!(o.reference, var(o, vi))
    delete(o.mscip, VarRef(vi.value))
    return nothing
end

## variable types (binary, integer)

MOI.supports_constraint(o::Optimizer, ::Type{SVF}, ::Type{<:VAR_TYPES}) = true

scip_vartype(::Type{MOI.ZeroOne}) = SCIP_VARTYPE_BINARY
scip_vartype(::Type{MOI.Integer}) = SCIP_VARTYPE_INTEGER

function MOI.add_constraint(o::Optimizer, func::SVF, set::S) where {S <: VAR_TYPES}
    allow_modification(o)
    vi = func.variable
    v = var(o, vi)
    infeasible = Ref{SCIP_Bool}()
    @SC SCIPchgVarType(o, v, scip_vartype(S), infeasible)
    # TODO: warn if infeasible[] == TRUE?
    if S <: MOI.ZeroOne
        # Need to adjust bounds for SCIP, which fails with an error otherwise.
        # Check for conflicts with existing bounds first:
        lb, ub = SCIPvarGetLbOriginal(v), SCIPvarGetUbOriginal(v)
        if lb == -SCIPinfinity(o) && ub == SCIPinfinity(o)
            @debug "Implicitly setting bounds [0,1] for binary variable at $(vi.value)!"
            @SC SCIPchgVarLb(o, v, 0.0)
            @SC SCIPchgVarUb(o, v, 1.0)
        else
            # Store old bounds for later recovery.
            o.binbounds[vi] = MOI.Interval(lb, ub)
            if ub > 1.0
                @debug "Tightening upper bound $ub to 1 for binary variable at $(vi.value)!"
                @SC SCIPchgVarUb(o, v, 1.0)
            end
            if lb < 0.0
                @debug "Tightening lower bound $lb to 0 for binary variable at $(vi.value)!"
                @SC SCIPchgVarLb(o, v, 0.0)
            end
        end
    end
    # use var index for cons index of this type
    i = func.variable.value
    return register!(o, CI{SVF, S}(i))
end

function MOI.delete(o::Optimizer, ci::CI{SVF,S}) where {S <: VAR_TYPES}
    allow_modification(o)

    vi = VI(ci.value)
    v = var(o, vi)

    # Reset bounds from constraint if this was a binary, see below.
    reset_bounds = SCIPvarGetType(v) == SCIP_VARTYPE_BINARY

    # don't actually delete any SCIP constraint, just reset type
    infeasible = Ref{SCIP_Bool}()
    @SC SCIPchgVarType(o, var(o, vi), SCIP_VARTYPE_CONTINUOUS, infeasible)
    # TODO: warn if infeasible[] == TRUE?

    # Can only change bounds after chaging the var type.
    if reset_bounds
        bounds = o.binbounds[vi]
        @SC SCIPchgVarLb(o, v, bounds.lower)
        @SC SCIPchgVarUb(o, v, bounds.upper)
    end

    # but do delete the constraint reference
    delete!(o.constypes[SVF, S], ConsRef(ci.value))

    return nothing
end

## variable bounds

MOI.supports_constraint(o::Optimizer, ::Type{SVF}, ::Type{<:BOUNDS}) = true

function MOI.add_constraint(o::Optimizer, func::SVF, set::S) where S <: BOUNDS
    allow_modification(o)
    vi = func.variable
    v = var(o, vi)
    inf = SCIPinfinity(o)

    newlb, newub = bounds(set)
    newlb = newlb == nothing ? -inf : newlb
    newub = newub == nothing ?  inf : newub

    # Check for existing bounds first.
    oldlb, oldub = SCIPvarGetLbOriginal(v), SCIPvarGetUbOriginal(v)
    if (oldlb != -inf && newlb != -inf || oldub != inf && newub != inf)
        # special case: implicit bounds for binaries
        if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
            # Store new bounds
            o.binbounds[vi] = MOI.Interval(newlb, newub)

            if newlb < 0.0
                @debug "Ignoring relaxed lower bound $newlb for binary variable at $(vi.value)!"
                newlb = oldlb
            end

            if newub > 1.0
                @debug "Ignoring relaxed upper bounds $newub for binary variable at $(vi.value)!"
                newub = oldub
            end

            if newlb != oldlb || newub != oldub
                @debug "Overwriting existing bounds [0.0,1.0] with [$newlb,$newub] for binary variable at $(vi.value)!"
            end
        else # general case (non-binary variable)
            error("Already have bounds [$oldlb,$oldub] for variable at $(vi.value)!")
        end
    end

    if newlb != -inf
        @SC SCIPchgVarLb(o, v, newlb)
    end

    if newub != inf
        @SC SCIPchgVarUb(o, v, newub)
    end

    # use var index for cons index of this type
    i = func.variable.value
    return register!(o, CI{SVF, S}(i))
end

function reset_bounds(o::Optimizer, v, lb, ub, ci::CI{SVF, S}) where S <: Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}}
    @SC SCIPchgVarLb(o, v, lb)
    @SC SCIPchgVarUb(o, v, ub)
end

function reset_bounds(o::Optimizer, v, lb, ub, ci::CI{SVF, MOI.GreaterThan{Float64}})
    @SC SCIPchgVarLb(o, v, lb)
end

function reset_bounds(o::Optimizer, v, lb, ub, ci::CI{SVF, MOI.LessThan{Float64}})
    @SC SCIPchgVarUb(o, v, ub)
end

function MOI.delete(o::Optimizer, ci::CI{SVF,S}) where S <: BOUNDS
    allow_modification(o)

    # Don't actually delete any SCIP constraint, just reset bounds
    vi = VI(ci.value)
    v = var(o, vi)
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        reset_bounds(o, v, 0.0, 1.0, ci)
    else
        inf = SCIPinfinity(o)
        reset_bounds(o, v, -inf, inf, ci)
    end

    # but do delete the constraint reference
    haskey(o.binbounds, vi) && delete!(o.binbounds, vi)
    delete!(o.constypes[SVF, S], ConsRef(ci.value))

    return nothing
end

function MOI.set(o::SCIP.Optimizer, ::MOI.ConstraintSet, ci::CI{SVF,S}, set::S) where {S <: BOUNDS}
    allow_modification(o)
    v = var(o, VI(ci.value)) # cons index is actually var index
    lb, ub = bounds(set)
    lb = lb == nothing ? -SCIPinfinity(o) : lb
    ub = ub == nothing ?  SCIPinfinity(o) : ub
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        o.binbounds[vi] = MOI.Interval(lb, ub)
        lb = max(lb, 0.0)
        ub = min(lb, 1.0)
    end
    @SC SCIPchgVarLb(o, v, lb)
    @SC SCIPchgVarUb(o, v, ub)
    return nothing
end

# TODO: is actually wrong for unbounded variables?
function MOI.is_valid(o::Optimizer, ci::CI{SVF,<:BOUNDS})
    return haskey(o.mscip.vars, VarRef(ci.value))
end

function MOI.get(o::Optimizer, ::MOI.ConstraintFunction, ci::CI{SVF, S}) where S <: BOUNDS
    return SVF(VI(ci.value))
end

function MOI.get(o::Optimizer, ::MOI.ConstraintSet, ci::CI{SVF, S}) where S <: BOUNDS
    vi = VI(ci.value)
    v = var(o, vi)
    lb, ub = SCIPvarGetLbOriginal(v), SCIPvarGetUbOriginal(v)
    if SCIPvarGetType(v) == SCIP_VARTYPE_BINARY
        bounds = o.binbounds[vi]
        lb, ub = bounds.lower, bounds.upper
    end
    return from_bounds(S, lb, ub)
end

# (partial) warm starts

MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{VI}) = true

function MOI.get(o::Optimizer, ::MOI.VariablePrimalStart, vi::VI)
    return if haskey(o.start, vi)
        o.start[vi]
    else
        nothing
    end
end

function MOI.set(o::Optimizer, ::MOI.VariablePrimalStart, vi::VI, value::Float64)
    o.start[vi] = value
    return
end

function MOI.set(o::Optimizer, ::MOI.VariablePrimalStart, vi::VI, value::Nothing)
    if haskey(o.start, vi)
        delete!(o.start, vi)
    end
    return
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

"""
    SemidiscretizationEulerAcoustics(semi_acoustics::SemiAcoustics, semi_euler::SemiEuler;
                                     source_region=x->true, weights=x->1.0)

!!! warning "Experimental code"
    This semidiscretization is experimental and may change in any future release.

Construct a semidiscretization of the acoustic perturbation equations that is coupled with
the compressible Euler equations via source terms for the perturbed velocity. Both
semidiscretizations have to use the same mesh and solvers with a shared basis. The coupling region
is described by a function `source_region` that maps the coordinates of a single node to `true` or
`false` depending on whether the point lies within the coupling region or not. A weighting function
`weights` that maps coordinates to weights is applied to the acoustic source terms.
Note that this semidiscretization should be used in conjunction with
[`EulerAcousticsCouplingCallback`](@ref) and only works in two dimensions.
"""
struct SemidiscretizationEulerAcoustics{SemiAcoustics, SemiEuler, Cache} <: AbstractSemidiscretization
  semi_acoustics::SemiAcoustics
  semi_euler::SemiEuler
  performance_counter::PerformanceCounter
  cache::Cache

  function SemidiscretizationEulerAcoustics{SemiAcoustics, SemiEuler, Cache}(
      semi_acoustics, semi_euler, cache) where {SemiAcoustics, SemiEuler, Cache}

    # Currently both semidiscretizations need to use a shared mesh
    @assert semi_acoustics.mesh == semi_euler.mesh

    # Check if both solvers use the same polynomial basis
    @assert semi_acoustics.solver.basis == semi_euler.solver.basis

    performance_counter = PerformanceCounter()
    new(semi_acoustics, semi_euler, performance_counter, cache)
  end
end

function SemidiscretizationEulerAcoustics(semi_acoustics::SemiAcoustics, semi_euler::SemiEuler;
                                          source_region=x->true, weights=x->1.0) where
    {Mesh, SemiAcoustics<:SemidiscretizationHyperbolic{Mesh, <:AbstractAcousticPerturbationEquations},
     SemiEuler<:SemidiscretizationHyperbolic{Mesh, <:AbstractCompressibleEulerEquations}}

  cache = create_cache(SemidiscretizationEulerAcoustics, source_region, weights,
                       mesh_equations_solver_cache(semi_acoustics)...)

  return SemidiscretizationEulerAcoustics{typeof(semi_acoustics), typeof(semi_euler), typeof(cache)}(
    semi_acoustics, semi_euler, cache)
end

function create_cache(::Type{SemidiscretizationEulerAcoustics}, source_region, weights,
                      mesh, equations::AcousticPerturbationEquations2D, dg::DGSEM, cache)

  coupled_element_ids = get_coupled_element_ids(source_region, equations, dg, cache)

  acoustic_source_terms = zeros(eltype(cache.elements), (ndims(equations), nnodes(dg), nnodes(dg),
                                                         length(coupled_element_ids)))

  acoustic_source_weights = precompute_weights(source_region, weights, coupled_element_ids,
                                               equations, dg, cache)

  return (; acoustic_source_terms, acoustic_source_weights, coupled_element_ids)
end

function get_coupled_element_ids(source_region, equations, dg::DGSEM, cache)
  coupled_element_ids = Vector{Int}(undef, 0)

  for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      x = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, element)
      if source_region(x)
        push!(coupled_element_ids, element)
        break
      end
    end
  end

  return coupled_element_ids
end

function precompute_weights(source_region, weights, coupled_element_ids, equations, dg::DGSEM, cache)
  acoustic_source_weights = zeros(eltype(cache.elements),
                                  (nnodes(dg), nnodes(dg), length(coupled_element_ids)))

  @threaded for k in 1:length(coupled_element_ids)
    element = coupled_element_ids[k]
    for j in eachnode(dg), i in eachnode(dg)
      x = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, element)
      acoustic_source_weights[i, j, k] = source_region(x) ? weights(x) : zero(weights(x))
    end
  end

  return acoustic_source_weights
end

function Base.show(io::IO, semi::SemidiscretizationEulerAcoustics)
  @nospecialize semi # reduce precompilation time

  print(io, "SemidiscretizationEulerAcoustics(")
  print(io,       semi.semi_acoustics)
  print(io, ", ", semi.semi_euler)
  print(io, ", cache(")
  for (idx, key) in enumerate(keys(semi.cache))
    idx > 1 && print(io, " ")
    print(io, key)
  end
  print(io, "))")
end

function Base.show(io::IO, mime::MIME"text/plain", semi::SemidiscretizationEulerAcoustics)
  @nospecialize semi # reduce precompilation time

  if get(io, :compact, false)
    show(io, semi)
  else
    summary_header(io, "SemidiscretizationEulerAcoustics")
    summary_line(io, "semidiscretization Euler", semi.semi_euler |> typeof |> nameof)
    show(increment_indent(io), mime, semi.semi_euler)
    summary_line(io, "semidiscretization acoustics", semi.semi_acoustics |> typeof |> nameof)
    show(increment_indent(io), mime, semi.semi_acoustics)
    summary_footer(io)
  end
end

# The acoustics semidiscretization is the main semidiscretization.
@inline function mesh_equations_solver_cache(semi::SemidiscretizationEulerAcoustics)
  return mesh_equations_solver_cache(semi.semi_acoustics)
end

@inline Base.ndims(semi::SemidiscretizationEulerAcoustics) = ndims(semi.semi_acoustics)
@inline Base.real(semi::SemidiscretizationEulerAcoustics) = real(semi.semi_acoustics)

# Computes the coefficients of the initial condition
@inline function compute_coefficients(t, semi::SemidiscretizationEulerAcoustics)
  compute_coefficients(t, semi.semi_acoustics)
end

@inline function compute_coefficients!(u_ode, t, semi::SemidiscretizationEulerAcoustics)
  compute_coefficients!(u_ode, t, semi.semi_acoustics)
end

@inline function calc_error_norms(func, u, t, analyzer, semi::SemidiscretizationEulerAcoustics,
                                  cache_analysis)
  calc_error_norms(func, u, t, analyzer, semi.semi_acoustics, cache_analysis)
end

function rhs!(du_ode, u_ode, semi::SemidiscretizationEulerAcoustics, t)
  @unpack semi_acoustics, cache = semi
  @unpack acoustic_source_terms, acoustic_source_weights, coupled_element_ids = cache

  du_acoustics = wrap_array(du_ode, semi_acoustics)

  time_start = time_ns()

  @trixi_timeit timer() "acoustics rhs!" rhs!(du_ode, u_ode, semi_acoustics, t)

  @trixi_timeit timer() "add acoustic source terms" add_acoustic_source_terms!(
    du_acoustics, acoustic_source_terms, acoustic_source_weights, coupled_element_ids,
    mesh_equations_solver_cache(semi_acoustics)...)

  runtime = time_ns() - time_start
  put!(semi.performance_counter, runtime)

  return nothing
end

function add_acoustic_source_terms!(du_acoustics, acoustic_source_terms, source_weights,
                                    coupled_element_ids, mesh::TreeMesh{2}, equations, dg::DGSEM,
                                    cache)

  @threaded for k in 1:length(coupled_element_ids)
    element = coupled_element_ids[k]

    for j in eachnode(dg), i in eachnode(dg)
      du_acoustics[1, i, j, element] += source_weights[i, j, k] * acoustic_source_terms[1, i, j, k]
      du_acoustics[2, i, j, element] += source_weights[i, j, k] * acoustic_source_terms[2, i, j, k]
    end
  end

  return nothing
end

end # @muladd


text:
"""
    $(SIGNATURES)

An in-memory DistributedFactorGraph based on LightGraphs.jl with parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
mutable struct LightDFG{T <: AbstractParams, V <: AbstractDFGVariable, F <:AbstractDFGFactor} <: AbstractDFG
    g::FactorGraph{Int, V, F}
    description::String
    userId::String
    robotId::String
    sessionId::String
    userData::Dict{Symbol, String}
    robotData::Dict{Symbol, String}
    sessionData::Dict{Symbol, String}
    addHistory::Vector{Symbol} #TODO: Discuss more - is this an audit trail?
    solverParams::T # Solver parameters
end

"""
    $(SIGNATURES)

Create an in-memory LightDFG with the following parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
function LightDFG{T,V,F}(g::FactorGraph{Int,V,F}=FactorGraph{Int,V,F}();
                           description::String="LightGraphs.jl implementation",
                           userId::String="User ID",
                           robotId::String="Robot ID",
                           sessionId::String="Session ID",
                           userData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           robotData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           sessionData::Dict{Symbol, String} = Dict{Symbol, String}(),
                           params::T=NoSolverParams()) where {T <: AbstractParams, V <:AbstractDFGVariable, F<:AbstractDFGFactor}

    LightDFG{T,V,F}(g, description, userId, robotId, sessionId, userData, robotData, sessionData, Symbol[], params)
end

# LightDFG{T}(; kwargs...) where T <: AbstractParams = LightDFG{T,DFGVariable,DFGFactor}(;kwargs...)
"""
    $(SIGNATURES)

Create an in-memory LightDFG with the following parameters:
- T: Solver parameters (defaults to `NoSolverParams()`)
- V: Variable type
- F: Factor type
"""
LightDFG{T}(g::FactorGraph{Int,DFGVariable,DFGFactor}=FactorGraph{Int,DFGVariable,DFGFactor}(); kwargs...) where T <: AbstractParams = LightDFG{T,DFGVariable,DFGFactor}(g; kwargs...)

Base.propertynames(x::LightDFG, private::Bool=false) =
    (:g, :description, :userId, :robotId, :sessionId, :nodeCounter, :labelDict, :addHistory, :solverParams)
        # (private ? fieldnames(typeof(x)) : ())...)

Base.getproperty(x::LightDFG,f::Symbol) = begin
    if f == :nodeCounter
        @error "Field nodeCounter depreciated. returning number of nodes"
        nv(x.g)
    elseif f == :labelDict
        @error "Field labelDict depreciated. Consider using exists(dfg,label) or getLabelDict(dfg) instead. Returning internals copy"
        #TODO: https://github.com/JuliaRobotics/DistributedFactorGraphs.jl/issues/111
        copy(x.g.labels.sym_int)
    else
        getfield(x,f)
    end
end


text:
"""
    MarkedTimeSeries <: TS

A composite type to store time data associated with marks.
"""
struct MarkedTimeSeries <: TS
    times::AbstractVector{<:SupportedTimeTypes}
    marks::AbstractVector{<:Real}

    function MarkedTimeSeries(times, marks)
        length(times) == length(marks) || error("times and marks must be the same length")
        times = vec(times)
        marks = vec(marks)
        return new(times, marks)
    end
end

"""
    MarkedTimeSeries(ts::TS)::MarkedTimeSeries

A method to see a point process as marked with zero marks.
"""
function MarkedTimeSeries(ts::TS)::MarkedTimeSeries
    if ts isa TimeSeries
        times = ts.times
        marks = fill(0, size(times))
        return MarkedTimeSeries(times, marks)
    end

    return ts
end

"""
    copy_above(mts::MarkedTimeSeries, mag::Real)::MarkedTimeSeries

Create a copy of mts keeping only events above mag.
"""
function copy_above(mts::MarkedTimeSeries, mag::Real)::MarkedTimeSeries
    idx = mts.marks .> mag
    return MarkedTimeSeries(mts.times[idx], mts.marks[idx])
end


text:
"""
`ExplicitRKTableau`

Holds a tableau which defines an explicit Runge-Kutta method.
"""
mutable struct ExplicitRKTableau{MType<:AbstractMatrix,VType<:AbstractVector,fsal} <: ODERKTableau
  A::MType
  c::VType
  α::VType
  αEEst::VType
  stages::Int
  order::Int
  adaptiveorder::Int #The lower order of the pair. Only used for adaptivity.
end
ExplicitRKTableau(A::MType,c::VType,α::VType,order;
                  adaptiveorder=0,αEEst=VType(),fsal=false) where {MType,VType} =
                  ExplicitRKTableau{MType,VType,fsal}(
                  A,c,α,αEEst,length(α),order,adaptiveorder)

"""
`ImplicitRKTableau`

Holds a tableau which defines an implicit Runge-Kutta method.
"""
mutable struct ImplicitRKTableau{MType<:AbstractMatrix,VType<:AbstractVector} <: ODERKTableau
  A::MType
  c::VType
  α::VType
  αEEst::VType
  stages::Int
  order::Int
  adaptiveorder::Int #The lower order of the pair. Only used for adaptivity.
end
ImplicitRKTableau(A::MType,c::VType,α::VType,order;
                  adaptiveorder=0,αEEst=VType()) where {MType,VType} =
                  ImplicitRKTableau{MType,VType}(
                  A,c,α,αEEst,length(α),order,adaptiveorder)


text:
function generate_multi_index(d::S,k::S) where {S <: Integer} # Recursive function

  if d < 1
    error("d must be positive")
  end

  if k < 0
    error("k must be non-negative")
  end

  kprime = k+1

  if d == 1
    mi = zeros(S,kprime)
    for i = 1:k+1
      if i <= kprime
        mi[i] = i-1
      end
    end
    return mi
  else
    mi_base = generate_multi_index(d-1,k)
    N = size(mi_base,1)
    mi = zeros(S,N*kprime,d)
    pos = 0
    @inbounds @views for j = 1:N
      for i = 1:k+1
        if prod(mi_base[j,:].+1)*i <= kprime
          pos += 1
          mi[pos,1:d-1] .= mi_base[j,:]
          mi[pos,d] = i-1
        end
      end
    end
    return mi[1:pos,:]
  end
end

function generate_multi_index(d::S,k::S,n::S) where {S <: Integer} # Recursive function

  if d < 1
    error("d must be positive")
  end

  if k < 0
    error("k must be non-negative")
  end

  if n < 1
    error("n must be positive")
  end

  if !isodd(n)
    error("n must be odd")
  end

  if n >= 2*k+1
    return generate_multi_index(d,k)
  end

  kprime = k+1
  #ki = Int((n-1)/2) + 1
  ki = div(n-1,2) + 1

  if d == 1
    mi = zeros(S,ki)
    for i = 1:ki
      if i <= kprime
        mi[i] = i-1
      end
    end
    return mi
  else
    mi_base = generate_multi_index(d-1,k,n)
    N = size(mi_base,1)
    mi = zeros(S,N*ki,d)
    pos = 0
    @inbounds @views for j = 1:N
      for i = 1:ki
        if prod(mi_base[j,:].+1)*i <= kprime
          pos += 1
          mi[pos,1:d-1] .= mi_base[j,:]
          mi[pos,d] = i-1
        end
      end
    end
    return mi[1:pos,:]
  end

end

function generate_multi_index(d::S,k::S,n::Array{S,1}) where {S <: Integer} # Recursive function

  if d < 1
    error("d must be positive")
  end

  if k < 0
    error("k must be non-negative")
  end

  if length(n) != d
    error("node number for each dimension must be given")
  end

  for i = 1:d
    if n[i] < 1
      error("n must be positive in all dimensions")
    end
    if !isodd(n[i])
      error("each element in n must be odd")
    end
  end

  kprime = k + 1
  #ki = Int((n[end]-1)/2) + 1
  ki = div(n[end]-1,2) + 1

  if d == 1
    mi = zeros(S,ki)
    for i = 1:ki
      if i <= kprime
        mi[i] = i-1
      end
    end
    return mi
  else
    mi_base = generate_multi_index(d-1,k,n[1:d-1])
    N = size(mi_base,1)
    mi = zeros(S,N*ki,d)
    pos = 0
    @inbounds @views for j = 1:N
      for i = 1:ki
        if prod(mi_base[j,:].+1)*i <= kprime
          pos += 1
          mi[pos,1:d-1] .= mi_base[j,:]
          mi[pos,d] = i-1
        end
      end
    end
    return mi[1:pos,:]
  end
end


text:
using Makie
using AbstractPlotting: textslider

λ = 0.1
g(x,t) = exp(-λ*π^2*t)*sin(π*x)

function main()
    t_slider, t_obs = textslider(-1:0.2:1,
                                 "t", start = 0)
    x = 0:0.001:1
    plotting = lift(t_obs) do t
        Point2f0.(x, g.(x, t))
    end

    s1 = lines(plotting)
    RecordEvents(hbox(s1, t_slider, parent = Scene()),
        "output",
    )
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module BlendModes

export BlendMode,
       BlendNormal,
       BlendMultiply,
       BlendScreen,
       BlendOverlay,
       BlendDarken,
       BlendLighten,
       BlendColorDodge,
       BlendColorBurn,
       BlendHardLight,
       BlendSoftLight,
       BlendDifference,
       BlendExclusion,
       BlendHue,
       BlendSaturation,
       BlendColor,
       BlendLuminosity

"""
    BlendMode{mode}

A type used for specifying the mixing mode of opaque colors. The `mode` should
be a symbol.
"""
struct BlendMode{mode} end

"""
    BlendNormal

The destination color is always the source color.
```
    Cdest = Csrc
```
"""
const BlendNormal = BlendMode{:normal}()

"""
    BlendMultiply

The source color is multiplied by the backdrop color.
```
    Cdest = Cb × Csrc
```
"""
const BlendMultiply = BlendMode{:multiply}()

"""
    BlendScreen

The complementary colors of the source and backdrop colors are multiplied, and
the destination color is the complementary color of the multiplicated color.
```
    Cdest = 1 - ((1 - Cb) × (1 - Csrc))
```
"""
const BlendScreen = BlendMode{:screen}()

"""
    BlendOverlay

This mode uses the [`multiply`](@ref BlendMultiply) or [`screen`](@ref
BlendScreen) mode, depending on the backdrop color. The overlay mode is the
inverse of the [`hard-light`](@ref BlendHardLight) mode.
```
    if Cb <= 0.5
        Cdest = Csrc × 2Cb
    else
        Cdest = 1 - ((1 - Csrc) × (1 - 2Cb))
    end
```
"""
const BlendOverlay = BlendMode{:overlay}()

"""
    BlendDarken

The darker values of the backdrop color and the source color are selected.
```
    Cdest = min(Cb, Csrc)
```
"""
const BlendDarken = BlendMode{:darken}()

"""
    BlendLighten

The lighter values of the backdrop color and the source color are selected.
```
    Cdest = max(Cb, Csrc)
```
"""
const BlendLighten = BlendMode{:lighten}()

"""
    BlendColorDodge

The destination color is the result of dividing the backdrop color by the
complementary color of the source color.
```
    if Cb == 0
        Cdest = 0
    elseif Csrc == 1
        Cdest = 1
    else
        Cdest = min(1, Cb / (1 - Csrc))
    end
```
"""
const BlendColorDodge = BlendMode{Symbol("color-dodge")}()

"""
    BlendColorBurn

The destination color is the result of dividing the complementary color of the
backdrop color by the source color.

ColorBlendModes uses the definition of W3C drafts as shown below. Note that
there is a variant, which returns `0` when `Cb == 1` and `Csrc == 0`.
```
    if Cb == 1
        Cdest = 1
    elseif Csrc == 0
        Cdest = 0
    else
        Cdest = 1 - min(1, (1 - Cb) / Csrc)
    end
```
"""
const BlendColorBurn = BlendMode{Symbol("color-burn")}()

"""
    BlendHardLight

This mode uses the [`multiply`](@ref BlendMultiply) or [`screen`](@ref
BlendScreen) mode, depending on the source color. The overlay mode is the
inverse of the [`overlay`](@ref BlendOverlay) mode.
```
    if Csrc <= 0.5
        Cdest = Cb × 2Csrc
    else
        Cdest = 1 - ((1 - Cb) × (1 - 2Csrc))
    end
```
"""
const BlendHardLight = BlendMode{Symbol("hard-light")}()

"""
    BlendSoftLight

The result is similar to the [`hard-light`](@ref BlendHardLight) mode, but
milder. This mode is also related to the [`overlay`](@ref BlendOverlay) mode.

ColorBlendModes uses the definition of W3C drafts as shown below. Note that
there are different definitions of the `soft-light`.
```
    if Csrc <= 0.5
        Cdest = Cb - (1 - 2Csrc) × Cb × (1 - Cb)
    else
        if Cb <= 0.25
            D = ((4Cb - 3) × Cb + 1) × 4Cb
        else
            D = sqrt(Cb)
        end
        Cdest = Cb + (2Csrc - 1) × (D - Cb)
    end
```
"""
const BlendSoftLight = BlendMode{Symbol("soft-light")}()

"""
    BlendDifference

The destination values are the subtraction of the darker values from the lighter
values of the backdrop and source colors.
```
    Cdest = abs(Csrc - Cb)
```
"""
const BlendDifference = BlendMode{:difference}()

"""
    BlendExclusion

The result is similar to the [`difference`](@ref BlendDifference) mode, but
milder.
```
    Cdest = Cb + Csrc - 2Cb × Csrc
```
"""
const BlendExclusion = BlendMode{:exclusion}()

"""
    BlendHue

The result is a color with the hue of the source color and the saturation and
luminosity of the backdrop color.
"""
const BlendHue= BlendMode{:hue}()

"""
    BlendSaturation

The result is a color with the saturation of the source color and the hue and
luminosity of the backdrop color.
"""
const BlendSaturation = BlendMode{:saturation}()

"""
    BlendColor

The result is a color with the hue and saturation of the source color and the
luminosity of the backdrop color.
"""
const BlendColor = BlendMode{:color}()

"""
    BlendLuminosity

The result is a color with the luminosity of the source color and the hue and
saturation of the backdrop color.
"""
const BlendLuminosity = BlendMode{:luminosity}()

end # module


text:
abstract type StochasticIntegratorRK{dType, tType, D, M, S} <: SDEIntegrator{dType, tType} end
abstract type StochasticIntegratorPRK{dType, tType, D, M, S} <: PSDEIntegrator{dType, tType} end

AbstractStochasticIntegratorRK{DT,TT,D,M,S} = Union{StochasticIntegratorRK{DT,TT,D,M,S}, StochasticIntegratorPRK{DT,TT,D,M,S}}

@inline parameters(integrator::AbstractStochasticIntegratorRK) = integrator.params
@inline equation(integrator::AbstractStochasticIntegratorRK, i::Symbol) = integrator.params.equs[i]
@inline equations(integrator::AbstractStochasticIntegratorRK) = integrator.params.equs
@inline timestep(integrator::AbstractStochasticIntegratorRK) = integrator.params.Δt
@inline tableau(integrator::AbstractStochasticIntegratorRK)  = integrator.params.tab

@inline nstages(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = S
@inline noisedims(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = M
@inline Base.ndims(::AbstractStochasticIntegratorRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = D
@inline Base.eltype(::AbstractStochasticIntegratorRK{DT}) where {DT} = DT

@inline eachstage(integrator::AbstractStochasticIntegratorRK) = 1:nstages(integrator)
@inline eachnoise(integrator::AbstractStochasticIntegratorRK) = 1:noisedims(integrator)
@inline eachdim(integrator::AbstractStochasticIntegratorRK) = 1:ndims(integrator)


text:
#Example of CellListMap application.

# Function to be evalulated for each pair: push pair if d<cutoff
function push_pair!(i,j,d2,pairs,cutoff)
  d = sqrt(d2)
  if d < cutoff
    push!(pairs,(i,j,d))
  end
  return pairs
end

# Reduction function
function reduce_pairs(pairs,pairs_threaded)
  pairs = pairs_threaded[1]
  for i in 2:nthreads()
    append!(pairs,pairs_threaded[i])
  end
  return pairs
end

# Initialize output
pairs = Tuple{Int,Int,Float64}[]

# Run pairwise computation
pairs = map_pairwise!(
  (x,y,i,j,d2,pairs) -> push_pair!(i,j,d2,pairs,cutoff),
  pairs,x,box,lc,
  reduce=reduce_pairs,
  parallel=parallel
)


text:
# Test cross-tabulation

# simple ClusteringResult subtype for testing
struct SimpleCluRes{T} <: ClusteringResult
    assignments::Vector{T}
    counts::Vector{Int}

    SimpleCluRes(assignments::AbstractVector{T}) where {T<:Integer} =
        new{T}(assignments,
               isempty(assignments) ? Vector{Int}() :
               counts(assignments, (1:maximum(assignments))))
end

@testset "counts() (contingency matrix)" begin

# StatsBase's counts() doesn't allow empty inputs
@test_throws ArgumentError counts(Int[], Int[])
# Clustering's counts()
@test counts(SimpleCluRes(Int[]), Int[]) == Matrix{Int}(undef, 0, 0)
@test_throws DimensionMismatch counts(SimpleCluRes([1]), Int[])
@test_throws DimensionMismatch counts(SimpleCluRes([1]), [2, 1])
#@test_throws ArgumentError counts([1, 1], [0, 1]) # doesn't throw as StatsBase.counts() is called
#@test_throws ArgumentError counts([1, -1], [2, 1]) # doesn't throw as StatsBase.counts() is called
@test_throws ArgumentError counts(SimpleCluRes([1, 1]), [0, 1])
@test_throws ArgumentError counts(SimpleCluRes([1, -1]), [2, 1])

# supports different Integer types
@test counts(SimpleCluRes(Int32[1, 2]), Int16[1, 2]) == Matrix{Int}(I, 2, 2)

@test counts([1, 2, 3], [1, 2, 3]) == Matrix{Int}(I, 3, 3)
@test counts([2, 3, 1], [1, 2, 3]) == [0 0 1; 1 0 0; 0 1 0]

# 3rd cluster in A missing
@test counts([2, 4, 1], [1, 2, 3]) == [0 0 1; 1 0 0; 0 0 0; 0 1 0]
@test counts(SimpleCluRes([2, 4, 1]), [1, 2, 3]) == [0 0 1; 1 0 0; 0 0 0; 0 1 0]

# 1st cluster in B missing (StatsBase.counts() and Clustering.counts() give different results)
@test counts([2, 3, 1], [2, 2, 3]) == [0 1; 1 0; 1 0]
@test counts(SimpleCluRes([2, 3, 1]), [2, 2, 3]) == [0 0 1; 0 1 0; 0 1 0]

@test counts([2, 2, 1], [1, 1, 1]) == reshape([1; 2], (2, 1))
@test counts([1, 1, 1], [2, 2, 1]) == [1 2]

@testset "with ClusteringResult objects" begin
    Random.seed!(34568)
    X = rand(3, 25)

    clu3 = kmeans(X, 3)
    clu5 = kmeans(X, 5)

    Y = rand(3, 20)
    cluY = kmeans(Y, 3)

    @test_throws DimensionMismatch counts(clu3, cluY)
    @test counts(clu3, clu5) == counts(clu5, clu3)'
    @test size(counts(clu3, clu5)) == (3, 5)
    @test sum(counts(clu3, clu5)) == 25
    @test counts(clu3, clu5) == counts(clu3, assignments(clu5))
    @test counts(clu3, clu5) == counts(assignments(clu3), clu5)
    @test counts(clu3, clu5) == counts(assignments(clu3), assignments(clu5))
end

end


text:
workspace()
using HoTRG_lemon.SimulatorModule
using HoTRG_lemon.SpinModule
using HoTRG_lemon.LatticeModule
using TensorOperations
using TensorMatrices_lemon

dimM = 5
trotterparameter = 0.01
trotteriteration = 50

isinginfo = SpinInfo("quantum_ising_2",1.0e-13,0.0)
trotterinfo = TrotterInfo(trotterparameter, trotteriteration)
ising = QuantumIsingModel(isinginfo, trotterinfo)
lattice = buildLattice("quantum_2d_square", ising)
simulator = buildSimulator(lattice,dimM)

fieldrange = linspace(1.0,3.5,20)
simulatorQuantum(fieldrange, simulator)
#tensorT = getTensorT(simulator)[1]
#tenmatMMd = getTenmatMMd(tensorT)
#Ul, lambdaVector, Uld = svd(tenmatMMd.matrix)
#trunUl = truncMatrixU(Ul,dimM)
#tensorU = matU2tenU(trunUl)


text:
module FreqTables
    using CategoricalArrays
    using Tables
    using NamedArrays
    using Missings

    include("freqtable.jl")

    export freqtable, proptable, prop, Name
end # module


text:
# ---
# title: JuliaRL\_SAC\_Pendulum
# cover: assets/JuliaRL_SAC_Pendulum.png
# description: SAC applied to Pendulum
# date: 2021-05-22
# author: "[Roman Bange](https://github.com/rbange)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses
using IntervalSets

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:SAC},
    ::Val{:Pendulum},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    rng = StableRNG(seed)
    inner_env = PendulumEnv(T = Float32, rng = rng)
    action_dims = inner_env.n_actions
    A = action_space(inner_env)
    low = A.left
    high = A.right
    ns = length(state(inner_env))
    na = 1

    env = ActionTransformedEnv(
        inner_env;
        action_mapping = x -> low + (x[1] + 1) * 0.5 * (high - low),
    )
    init = glorot_uniform(rng)

    create_policy_net() = NeuralNetworkApproximator(
        model = GaussianNetwork(
            pre = Chain(
                Dense(ns, 30, relu, init = init),
                Dense(30, 30, relu, init = init),
            ),
            μ = Chain(Dense(30, na, init = init)),
            logσ = Chain(Dense(30, na, x -> clamp.(x, typeof(x)(-10), typeof(x)(2)), init = init)),
        ),
        optimizer = ADAM(0.003),
    )

    create_q_net() = NeuralNetworkApproximator(
        model = Chain(
            Dense(ns + na, 30, relu; init = init),
            Dense(30, 30, relu; init = init),
            Dense(30, 1; init = init),
        ),
        optimizer = ADAM(0.003),
    )

    agent = Agent(
        policy = SACPolicy(
            policy = create_policy_net(),
            qnetwork1 = create_q_net(),
            qnetwork2 = create_q_net(),
            target_qnetwork1 = create_q_net(),
            target_qnetwork2 = create_q_net(),
            γ = 0.99f0,
            τ = 0.005f0,
            α = 0.2f0,
            batch_size = 64,
            start_steps = 1000,
            start_policy = RandomPolicy(Space([-1.0..1.0 for _ in 1:na]); rng = rng),
            update_after = 1000,
            update_freq = 1,
            automatic_entropy_tuning = true,
            lr_alpha = 0.003f0,
            action_dims = action_dims,
            rng = rng,
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = 10000,
            state = Vector{Float32} => (ns,),
            action = Vector{Float32} => (na,),
        ),
    )

    stop_condition = StopAfterStep(10_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "# Play Pendulum with SAC")
end

#+ tangle=false
using Plots
pyplot() #hide
ex = E`JuliaRL_SAC_Pendulum`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_SAC_Pendulum.png") #hide

# ![](assets/JuliaRL_SAC_Pendulum.png)


text:
using HorizonSideRobots
function main(r::Robot)
    ugol(Nord,Ost)
end

function ugol(west1,west2)

    while (isborder(r,west1)==false) || (isborder(r,west2)==false)
        while isborder(r,west1)==false
            moves!(r,west1)
        end
        while isborder(r,west2)==false
            moves!(r,west2)
        end
        #move!(r,west2)
        #move!(r,west1)
        #move!(r,inverse(west1))
    end
end

function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

r=Robot("5.sit", animate=true)

max_stars_repo_name: Savko-fokin/mirea-progs

max_issues_repo_name: Savko-fokin/mirea-progs

max_forks_repo_name: Savko-fokin/mirea-progs


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
struct F{p, T<:Integer} <: Integer
    a::T
    F{p, T}(a::Integer) where {p, T<:Integer} = new{p, T}(mod(T(a), p))
end
F{p}(a::Integer) where p = F{p, typeof(a)}(mod(a, p))
F{p, S}(x::F{p, T}) where {S<:Integer, p, T<:Integer} = F{p, S}(x.a)

Base.promote_rule(::Type{F{p, T}}, ::Type{S}) where {p, T<:Integer, S<:Integer} =
    F{p, promote_type(T, S)}

Base.zero(::Type{F{p, T}}) where {p, T<:Integer} = F{p}(mod(zero(T), p))
Base.one(::Type{F{p, T}}) where {p, T<:Integer} = F{p}(mod(one(T), p))
for op in (:-, :+)
    @eval Base.$op(x::F{p}) where p = F{p}(mod($op(x.a), p))
end
for op in (:-, :+, :*)
    @eval Base.$op(x::F{p}, y::F{p}) where p = F{p}(mod($op(x.a, y.a), p))
end
Base.inv(x::F{p}) where p = F{p}(invmod(x.a, p))
Base.:/(x::F{p}, y::F{p}) where p = x * inv(y)
Base.:\(x::F{p}, y::F{p}) where p = inv(x) * y
Base.:(==)(x::F{p}, y::F{p}) where p = x.a == y.a
Base.:<(x::F{p}, y::F{p}) where p = x.a < y.a

Base.show(io::IO, x::F{p}) where p = print(io, "F", p, '(', x.a, ')')

F7 = F{7, BigInt}
x, y = F7(10), F7(-2)
@show(x, y, zero(x), one(x), +x, -x, x + y, x - y, x * y, x / y, x \ y, x^3, x^-5, x == F7(3), x == 3)
println()
A = F7[1 2; 3 4]
@show(A, 4A, A/4)
println()
using LinearAlgebra
L, U = lu(A, NoPivot())
@show(L, U, L * U, det(A), inv(A), inv(A) * A, A * inv(A));

# %%
"""See https://docs.julialang.org/en/v1/manual/interfaces/"""
Base.iterate(Fp::Type{F{p, T}}) where {p, T<:Integer} = (zero(Fp), zero(T))
function Base.iterate(Fp::Type{F{p, T}}, state) where {p, T<:Integer}
    nextstate = state + 1
    nextstate < p ? (Fp(nextstate), nextstate) : nothing
end
Base.IteratorSize(Fp::Type{F{p, T}}) where {p, T<:Integer} = Base.HasLength()
Base.length(Fp::Type{F{p, T}}) where {p, T<:Integer} = p
Base.eltype(Fp::Type{F{p, T}}) where {p, T<:Integer} = Fp

squares(Fp) = Fp[x^2 for x in Fp]
squareroots(k, Fp) = Fp[x for x in Fp if x^2 == k]
@show(collect(F7), squares(F7), squareroots.(0:6, Ref(F7)));

# %%
using GaloisFields, LinearAlgebra
GF7 = @GaloisField 7
B = GF7[1 2; 3 4]
det(B)

# %%
inv(B)

# %%
lu(B)

# %%
using AbstractAlgebra

@show GF7 = GF(7)
x, y = GF7(10), GF7(-2)
@show(x, y, zero(x), one(x), -x, x + y, x - y, x * y, x^3, x^-5, x == GF7(3), x == 3)
println()

@show C = GF7[1 2; 3 4]
@show P, x = PolynomialRing(GF7, "x")
@show(det(C), inv(C), lu(C), charpoly(P, C))
println()

squares(Fp) = [x^2 for x in Fp]
squareroots(k, Fp) = [x for x in Fp if x^2 == k]
@show(collect(GF7), squares(GF7), squareroots.(0:6, Ref(GF7)));

# %%

max_stars_repo_path: 0010/minimal def of a prime field.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0010/minimal def of a prime field.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0010/minimal def of a prime field.jl

max_forks_repo_name: genkuroki/public


text:
const HEADERS_JSON = [
    "Access-Control-Allow-Origin" => "*",
    "Access-Control-Allow-Headers" => "*",
    "Access-Control-Allow-Methods" => "*",
    "Content-Type" => "application/json"
]

const HEADERS_TEXT = [
    "Access-Control-Allow-Origin" => "*",
    "Access-Control-Allow-Headers" => "*",
    "Access-Control-Allow-Methods" => "*",
    "Content-Type" => "text/html"
]

function printroutes(app::App, ip::IPAddr, port::Int)
    if ip == Sockets.localhost
        println("""\n
          App running at:
          http://localhost:$port/

          Method           Route
          ==============================================
        """)
    else
        println("""\n
          App running at:
          http://$ip:$port/

          Method           Route
          ==============================================
        """)
    end

    for (method, routes) in app.routelist
        for route in routes
            println("  $(rpad(method, 17))", route)
        end
    end
end

function _body(::JSONParser, req::HTTP.Request)
    body = IOBuffer(HTTP.payload(req))
    if eof(body)
        return JSON3.Object()
    else
        return JSON3.read(body)
    end
end

function _body(::TextParser, req::HTTP.Request)
    body = HTTP.payload(req)
    return String(body)
end

_response(res_body) = HTTP.Response(200, HEADERS_TEXT, body = string(res_body))

_response(res_body::AbstractString) = HTTP.Response(200, HEADERS_TEXT, body = res_body)

_response(res_body::Union{Dict, NamedTuple}) = HTTP.Response(200, HEADERS_JSON, body = JSON3.write(res_body))

serve(app::App) = serve(app, Sockets.localhost, 8081)

serve(app::App, port::Int) = serve(app, Sockets.localhost, port)

serve(app::App, ip::IPAddr) = serve(app, ip, 8081)

serve(app::App, ip::AbstractString) = serve(app, parse(IPAddr, ip), 8081)

serve(app::App, ip::AbstractString, port::Int) = serve(app, parse(IPAddr, ip), port)

function serve(app::App, ip::IPAddr, port::Int)
    printroutes(app, ip, port)

    HTTP.serve(ip, port) do req::HTTP.Request
        req_body = _body(app.reqparser, req)
        res_body = HTTP.handle(app.router, req, req_body)
        return _response(res_body)
    end
end


text:
using JuMP: @variable, Model

@testset "JuMP variables" begin
    @polyvar x y z
    m = Model()
    @variable(m, q[1:2])

    p = @inferred q[2] * y + q[1] * x
    @test @inferred(q' * [x, y]) == p
    @test @inferred([x, y]' * q) == p

    @test @inferred((x + y) * p) isa Polynomial{JuMP.AffExpr}

    m = Model()
    @variable(m, Q[1:2, 1:2])
    @test @inferred([x, y]' * Q * [x, y]) == Q[1, 1] * x^2 + Q[2, 1] * x * y + Q[1, 2] * x * y + Q[2, 2] * y^2
end


text:
module PolyesterForwardDiff

using Polyester
import ForwardDiff

const DiffResult = ForwardDiff.DiffResults.DiffResult

function cld_fast(a::A,b::B) where {A,B}
    T = promote_type(A,B)
    cld_fast(a%T,b%T)
end
function cld_fast(n::T, d::T) where {T}
    x = Base.udiv_int(n, d)
    x += n != d*x
end

store_val!(r::Base.RefValue{T}, x::T) where {T} = (r[] = x)
store_val!(r::Ptr{T}, x::T) where {T} = Base.unsafe_store!(r, x)

function evaluate_chunks!(f::F, (r,Δx,x), start, stop, ::ForwardDiff.Chunk{C}) where {F,C}
    cfg = ForwardDiff.GradientConfig(f, x, ForwardDiff.Chunk{C}(), nothing)
    N = length(x)
    last_stop = cld_fast(N, C)
    is_last = last_stop == stop
    stop -= is_last

    xdual = cfg.duals
    seeds = cfg.seeds
    ForwardDiff.seed!(xdual, x)
    for c ∈ start:stop
        i = (c-1) * C + 1
        ForwardDiff.seed!(xdual, x, i, seeds)
        ydual = f(xdual)
        ForwardDiff.extract_gradient_chunk!(Nothing, Δx, ydual, i, C)
        ForwardDiff.seed!(xdual, x, i)
    end
    if is_last
        lastchunksize = C + N - last_stop*C
        lastchunkindex = N - lastchunksize + 1
        ForwardDiff.seed!(xdual, x, lastchunkindex, seeds, lastchunksize)
        _ydual = f(xdual)
        ForwardDiff.extract_gradient_chunk!(Nothing, Δx, _ydual, lastchunkindex, lastchunksize)
        store_val!(r, ForwardDiff.value(_ydual))
    end
end

function threaded_gradient!(f::F, Δx::AbstractVector, x::AbstractVector, ::ForwardDiff.Chunk{C}) where {F,C}
    N = length(x)
    d = cld_fast(N, C)
    r = Ref{eltype(Δx)}()
    batch((d,min(d,num_threads())), r, Δx, x) do rΔxx,start,stop
        evaluate_chunks!(f, rΔxx, start, stop, ForwardDiff.Chunk{C}())
    end
    r[]
end

#### in-place jac, out-of-place f ####

function evaluate_jacobian_chunks!(f::F, (Δx,x), start, stop, ::ForwardDiff.Chunk{C}) where {F,C}
    cfg = ForwardDiff.JacobianConfig(f, x, ForwardDiff.Chunk{C}(), nothing)

    # figure out loop bounds
    N = length(x)
    last_stop = cld_fast(N, C)
    is_last = last_stop == stop
    stop -= is_last

    # seed work arrays
    xdual = cfg.duals
    ForwardDiff.seed!(xdual, x)
    seeds = cfg.seeds

    # handle intermediate chunks
    for c ∈ start:stop
        # compute xdual
        i = (c-1) * C + 1
        ForwardDiff.seed!(xdual, x, i, seeds)

        # compute ydual
        ydual = f(xdual)

        # extract part of the Jacobian
        Δx_reshaped = ForwardDiff.reshape_jacobian(Δx, ydual, xdual)
        ForwardDiff.extract_jacobian_chunk!(Nothing, Δx_reshaped, ydual, i, C)
        ForwardDiff.seed!(xdual, x, i)
    end

    # handle the last chunk
    if is_last
        lastchunksize = C + N - last_stop*C
        lastchunkindex = N - lastchunksize + 1

        # compute xdual
        ForwardDiff.seed!(xdual, x, lastchunkindex, seeds, lastchunksize)

        # compute ydual
        _ydual = f(xdual)

        # extract part of the Jacobian
        _Δx_reshaped = ForwardDiff.reshape_jacobian(Δx, _ydual, xdual)
        ForwardDiff.extract_jacobian_chunk!(Nothing, _Δx_reshaped, _ydual, lastchunkindex, lastchunksize)
    end
end

function threaded_jacobian!(f::F, Δx::AbstractArray, x::AbstractArray, ::ForwardDiff.Chunk{C}) where {F,C}
    N = length(x)
    d = cld_fast(N, C)
    batch((d,min(d,num_threads())), Δx, x) do Δxx,start,stop
        evaluate_jacobian_chunks!(f, Δxx, start, stop, ForwardDiff.Chunk{C}())
    end
    return Δx
end

# # #### in-place jac, in-place f ####

function evaluate_f_and_jacobian_chunks!(f!::F, (y,Δx,x), start, stop, ::ForwardDiff.Chunk{C}) where {F,C}
    cfg = ForwardDiff.JacobianConfig(f!, y, x, ForwardDiff.Chunk{C}(), nothing)

    # figure out loop bounds
    N = length(x)
    last_stop = cld_fast(N, C)
    is_last = last_stop == stop
    stop -= is_last

    # seed work arrays
    ydual, xdual = cfg.duals
    ForwardDiff.seed!(xdual, x)
    seeds = cfg.seeds
    Δx_reshaped = ForwardDiff.reshape_jacobian(Δx, ydual, xdual)

    # handle intermediate chunks
    for c ∈ start:stop
        # compute xdual
        i = (c-1) * C + 1
        ForwardDiff.seed!(xdual, x, i, seeds)

        # compute ydual
        f!(ForwardDiff.seed!(ydual, y), xdual)

        # extract part of the Jacobian
        ForwardDiff.extract_jacobian_chunk!(Nothing, Δx_reshaped, ydual, i, C)
        ForwardDiff.seed!(xdual, x, i)
    end

    # handle the last chunk
    if is_last
        lastchunksize = C + N - last_stop*C
        lastchunkindex = N - lastchunksize + 1

        # compute xdual
        ForwardDiff.seed!(xdual, x, lastchunkindex, seeds, lastchunksize)

        # compute ydual
        f!(ForwardDiff.seed!(ydual, y), xdual)

        # extract part of the Jacobian
        ForwardDiff.extract_jacobian_chunk!(Nothing, Δx_reshaped, ydual, lastchunkindex, lastchunksize)
        map!(ForwardDiff.value, y, ydual)
    end
end

function threaded_jacobian!(f!::F, y::AbstractArray, Δx::AbstractArray, x::AbstractArray, ::ForwardDiff.Chunk{C}) where {F,C}
    N = length(x)
    d = cld_fast(N, C)
    batch((d,min(d,num_threads())), y, Δx, x) do yΔxx,start,stop
        evaluate_f_and_jacobian_chunks!(f!, yΔxx, start, stop, ForwardDiff.Chunk{C}())
    end
    Δx
end

end

max_stars_repo_name: JuliaDiff/PolyesterForwardDiff

max_issues_repo_name: ChrisRackauckas/PolyesterForwardDiff

max_forks_repo_name: JuliaDiff/PolyesterForwardDiff


text:
function bondtensor(::Type{TT}, J) where TT
    TT.([J -J; -J J])
end

function vertextensor(::Type{TT}, h) where TT
    TT.([h, -h])
end

function copytensor(::Type{T}) where T
    Diagonal([one(T), zero(T), zero(T), one(T)])
end

function resettensor(::Type{T}) where T
    [one(T) one(T); zero(T) zero(T)]
end

"""
    hypercubicI([T], ndim::Int, D::Int)

Get a `ndim`-dimensional identity hypercubic, with bound dimension `D`.
"""
function hypercubicI(::Type{T}, ndim::Int, D::Int) where T
    res = zeros(T, fill(D, ndim)...)
    @inbounds for i=1:D
        res[fill(i,ndim)...] = one(T)
    end
    return res
end

hypercubicI(ndim::Int, D::Int) = hypercubicI(Float64, ndim, D)

Gh(vertex_tensor::Vector{T}) where T = tropicalblock(Diagonal(vertex_tensor))
Gvb(bond_tensor::Matrix{T}) where T = tropicalblock(Diagonal([bond_tensor...]))
Ghb(bond_tensor::Matrix{T}) where T = tropicalblock(transpose(bond_tensor))
function G16(::Type{TT}, Js) where TT<:TropicalTypes
    xs = map(x->_bondtensor(TT, x), Js)
    mat = zeros(TT, 2, 2, 2, 2, 2, 2, 2, 2)
    for a=1:2, α=1:2, β=1:2, b=1:2, γ=1:2, c=1:2, d=1:2, δ=1:2
        mat[a,b,c,d,α,β,γ,δ] += xs[1][a,α] * xs[2][a,β] * xs[3][a,γ] * xs[4][a,δ] *
        xs[5][b,α] * xs[6][b,β] * xs[7][b,γ] * xs[8][b,δ] * xs[9][c,α] * xs[10][c,β] *
        xs[11][c,γ] * xs[12][c,δ] * xs[13][d,α] * xs[14][d,β] * xs[15][d,γ] * xs[16][d,δ]
    end
    tropicalblock(reshape(mat, 16, 16))
end

function _bondtensor(::Type{TT}, J) where TT
    TT.([J -J; -J J])
end

"""
    Gcp(T)

copy state of qubit 2 -> 1.
"""
function Gcp(::Type{TT}) where TT
    tropicalblock(copytensor(TT))
end

function Greset(::Type{TT}) where TT
    tropicalblock([one(TT) one(TT); zero(TT) zero(TT)])
end

function Gcut(::Type{TT}) where TT
    tropicalblock([one(TT) one(TT); one(TT) one(TT)])
end


text:
"""
$(TYPEDEF)
"""
abstract type AbstractElementGeometry <: AbstractExtendableGridApexType end

"""
$(TYPEDSIGNATURES)

List supported element geometries.
"""
elementgeometries()=AbstractTrees.print_tree(AbstractElementGeometry,5,indicate_truncation=false)

"""
$(TYPEDEF)
"""
abstract type AbstractElementGeometry0D <: AbstractElementGeometry end
"""
$(TYPEDEF)
"""
abstract type AbstractElementGeometry1D <: AbstractElementGeometry end
"""
$(TYPEDEF)
"""
abstract type AbstractElementGeometry2D <: AbstractElementGeometry end
"""
$(TYPEDEF)
"""
abstract type AbstractElementGeometry3D <: AbstractElementGeometry end
"""
$(TYPEDEF)
"""
abstract type AbstractElementGeometry4D <: AbstractElementGeometry end

"""
$(TYPEDEF)
"""
abstract type Vertex0D <: AbstractElementGeometry0D end

"""
$(TYPEDEF)
"""
abstract type Edge1D <: AbstractElementGeometry1D end

"""
$(TYPEDEF)
"""
abstract type Polygon2D <: AbstractElementGeometry2D end
"""
$(TYPEDEF)
"""
abstract type Triangle2D <: Polygon2D end
"""
$(TYPEDEF)
"""
abstract type Quadrilateral2D <: Polygon2D end
"""
$(TYPEDEF)
"""
abstract type Pentagon2D <: Polygon2D end
"""
$(TYPEDEF)
"""
abstract type Hexagon2D <: Polygon2D end
"""
$(TYPEDEF)
"""
abstract type Parallelogram2D <: Quadrilateral2D end

"""
$(TYPEDEF)
"""
abstract type Circle2D <: AbstractElementGeometry2D end

"""
$(TYPEDEF)
"""
abstract type Polyhedron3D <: AbstractElementGeometry3D end
"""
$(TYPEDEF)
"""
abstract type Tetrahedron3D <: Polyhedron3D end
"""
$(TYPEDEF)
"""
abstract type Hexahedron3D <: Polyhedron3D end
"""
$(TYPEDEF)
"""
abstract type Parallelepiped3D <: Hexahedron3D end
"""
$(TYPEDEF)
"""
abstract type Prism3D <: Polyhedron3D end
"""
$(TYPEDEF)
"""
abstract type TrianglePrism3D <: Prism3D end
"""
$(TYPEDEF)
"""
abstract type Sphere3D <: AbstractElementGeometry3D end

"""
$(TYPEDEF)
"""
abstract type Polychoron4D <: AbstractElementGeometry4D end
"""
$(TYPEDEF)
"""
abstract type HyperCube4D <: AbstractElementGeometry4D end

"""
$(TYPEDSIGNATURES)
"""
dim_element(::Type{<:AbstractElementGeometry0D})=0
"""
$(TYPEDSIGNATURES)
"""
dim_element(::Type{<:AbstractElementGeometry1D})=1
"""
$(TYPEDSIGNATURES)
"""
dim_element(::Type{<:AbstractElementGeometry2D})=2
"""
$(TYPEDSIGNATURES)
"""
dim_element(::Type{<:AbstractElementGeometry3D})=3
"""
$(TYPEDSIGNATURES)
"""
dim_element(::Type{<:AbstractElementGeometry4D})=4


text:
type TDPrimes{T<:Integer}
    plim::T
end

function Base.start{T<:Integer}(pl::TDPrimes{T})
    2ones(T, 1)
end

function Base.done{T<:Integer}(pl::TDPrimes{T}, p::Array{T,1})
    p[end] > pl.plim
end

function Base.next{T<:Integer}(pl::TDPrimes{T}, p::Array{T,1})
    pr = p[end]
    for i in (pr+1):(pl.plim)
        ispr = true
        for j in p
            if i%j == 0
                ispr = false
                break
            end
        end
        if ispr
            push!(p, i)
            return (pr, p)
        end
    end
    push!(p, typemax(T))
    return (pr, p)
end

n = 100
print("The primes <= ", n, " are:\n    ")

for i in TDPrimes(n)
    print(i, " ")
end
println()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module LogDensities

using SparseQuadratureGrids, StaticArrays, ConstrainedParameters, DiffBase, ForwardDiff, Optim

import  Base: show, getindex, setindex!, size, IndexStyle, convert, length
import  ConstrainedParameters: type_length, param_type_length, construct, log_jacobian, update!
import  Optim:  Options,
                initial_state,
                update_state!,
                NewtonTrustRegionState,
                solve_tr_subproblem!,
                assess_convergence,
                value,
                gradient
#import  DiffBase: MutableDiffResult

export  Data,
        parameters,
        CovarianceMatrix,
        PositiveVector,
        ProbabilityVector,
        RealVector,
        Model,
        ModelParam,
        construct,
        log_density,
        log_jacobian!,
        quad_form,
        inv_det,
        inv_root_det,
        root_det,
        log_root_det,
        trace_inverse,
        lpdf_InverseWishart,
        lpdf_normal,
        logit,
        logistic,
        update!,
        type_length,
        ModelRank,
        StaticRank,
        DynamicRank,
        Full,
        FixedRank,
        LDR,
        Dynamic,
        MarginalBuffer,
        optimize!,
        ModelDiff

include("model_utilities.jl")
include("parameter_methods.jl")
include("model_derivatives.jl")
include("model.jl")

end # module


text:
struct MultiDiscreteNonParametricSampler{T <: Real, S <: AbstractVector{T}, A <: AliasTable} <: Sampleable{Multivariate,Discrete}
    support::Vector{S}
    probabilities::S
    aliastable::A

    function MultiDiscreteNonParametricSampler{T,S}(support::Vector{S}, probs::AbstractVector{<:Real}) where {T <: Real, S<:AbstractVector{T}}
        isempty(support) && "Empty support given."
        aliastable = AliasTable(probs)
        new{T,S,typeof(aliastable)}(support, convert(S, probs), aliastable)
    end
end

function MultiDiscreteNonParametricSampler(support::Vector{S}, probs::AbstractVector{<:Real}) where {T <: Real, S <: AbstractVector{T}}
    return MultiDiscreteNonParametricSampler{T,S}(support, probs)
end

function Base.length(s::MultiDiscreteNonParametricSampler)
    return length(s.support[1])
end
Base.eltype(::Type{<:MultiDiscreteNonParametricSampler}) = Float64

function Distributions._rand!(rng::AbstractRNG, s::MultiDiscreteNonParametricSampler, x::AbstractVector)
    @inbounds x .= s.support[rand(rng, s.aliastable)]
end
"""
    SMPSSampler

Sampler object for SMPS scenarios. Obtained by reading from a model defined in SMPS format.

See also: [`SMPSScenario`](@ref)
"""
struct SMPSSampler{T <: AbstractFloat, M <: AbstractMatrix} <: AbstractSampler{SMPSScenario{T,M}}
    template::LPData{T,M}
    technology::UnitRange{Int}
    recourse::UnitRange{Int}
    random_variables::Dict{RowCol, Sampleable}
    random_vectors::Dict{Vector{RowCol}, Sampleable}
    inclusions::Dict{RowCol, Symbol}
end

function finite_support(sampler::SMPSSampler{T}) where T <: AbstractFloat
    return all(ran_var -> ran_var isa DiscreteNonParametric, values(sampler.random_variables)) &&
        all(ran_vec -> ran_vec isa MultiDiscreteNonParametricSampler, values(sampler.random_vectors))
end

function support_size(sampler)
    ran_var_sizes = map(values(sampler.random_variables)) do ran_var
        return Float64(length(ran_var.support))
    end
    ran_vec_sizes = map(values(sampler.random_vectors)) do ran_vec
        return Float64(length(ran_vec))
    end
    return reduce(*, vcat(ran_var_sizes, ran_vec_sizes))
end

function SMPSSampler(sto::RawStoch{N}, stage::SMPSStage{T}) where {N, T <: AbstractFloat}
    2 <= stage.id <= N + 1 || error("$(stage.id) not in range 2 to $(N + 1).")
    random_variables = Dict{RowCol, Sampleable}()
    random_vectors = Dict{Vector{RowCol}, Sampleable}()
    inclusions = Dict{RowCol, Symbol}()
    for ran_var in sto.random_variables[stage.id - 1]
        inclusions[ran_var.rowcol] = ran_var.inclusion
        if ran_var isa IndepDiscrete
            random_variables[ran_var.rowcol] =
                DiscreteNonParametric(ran_var.support, ran_var.probabilities)
        elseif ran_var isa IndepDistribution
            if ran_var.distribution == UNIFORM
                random_variables[ran_var.rowcol] =
                    Uniform(first(ran_var.parameters), second(ran_var.parameters))
            elseif ran_var.distribution == NORMAL
                random_variables[ran_var.rowcol] =
                    Normal(first(ran_var.parameters), second(ran_var.parameters))
            elseif ran_var.distribution == GAMMA
                random_variables[ran_var.rowcol] =
                    Gamma(first(ran_var.parameters), second(ran_var.parameters))
            elseif ran_var.distribution == BETA
                random_variables[ran_var.rowcol] =
                    Beta(first(ran_var.parameters), second(ran_var.parameters))
            elseif ran_var.distribution == LOGNORM
                random_variables[ran_var.rowcol] =
                    LogNormal(first(ran_var.parameters), second(ran_var.parameters))
            end
        end
    end
    for ran_vec in sto.random_vectors[stage.id - 1]
        if ran_vec isa BlockDiscrete
            isempty(ran_vec.support[1]) && error("Block $block has empty support.")
            rowcols = Vector{RowCol}()
            support = Vector{Vector{T}}()
            push!(support, Vector{T}())
            for (rowcol, val) in ran_vec.support[1]
                push!(rowcols, rowcol)
                push!(support[end], val)
                inclusions[rowcol] = ran_vec.inclusion
            end
            for remaining in ran_vec.support[2:end]
                push!(support, Vector{T}())
                for (rowcol, val) in ran_vec.support[1]
                    if haskey(remaining, rowcol)
                        push!(support[end], remaining[rowcol])
                    else
                        push!(support[end], val)
                    end
                end
            end
            random_vectors[rowcols] = MultiDiscreteNonParametricSampler(support, ran_vec.probabilities)
        end
    end
    return SMPSSampler(stage.uncertain,
                       stage.technology,
                       stage.recourse,
                       random_variables,
                       random_vectors,
                       inclusions)
end

function (sampler::SMPSSampler{T})() where T <: AbstractFloat
    # Collect samples
    samples = Dict{RowCol,T}()
    for (rowcol, ran_var) in sampler.random_variables
        samples[rowcol] = rand(ran_var)
    end
    for (rowcols, ran_vec) in sampler.random_vectors
        block_sample = rand(ran_vec)
        for (idx, rowcol) in enumerate(rowcols)
            samples[rowcol] = block_sample[idx]
        end
    end
    return create_scenario(sampler, samples)
end

function full_support(sampler::SMPSSampler{T}) where T <: AbstractFloat
    # Sanity check
    finite_support(sampler) || error("Sampler does not have finite support, cannot return full support.")
    # Warn if support is large
    nscenarios = support_size(sampler)
    if nscenarios > 1e5
        @warn "Generating full support of $nscenarios scenarios."
    end
    # Collect all realizations
    all_realizations = Vector{Vector{Tuple{T,Dict{RowCol,T}}}}()
    # Collect the full support of every random variable
    for (rowcol, ran_var) in sampler.random_variables
        push!(all_realizations, Vector{Tuple{T,Dict{RowCol,T}}}())
        for (π,x) in zip(ran_var.p, ran_var.support)
            push!(all_realizations[end], (π, Dict(rowcol => x)))
        end
    end
    # Collect the full support of every random vector
    for (rowcols, ran_vec) in sampler.random_vectors
        push!(all_realizations, Vector{Tuple{T,Dict{RowCol,T}}}())
        for (π,vec) in zip(ran_vec.probabilities, ran_vec.support)
            d = Dict{RowCol,T}()
            for (rowcol,x) in zip(rowcols,vec)
                d[rowcol] = x
            end
            push!(all_realizations[end], (π, d))
        end
    end
    # Generate scenario for every realization combination
    return mapreduce(vcat, Base.Iterators.product(all_realizations...)) do realization
        π = 1.0
        samples = Dict{RowCol,T}()
        for (ρ, outcome) in realization
            samples = merge(samples, outcome)
            π *= ρ
        end
        return create_scenario(sampler, samples; π = π)
    end
end

function create_scenario(sampler::SMPSSampler{T}, samples::Dict{RowCol,T}; π::AbstractFloat = 1.0) where T <: AbstractFloat
    # Prepare scenario data
    Δq  = copy(sampler.template.c₁)
    A   = copy(sampler.template.A)
    Δd₁ = copy(sampler.template.d₁)
    ΔC  = copy(sampler.template.C)
    Δd₂ = copy(sampler.template.d₂)
    Δh  = copy(sampler.template.b)
    # Fill scenario data
    for (rowcol, ξ) in samples
        (row, col) = rowcol
        (i,j,type) = sampler.template.indexmap[rowcol]
        if type == OBJ
            if sampler.inclusions[rowcol] == MULTIPLY
                Δq[j] += abs(Δq[j]) * ξ
            else
                Δq[j] += ξ
            end
        elseif type == EQ
            if col == RHS
                if sampler.inclusions[rowcol] == MULTIPLY
                    Δh[i] += abs(Δh[i]) * ξ
                else
                    Δh[i] += ξ
                end
            else
                if sampler.inclusions[rowcol] == MULTIPLY
                    A[i,j] += abs(A[i,j]) * ξ
                else
                    A[i,j] += ξ
                end
            end
        elseif type == LEQ
            if col == RHS
                if sampler.inclusions[rowcol] == MULTIPLY
                    Δd₂[i] += abs(Δd₂[i]) * ξ
                else
                    Δd₂[i] += ξ
                end
            else
                if sampler.inclusions[rowcol] == MULTIPLY
                    ΔC[i,j] += abs(ΔC[i,j]) * ξ
                else
                    ΔC[i,j] += ξ
                end
            end
        elseif type == GEQ
            if col == RHS
                if sampler.inclusions[rowcol] == MULTIPLY
                    Δd₁[i] += abs(d₁[i]) * ξ
                else
                    Δd₁[i] += ξ
                end
            else
                if sampler.inclusions[rowcol] == MULTIPLY
                    ΔC[i,j] += abs(ΔC[i,j]) * ξ
                else
                    ΔC[i,j] += ξ
                end
            end
        elseif type == RANGE
            if sampler.inclusions[rowcol] == MULTIPLY
                Δd₁[i] += abs(d₁) * ξ
                Δd₂[i] += abd(d₂) * ξ
            else
                Δd₁[i] += ξ
                Δd₂[i] += ξ
            end
        end
    end
    ΔT = A[:,sampler.technology]
    ΔW = A[:,sampler.recourse]
    return SMPSScenario(Probability(π), Δq, ΔT, ΔW, Δh, ΔC, Δd₁, Δd₂)
end


text:
macro isbn_str(s)
    ISBN(s)
end

struct ISBN <: AbstractString
    s::AbstractString

    ISBN(s) = verify(s) ? new(replace(s, "-", "")) : throw(ArgumentError("invalid ISBN string: $s"))
end

string(s::ISBN) = s.s

function verify(s::AbstractString)
    s = replace(s, "-", "")
    chars = split(s, "")

    length(chars) == 10 || return false

    if last(chars) == "X"
        chars[end] = "10"
    end

    all(c -> all(isdigit, c), chars) || return false

    sum(parse(Int, c) * i for (c, i) in zip(chars, 10:-1:1)) % 11 == 0
end

Base.isvalid(::Type{ISBN}, s::AbstractString) = verify(s)


text:
using SHA

# Script to download the leap seconds list from NIST, and format it into julia
# code for use in Chrono.jl.

#-------------------------------------------------------------------------------
# Download leap-seconds.list
# Official URL (?)
url = "ftp://time.nist.gov/pub/leap-seconds.list"
# IETF has a mirror if the NIST site is unresponsive, which it commonly is.
# The mirror is not as up to date though.
# url = "https://www.ietf.org/timezones/data/leap-seconds.list"
local_name = "leap-seconds.list"
if length(ARGS) < 1 || ARGS[1] != "false"
    download(url, local_name)
end
lines = readlines(local_name)

#-------------------------------------------------------------------------------
# Parse file.
rawdata = [] # for hash computation
last_update_ntp = 0
expiry_ntp = 0
leap_epochs_ntp = Int[]
leap_secs = Int[]
expected_sha = ""
for line in lines
    # Mostly, lines starting with # are comments.  However, due to
    # terribleness, lines starting with "#@", "#$" or "#h" also contain valid
    # data which needs to be parsed for hash verification.
    if startswith(line, "#\$")
        datastr = match(r"^#\$\s*(\d+)\s*$", line)[1]
        last_update_ntp = parse(Int64, datastr)
        push!(rawdata, datastr)
    elseif startswith(line, "#@")
        datastr = match(r"^#@\s*(\d+)\s*$", line)[1]
        expiry_ntp = parse(Int64, datastr)
        push!(rawdata, datastr)
    elseif startswith(line, "#h")
        # The specified hash method for leap-seconds.list is super broken,
        # since it removes semantically meaningful whitespace.  At least we can
        # guard against transmission errors...
        #
        # There's also a bug in writing out the expected hash!  Work around it
        # by parsing and reformatting, ugh.
        pieces = map(x->@sprintf("%08x", parse(Int,x,16)), split(line[3:end]))
        expected_sha = join(pieces)
    elseif !startswith(line, "#")
        m = match(r"(\d+)\s*(\d+)\s*#.*", line)
        push!(leap_epochs_ntp, parse(Int64, m[1]))
        push!(leap_secs, parse(Int, m[2]))
        push!(rawdata, m[1])
        push!(rawdata, m[2])
    end
end

actual_sha = bytes2hex(SHA.sha1(join(rawdata)))
@show expected_sha actual_sha
@assert expected_sha == actual_sha
@assert all(diff(leap_secs) .== 1)

#-------------------------------------------------------------------------------
# Generate julia code
function format_ntp_secs_as_date(ntp_epoch_secs)
    d = Dates.Second(ntp_epoch_secs) + Dates.DateTime(1900)
    @assert Dates.hour(d) == 0 && Dates.minute(d) == 0 &&
            Dates.second(d) == 0 && Dates.millisecond(d) == 0
    "Date($(Dates.year(d)), $(Dates.month(d)), $(Dates.day(d)))"
end

genname = "leap_seconds.jl"
open(genname, "w") do file
    write(file,
        """
        # Automatically generated by update_leap_seconds.jl - do not edit!

        const leap_index_offset = $(leap_secs[1]-1) # For use with searchsortedlast()

        const utc_offsets = [
            typemin(Date),  # For use with searchsortedlast()
            $(join(map(format_ntp_secs_as_date, leap_epochs_ntp[2:end]), ",\n    "))
        ]

        const leap_secs_expiry_date = $(format_ntp_secs_as_date(expiry_ntp))
        """
    )
end

print(readstring(genname))