text:
xlocations(ex) = Expr(:call, :($YaoLocations.Locations), ex)
xctrl_locations(ex) = Expr(:call, :($YaoLocations.CtrlLocations), ex)

"""
    @gate <locs> => <gate>

Syntax sugar for `apply(gate, locs)`, must be used inside `@device`.
See also [`@device`](@ref).

!!! tips
    You don't have to write `@gate` in most cases inside `@device`.
    But in case there is ambuigity, you can annotate the expression
    with `@gate` explicitly.
"""
macro gate(ex::Expr)
    @match ex begin
        :($locs => $gate) => esc(xcall(GlobalRef(Intrinsics, :apply), gate, xlocations(locs)))
        _ => error("syntax: invalid syntax, expect @gate <locs> => <gate>")
    end
end

"""
    @ctrl <ctrl_locs> <locs> => <gate>

Syntax sugar for `apply(gate, locs, ctrl_locs)`, must be used inside `@device`.
See also [`@device`](@ref).
"""
macro ctrl(ctrl_locs, ex::Expr)
    @match ex begin
        :($locs => $gate) => esc(
            xcall(GlobalRef(Intrinsics, :apply), gate, xlocations(locs), xctrl_locations(ctrl_locs)),
        )
        _ => error("syntax: invalid syntax, expect @ctrl <ctrl_locs> <locs> => <gate>")
    end
end

"""
    @measure locs

Syntax sugar for `measure(locs)`.
"""
macro measure(locs)
    esc(xcall(GlobalRef(Intrinsics, :measure), xlocations(locs)))
end

"""
    @barrier locs

Syntax sugar for `barrier(locs)`.
"""
macro barrier(locs)
    esc(xcall(GlobalRef(Intrinsics, :barrier), xlocations(locs)))
end

"""
    @device <function def>

Annotate a Julia function as YaoLang device kernel.
"""
macro device(ex)
    esc(device_m(__module__, ex))
end

"""
    routine_stub(routine, args...)

A function to store the raw `CodeInfo` of defined routines
to support multiple dispatch. We don't directly define `main`
`gate`, `ctrl` etc. because we can't support multiple dispatch
on them directly.
"""
function routine_stub end

function device_m(mod::Module, ex)
    is_function(ex) || error("expect a function definition")
    jlfn = JLFunction(ex)
    isnothing(jlfn.kwargs) || error("kwargs is not supported")

    return codegen_routine(jlfn)
end

function codegen_routine(jlfn::JLFunction)
    typename =
        isnothing(jlfn.name) ? gensym(:routine) :
        Meta.isexpr(jlfn.name, :(::)) ? jlfn.name.args[end] : Symbol("#", jlfn.name, "#")

    return quote
        $(codegen_routine_type(jlfn, typename))
        $(codegen_device(jlfn, typename))
        $(codegen_inference_limit_heuristics(jlfn, typename))
        $(codegen_routine_stub(jlfn, typename))
        $(codegen_binding(jlfn, typename))
    end
end

function codegen_routine_type(def::JLFunction, typename)
    def.name isa Symbol || isnothing(def.name) || return
    name = isnothing(def.name) ? typename : def.name

    jlstruct =
        JLStruct(; name = typename, supertype = :($YaoCompiler.GenericRoutine{$(QuoteNode(name))}))
    return codegen_ast(jlstruct)
end

function codegen_binding(def::JLFunction, typename)
    Meta.isexpr(def.name, :(::)) && return

    if isnothing(def.name)
        return :(Core.@__doc__ $typename())
    else
        return :(Core.@__doc__ const $(def.name) = $typename())
    end
end

function codegen_device(def::JLFunction, typename)
    self = gensym(:self)
    args = name_only.(def.args)

    jlfn = JLFunction(;
        name = :($self::$typename),
        args = def.args,
        whereparams = def.whereparams,
        rettype = def.rettype,
        line = def.line,
        body = quote
            $YaoCompiler.Operation($self, $(xtuple(args...)))
        end,
    )

    return codegen_ast(jlfn)
end

function inference_limit_heuristics end

function codegen_inference_limit_heuristics(def::JLFunction, typename)
    self = gensym(:self)
    args = name_only.(def.args)

    jlfn = JLFunction(;
        name = :($YaoCompiler.inference_limit_heuristics),
        args = [:($self::$typename), def.args...],
        whereparams = def.whereparams,
        rettype = def.rettype,
        line = def.line,
        body = quote
            $YaoCompiler.Operation($self, $(xtuple(args...)))
        end,
    )
    return codegen_ast(jlfn)
end

function codegen_routine_stub(def::JLFunction, typename)
    self = @match def.name begin
        :($name::$type) => name
        :(::$type) => gensym(:self)
        _ => gensym(:self)
    end

    jlfn = JLFunction(;
        name = :($YaoCompiler.routine_stub),
        args = [:($self::$typename), def.args...],
        whereparams = def.whereparams,
        rettype = def.rettype,
        line = def.line,
        body = transpile_gate_syntax(transpile_intrinsic(def.body)),
    )

    return codegen_ast(jlfn)
end

function is_syntax_macro(ex)
    @match ex begin
        Symbol("@gate") => true
        Symbol("@ctrl") => true
        Symbol("@measure") => true
        Symbol("@barrier") => true
        Expr(:., :YaoCompiler, QuoteNode(name)) => is_syntax_macro(name)
        Expr(:., YaoCompiler, QuoteNode(name)) => is_syntax_macro(name)
        GlobalRef(YaoCompiler, name) => is_syntax_macro(name)
        _ => false
    end
end

# NOTE: locs => gate is only a syntax sugar for gate(gate, locs)
function transpile_intrinsic(ex)
    @match ex begin
        Expr(:call, :apply, gate, locs) =>
            xcall(GlobalRef(Intrinsics, :apply), gate, xlocations(locs))
        Expr(:call, :apply, gate, locs, ctrl) =>
            xcall(GlobalRef(Intrinsics, :apply), gate, xlocations(locs), xctrl_locations(ctrl))
        Expr(:call, :measure, locs) => xcall(GlobalRef(Intrinsics, :measure), xlocations(locs))
        Expr(:call, :barrier, locs) => xcall(GlobalRef(Intrinsics, :barrier), xlocations(locs))
        Expr(:call, :expect, locs) => xcall(GlobalRef(Intrinsics, :expect), xlocations(locs))
        # check misused keyword error
        Expr(:call, :apply, _...) => error("syntax: apply is a preserved intrinsic function")
        Expr(:call, :measure, _...) => error("syntax: apply is a preserved intrinsic function")
        Expr(:call, :barrier, _...) => error("syntax: apply is a preserved intrinsic function")
        Expr(:call, :expect, _...) => error("syntax: apply is a preserved intrinsic function")
        Expr(head, args...) => Expr(head, map(transpile_intrinsic, args)...)
        _ => ex
    end
end

function transpile_gate_syntax(ex)
    @match ex begin
        # this only treat => syntax in block/let/if/for etc. as gate stmt
        :($locs => $gate) => xcall(GlobalRef(Intrinsics, :apply), gate, xlocations(locs))
        # this will appear in anonymous function definition
        # TODO: disambuigity this and function contains only single line
        # @device function circuit(theta, phi)
        #     1 => X
        # end
        # Expr(:block, stmt1, line::LineNumberNode, stmt2) => ex
        Expr(:macrocall, Symbol("@device"), _...) => error("syntax: cannot have nested @device")
        Expr(:macrocall, name, args...) => begin
            if is_syntax_macro(name)
                return ex
            else
                # we force top scope locs=>gate to be treated as gate stmt inside
                # all @device region including non-YaoCompiler macros to make things
                # like @inbounds etc. work
                return Expr(:macrocall, name, map(transpile_gate_syntax, args)...)
            end
        end
        # we only white list other syntax here to be safe
        Expr(:block, args...) ||
            Expr(:if, args...) ||
            Expr(:elseif, args...) ||
            Expr(:let, args...) ||
            Expr(:for, args...) ||
            Expr(:try, args...) => Expr(ex.head, map(transpile_gate_syntax, args)...)

        Expr(:function, call, body) => Expr(:function, call, transpile_gate_syntax(body))
        _ => ex
    end
end

@generated function Intrinsics.apply(::AbstractRegister, op::Operation{P,Args}) where {P,Args}
    ci, nargs = obtain_codeinfo(op)
    new = NewCodeInfo(ci)
    register = insert!(new.slots, 2, Symbol("#register#"))
    device = insert!(new.slots, 3, Symbol("#op#"))
    unpack_device!(new, device, nargs)
    for (v, stmt) in new
        @switch stmt begin
            @case Expr(:call, GlobalRef(&Intrinsics, name), args...)
            new[v] = xcall(GlobalRef(Intrinsics, name), register, args...)
            @case Expr(:(=), slot, Expr(:call, GlobalRef(&Intrinsics, name), args...))
            new[v] = Expr(:(=), slot, Expr(:call, GlobalRef(Intrinsics, name), register, args...))
            @case _
            nothing
        end
    end
    code_info = finish(new)
    if code_info.method_for_inference_limit_heuristics === nothing
        method = first(methods(inference_limit_heuristics, Tuple{P,Args.parameters...}))
        code_info.method_for_inference_limit_heuristics = method
    end
    return code_info
end

function _update_slot_stmt(f, new, v, stmt)
    @switch stmt begin
        @case Expr(:call, _...)
        new[v] = f(stmt)
        @case Expr(:(=), slot, rhs)
        new[v] = Expr(:(=), slot, f(rhs))
        @case _
        nothing
    end
end

@generated function Intrinsics.apply(
    ::AbstractRegister,
    op::Operation{P,Args},
    ::Locations,
) where {P,Args}
    ci, nargs = obtain_codeinfo(op)
    new = NewCodeInfo(ci)
    register = insert!(new.slots, 2, Symbol("#register#"))
    device = insert!(new.slots, 3, Symbol("#op#"))
    glob_locs = insert!(new.slots, 4, Symbol("#locs#"))
    unpack_device!(new, device, nargs)

    for (v, stmt) in new
        _update_slot_stmt(new, v, stmt) do stmt
            @switch stmt begin
                @case Expr(:call, GlobalRef(&Intrinsics, :apply), gate, locs)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                return xcall(GlobalRef(Intrinsics, :apply), register, gate, new_locs)

                @case Expr(:call, GlobalRef(&Intrinsics, :apply), gate, locs, ctrl)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                new_ctrl = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, ctrl))
                return xcall(GlobalRef(Intrinsics, :apply), register, gate, new_locs, new_ctrl)

                @case Expr(:call, GlobalRef(&Intrinsics, :measure), locs)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                return xcall(GlobalRef(Intrinsics, :measure), register, new_locs)

                @case Expr(:call, GlobalRef(&Intrinsics, :barrier), locs)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                return xcall(GlobalRef(Intrinsics, :barrier), register, new_locs)

                @case Expr(:call, GlobalRef(&Intrinsics, :expect), locs)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                return xcall(GlobalRef(Intrinsics, :expect), register, new_locs)
                @case _
                return stmt
            end
        end
    end
    code_info = finish(new)
    if code_info.method_for_inference_limit_heuristics === nothing
        method = first(methods(inference_limit_heuristics, Tuple{P,Args.parameters...}))
        code_info.method_for_inference_limit_heuristics = method
    end
    return code_info
end

@generated function Intrinsics.apply(
    ::AbstractRegister,
    op::Operation{P,Args},
    ::Locations,
    ::CtrlLocations,
) where {P,Args}
    ci, nargs = obtain_codeinfo(op)
    new = NewCodeInfo(ci)
    register = insert!(new.slots, 2, Symbol("#register#"))
    device = insert!(new.slots, 3, Symbol("#op#"))
    glob_locs = insert!(new.slots, 4, Symbol("#locs#"))
    glob_ctrl = insert!(new.slots, 5, Symbol("#ctrl#"))
    unpack_device!(new, device, nargs)

    for (v, stmt) in new
        _update_slot_stmt(new, v, stmt) do stmt
            @switch stmt begin
                @case Expr(:call, GlobalRef(&Intrinsics, :apply), gate, locs)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                return xcall(GlobalRef(Intrinsics, :apply), register, gate, new_locs, glob_ctrl)

                @case Expr(:call, GlobalRef(&Intrinsics, :apply), gate, locs, ctrl)
                new_locs = push!(new, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                new_ctrl = push!(new, xcall(GlobalRef(Base, :getindex), glob_locs, ctrl))
                new_ctrl =
                    push!(new, xcall(GlobalRef(YaoLocations, :merge_locations), new_ctrl, glob_ctrl))
                return xcall(GlobalRef(Intrinsics, :apply), register, gate, new_locs, new_ctrl)

                @case Expr(:call, GlobalRef(&Intrinsics, :measure), locs)
                return :(error("cannot apply quantum control on measurement"))

                @case Expr(:call, GlobalRef(&Intrinsics, :barrier), locs)
                new_locs = insert!(new, v, xcall(GlobalRef(Base, :getindex), glob_locs, locs))
                return xcall(GlobalRef(Intrinsics, :barrier), register, new_locs)

                @case Expr(:call, GlobalRef(&Intrinsics, :expect), locs)
                return :(error("cannot apply quantum control on measurement (expectation)"))
                @case _
                return stmt
            end
        end
    end
    code_info = finish(new)
    if code_info.method_for_inference_limit_heuristics === nothing
        method = first(methods(inference_limit_heuristics, Tuple{P,Args.parameters...}))
        code_info.method_for_inference_limit_heuristics = method
    end
    return code_info
end

function obtain_codeinfo(::Type{Operation{P,Args}}) where {P,Args}
    nargs = length(Args.parameters)
    tt = Tuple{P,Args.parameters...}
    ms = methods(routine_stub, tt)
    @assert length(ms) == 1
    method = first(ms)
    method_args = Tuple{typeof(routine_stub),tt.parameters...}
    mi = Core.Compiler.specialize_method(method, method_args, Core.svec())
    ci = Core.Compiler.retrieve_code_info(mi)

    name = routine_name(P)
    linetable = Any[]
    for lineinfo in ci.linetable
        push!(
            linetable,
            Core.LineInfoNode(
                lineinfo.module,
                name,
                lineinfo.file,
                lineinfo.line,
                lineinfo.inlined_at,
            ),
        )
    end
    ci.linetable = linetable
    ci.edges = Core.MethodInstance[mi]
    return ci, nargs
end

function unpack_device!(new::NewCodeInfo, op::NewSlotNumber, nargs::Int)
    # %parent = op.parent
    parent = push!(new, Expr(:call, GlobalRef(Base, :getfield), op, QuoteNode(:parent)))

    if nargs > 0
        # %args = op.args
        args = push!(new, Expr(:call, GlobalRef(Base, :getfield), op, QuoteNode(:args)))
    end

    # %self = %parent
    push!(new, Expr(:(=), SlotNumber(2), parent))

    for k in 1:nargs
        push!(new, Expr(:(=), SlotNumber(k + 2), xcall(GlobalRef(Base, :getindex), args, k)))
    end
    return new
end


text:
using Poly
using Test
using StaticArrays

@testset "correctness tests" begin
    @testset "basic non-macro test" begin
        # multiply arrays and double
        instructions = [Instruction(:mult, :(out[i, j] += A[i, k] * B[k, j]), Set([:i, :j, :k]), :()),
                        Instruction(:double, :(out[i, j] *= 2), Set([:i, :j]), :())]

        dom_k = Domain(:k, 1, :(r), 1, Set([:j]), [instructions[1]])
        dom_j = Domain(:j, 1, :(m), 1, Set([:i]), [dom_k, instructions[2]])
        dom_i = Domain(:i, 1, :(n), 1, Set(), [dom_j])
        domains = [dom_i, dom_j, dom_k]

        kern = compile(LoopKernel(instructions, domains, [:out, :A, :B, :n, :r, :m], [Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}], []))

        A = rand(10, 10)
        B = rand(10, 10)
        out = zeros(10, 10)
        n = size(out, 1)
        r = size(A, 2)
        m = size(out, 2)

        kern(out=out, A=A, B=B, n=n, r=r, m=m)

        @test isapprox(out, A*B*2)
    end

    @testset "basic non-macro test 2" begin
        # add arrays with global offset
        offset = 4
        instructions = [Instruction(:add, :(out[i, j] = A[i, j] + B[i, j] + offset), Set([:i, :j]), :())]

        dom_j = Domain(:j, 1, :(m), 1, Set([:i]), [instructions[1]])
        dom_i = Domain(:i, 1, :(n), 1, Set(), [dom_j])
        domains = [dom_i, dom_j]

        kern = compile(LoopKernel(instructions, domains, [:out, :A, :B, :n, :m, :offset], [Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}], []))

        A = rand(10, 10)
        B = rand(10, 10)
        out = zeros(10, 10)
        n = size(out, 1)
        m = size(out, 2)

        kern(out=out, A=A, B=B, n=n, m=m, offset=offset)

        @test isapprox(out, A+B.+4)
    end

    @testset "simple macro test" begin
        # test macros simple
        A = rand(10)
        out = zeros(10)
        n = size(out, 1)

        @poly_loop for i = 1:n
            out[i] = A[i]*2
        end

        @test isapprox(out, A*2)
    end

    @testset "nested macro test" begin
        # test macros nested loops
        A = rand(10, 10)
        B = rand(10, 10)
        out = zeros(10, 10)
        n = size(out, 1)
        m = size(out, 2)
        r = size(A, 2)

        @poly_loop for i = 1:n
            for j = 1:m
                for k = 1:r
                    out[i, j] += A[i, k] * B[k, j]
                end
                out[i, j] *= 2
            end
        end

        @test isapprox(out, A*B*2)
    end

    @testset "macro stride test" begin
        # test loops with non-1 stride
        arr = ones(6)
        @poly_loop for i = 1:3:6
            arr[i] += 1
        end
        @test arr == [2, 1, 1, 2, 1, 1]
    end

    @testset "macro stride 2D test" begin
        # test loops with non-1 stride
        arr = ones(6, 6)
        @poly_loop for i = 1:2:6
            for j = 1:2:6
                arr[i, j] += 1
            end
        end

        arr2 = ones(6, 6)
        for i = 1:2:6
            for j = 1:2:6
                arr2[i, j] += 1
            end
        end

        @test arr == arr2
    end

    @testset "inter-loop dependence test" begin
        # test loop iterators depending on other iterators
        arr = zeros(3, 3)
        expected = [1 1 1;
                    0 1 1;
                    0 0 1]
        @poly_loop for i = 1:3
            for j = i:3
                arr[i, j] = 1
            end
        end

        @test arr == expected
    end

    @testset "inter-loop dependence test 2" begin
        # test loop iterators depending on other iterators
        arr = zeros(3, 3)
        expected = [1 0 0;
                    1 1 0;
                    1 1 1]
        @poly_loop for i = 1:3
            for j = 1:i
                arr[i, j] = 1
            end
        end

        @test arr == expected
    end

    @testset "inter-loop dependence test 3" begin
        # test loop iterators depending on other iterators
        arr = zeros(3, 3)
        expected = [6 6 6;
                    17 5 5;
                    9 8 3]
        @poly_loop for i = 1:3
            for j = 1:3
                for k = i:3
                    arr[i, j] += arr[j, k] + k
                end
            end
        end

        @test arr == expected
    end

    @testset "dependence test" begin
        # simple @depends_on test
        count = 0
        @poly_loop for i = 1:10
            count += 1
        end

        @test count == 10
    end

    @testset "min/max test" begin
        # checks that min and max functions work
        n = 10
        m = 12
        arr = zeros(n, m)

        @poly_loop for i=1:max(n, 8)
            for j=1:min(15, m)
                arr[i, j] = 1
            end
        end

        arr2 = zeros(n, m)

        for i=1:max(n, 8)
            for j=1:min(15, m)
                arr2[i, j] = 1
            end
        end

        @test arr == arr2
    end

    @testset "floor test" begin
        # tests that floor function works
        n = 10
        arr = zeros(n, n)

        @poly_loop for i=1:floor(n/2)
            for j=1:floor(n/2)
                arr[i, j] = 1
            end
        end

        arr2 = zeros(n, n)

        for i=1:Int(floor(n/2))
            for j=1:Int(floor(n/2))
                arr2[i, j] = 1
            end
        end

        @test arr == arr2
    end

    @testset "threading test" begin
        n = 128
        r = 128
        m = 128
        A = rand(n, r)
        B = rand(r, m)
        out = zeros(n, m)
        n = size(out, 1)
        m = size(out, 2)
        r = size(A, 2)

        @poly_loop thread=true for i = 1:n
            for j = 1:m
                for k = 1:r
                    out[i, j] += A[i, k] * B[k, j]
                end
            end
        end

        @test isapprox(out, A*B)
    end

    @testset "if block test" begin
        # tests that an if block works
        n = 10
        arr = zeros(n)

        @poly_loop for i=1:n
            if i % 2 == 0
                arr[i] = 1
            end
        end

        arr2 = zeros(n)

        for i=1:n
            if i % 2 == 0
                arr2[i] = 1
            end
        end

        @test arr == arr2
    end

    @testset "if/elseif block test" begin
        # tests that an if/elseif block works
        n = 10
        arr = zeros(n, n)

        @poly_loop for i=1:n
            for j=1:n
                if i % 2 == 0
                    arr[i, j] = 1
                elseif j % 2 == 0
                    arr[i, j] = 2
                end
            end
        end

        arr2 = zeros(n, n)

        for i=1:n
            for j=1:n
                if i % 2 == 0
                    arr2[i, j] = 1
                elseif j % 2 == 0
                    arr2[i, j] = 2
                end
            end
        end

        @test arr == arr2
    end

    @testset "if/elseif/else block test" begin
        # tests that an if/elseif/else block works
        n = 10
        arr = zeros(n, n)

        @poly_loop for i=1:n
            for j=1:n
                if i % 2 == 0
                    arr[i, j] = 1
                elseif (j + n) % 2 == 0
                    arr[i, j] = 2
                else
                    arr[i, j] = 3
                end
            end
        end

        arr2 = zeros(n, n)

        for i=1:n
            for j=1:n
                if i % 2 == 0
                    arr2[i, j] = 1
                elseif (j + n) % 2 == 0
                    arr2[i, j] = 2
                else
                    arr2[i, j] = 3
                end
            end
        end

        @test arr == arr2
    end

    @testset "if/else block test" begin
        # tests that an if/else block works
        n = 10
        arr = zeros(n, n)

        @poly_loop for i=1:n
            for j=1:n
                if n % 2 == 0
                    arr[i, j] = 1
                else
                    arr[i, j] = 2
                end
            end
        end

        arr2 = zeros(n, n)

        for i=1:n
            for j=1:n
                if n % 2 == 0
                    arr2[i, j] = 1
                else
                    arr2[i, j] = 2
                end
            end
        end

        @test arr == arr2
    end

    @testset "interpolation test 1" begin
        # checks that interpolation works
        c = 2
        n = 10
        arr = zeros(n)

        @poly_loop for i=1:$c:n
            arr[i] = 1
        end

        arr2 = zeros(n)
        for i=1:c:n
            arr2[i] = 1
        end

        @test arr == arr2
    end

    @testset "interpolation test 2" begin
        # checks that multiple 2d interpolation works
        c = 2
        d = 3
        n = 10
        m = 9
        arr = zeros(n, m)

        @poly_loop for i=1:$c:$n
            for j=1:$d:$m
                arr[i, j] = 1
            end
        end

        arr2 = zeros(n, m)
        for i=1:c:n
            for j=1:d:m
                arr2[i, j] = 1
            end
        end

        @test arr == arr2
    end

    @testset "non-uniform dependence" begin
        n = 64
        arr = ones(n, n)
        @poly_loop tile=0 for i = 1:n
            for j = 3:n
                arr[i, j] = arr[j, i] + arr[i, j-1]
            end
        end

        arr2 = ones(n, n)
        for i = 1:n
            for j = 3:n
                arr2[i, j] = arr2[j, i] + arr2[i, j-1]
            end
        end

        @test isapprox(arr, arr2)
    end

    @testset "basic matmul" begin
        # test macros basic matmul
        n = 128
        r = 128
        m = 128
        A = rand(n, r)
        B = rand(r, m)
        out = zeros(n, m)
        n = size(out, 1)
        m = size(out, 2)
        r = size(A, 2)

        @poly_loop for i = 1:n
            for j = 1:m
                for k = 1:r
                    out[i, j] += A[i, k] * B[k, j]
                end
            end
        end

        @test isapprox(out, A*B)
    end

    @testset "1d jacobi" begin
        N = 50
        T = 10
        a = rand(N)
        aa = copy(a)
        b = zeros(N)

        @poly_loop tile=0 for t = 1:T-1
            for i = 3:N-2
                b[i] = (a[i - 1] + a[i] + a[i + 1])/3
            end
            for ii = 3:N-2
                a[ii] = b[ii]
            end
        end

        bb = zeros(N)
        for t = 1:T-1
            for i = 3:N-2
                bb[i] = (aa[i - 1] + aa[i] + aa[i + 1])/3
            end
            for ii = 3:N-2
                aa[ii] = bb[ii]
            end
        end

        @test isapprox(a, aa)
    end

    @testset "lu decomposition" begin
        n = 128
        A = rand(n, n)*10
        L = zeros(n, n)
        U = zeros(n, n)
        P = zeros(n, n)
        for i = 1:n
            P[i, i] = 1.0
        end
        for j=1:n
            i = findmax(A[:, j])[2]
            if j != i
                P[j, :], P[i, :] = P[i, :], P[j, :]
            end
        end
        PA = P*A

        @poly_loop for j=1:n
            L[j, j] = 1.0
            for i=1:j
                s1 = 0.0
                for k=1:i
                    s1 += U[k, j] * L[i, k]
                end
                U[i, j] = PA[i, j] - s1
            end

            for ii=j+1:n
                s2 = 0.0
                for kk=1:j
                    s2 += U[kk, j] * L[ii, kk]
                end
                L[ii, j] = (PA[ii, j] - s2) / U[j, j]
            end
        end

        @test isapprox(PA, L*U)
    end

    @testset "tiled matrix multiplication" begin
        # test macros complicated (tiled matrix multiplication)
        N = 128

        A = rand(N, N)
        B = rand(N, N)
        C = zeros(N, N)

        TILE_DIM = 32
        tile1 = @MArray zeros(TILE_DIM, TILE_DIM)
        tile2 = @MArray zeros(TILE_DIM, TILE_DIM)

        @poly_loop tile=0 for gj = 0:$TILE_DIM:N-1
            for gi = 0:$TILE_DIM:N-1
                for t = 0:$TILE_DIM:N-1
                    for j = 1:$TILE_DIM
                        for i = 1:$TILE_DIM
                            tile1[i, j] = A[gi + i, t + j]
                            tile2[i, j] = B[t + i, gj + j]
                        end
                    end
                    for jj in 1:$TILE_DIM
                        for k = 1:$TILE_DIM
                            for ii = 1:$TILE_DIM
                                C[gi + ii, gj + jj] += tile1[ii, k] * tile2[k, jj]
                            end
                        end
                    end
                end
            end
        end

        @test isapprox(C, A*B)
    end

end

# uncomment to run performance testing (takes a while, so commented out for github)
# include("performancetests.jl")

ext: jl


text:
@with_kw mutable struct RandomSearchSolver
    n_iterations::Int64 = 100
    episode_length::Int64 = 10
    show_progress::Bool = true
end

mutable struct RandomSearchPlanner{P<:Union{MDP,POMDP}}
    solver::RandomSearchSolver
    mdp::P
end

# No work performed.
POMDPs.solve(solver::RandomSearchSolver, mdp::Union{POMDP,MDP}) = RandomSearchPlanner(solver, mdp)

# Online work performed here.
function POMDPs.action(planner::RandomSearchPlanner, s)
    # Pull out variables
    mdp::ASTMDP = planner.mdp
    sim::GrayBox.Simulation = mdp.sim
    n_iterations::Int64 = planner.solver.n_iterations
    depth::Int64 = planner.solver.episode_length
    planner.solver.show_progress ? progress = Progress(n_iterations) : nothing

    # Run simulation
    for i in 1:n_iterations
        BlackBox.initialize!(sim)
        AST.go_to_state(mdp, s)
        AST.rollout(mdp, s, depth) # Evaluate system.
        planner.solver.show_progress ? next!(progress) : nothing
    end

    # Pass back action trace if recording is on (i.e. top_k)
    if mdp.params.top_k > 0
        return get_top_path(mdp)
    end
end

# Search using the planner from the initial AST state.
# Pass back best action trace.
function AST.search!(planner::RandomSearchPlanner)
    mdp::ASTMDP = planner.mdp
    Random.seed!(mdp.params.seed) # Determinism
    s = AST.initialstate(mdp)
    return action(planner, s)
end


text:
struct VenusType <: abstractPlanet end
Venus = Planet(VenusType)

#Data, as previously used in AstroDynBase
mu(::Type{VenusType}) = 3.24859e5
j2(::Type{VenusType}) = 6e-5
mean_radius(::Type{VenusType}) = 6051.8
equatorial_radius(::Type{VenusType}) = 6051.8
polar_radius(::Type{VenusType}) = 6051.8

naif_id(::Type{VenusType}) = 299
parent(::Type{VenusType}) = VenusBarycenter
frame(::Type{VenusType}) = "IAU_VENUS"


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
Table

The object containing the flatbuffer and positional information specific to the table.
The `vtable` containing the offsets for specific members precedes `pos`.
The actual values in the table follow `pos` offset and size of the vtable.

- `bytes::Vector{UInt8}`: the flatbuffer itself
- `pos::Integer`:  the base position in `bytes` of the table
"""
abstract type Table end
abstract type Struct end

const TableOrStruct = Union{Table, Struct}

bytes(x::TableOrStruct) = getfield(x, :bytes)
pos(x::TableOrStruct) = getfield(x, :pos)

(::Type{T})(b::Builder) where {T <: TableOrStruct} = T(b.bytes[b.head+1:end], get(b, b.head, Int32))

getrootas(::Type{T}, bytes::Vector{UInt8}, offset) where {T <: Table} = init(T, bytes, offset + readbuffer(bytes, offset, UOffsetT))
init(::Type{T}, bytes::Vector{UInt8}, pos::Integer) where {T <: TableOrStruct} = T(bytes, pos)

const TableOrBuilder = Union{Table, Struct, Builder}

Base.get(t::TableOrBuilder, pos, ::Type{T}) where {T} = readbuffer(bytes(t), pos, T)
Base.get(t::TableOrBuilder, pos, ::Type{T}) where {T <: Enum} = T(get(t, pos, basetype(T)))

"""
`offset` provides access into the Table's vtable.

Deprecated fields are ignored by checking against the vtable's length.
"""
function offset(t::Table, vtableoffset)
    vtable = pos(t) - get(t, pos(t), SOffsetT)
    return vtableoffset < get(t, vtable, VOffsetT) ? get(t, vtable + vtableoffset, VOffsetT) : VOffsetT(0)
end

"`indirect` retrieves the relative offset stored at `offset`."
indirect(t::Table, off) = off + get(t, off, UOffsetT)

getvalue(t, o, ::Type{Nothing}) = nothing
getvalue(t, o, ::Type{T}) where {T <: Scalar} = get(t, pos(t) + o, T)
getvalue(t, o, ::Type{T}) where {T <: Enum} = T(get(t, pos(t) + o, enumtype(T)))

function Base.String(t::Table, off)
    off += get(t, off, UOffsetT)
    start = off + sizeof(UOffsetT)
    len = get(t, off, UOffsetT)
    return unsafe_string(pointer(bytes(t), start + 1), len)
end

function bytevector(t::Table, off)
    off += get(t, off, UOffsetT)
    start = off + sizeof(UOffsetT)
    len = get(t, off, UOffsetT)
    return view(bytes(t), (start + 1):(start + len + 1))
end

"""
`vectorlen` retrieves the length of the vector whose offset is stored at
`off` in this object.
"""
function vectorlen(t::Table, off)
    off += pos(t)
    off += get(t, off, UOffsetT)
    return Int(get(t, off, UOffsetT))
end

"""
`vector` retrieves the start of data of the vector whose offset is stored
at `off` in this object.
"""
function vector(t::Table, off)
    off += pos(t)
    x = off + get(t, off, UOffsetT)
    # data starts after metadata containing the vector length
    return x + sizeof(UOffsetT)
end

struct Array{T, S, TT} <: AbstractVector{T}
    _tab::TT
    pos::Int64
    data::Vector{S}
end

function Array{T}(t::Table, off) where {T}
    a = vector(t, off)
    S = T <: Table ? UOffsetT : T <: Struct ? NTuple{structsizeof(T), UInt8} : T
    ptr = convert(Ptr{S}, pointer(bytes(t), a + 1))

    if S === T # native type
        return reinterpret(S, bytes(t)[a+1:a+sizeof(S)*vectorlen(t, off)])
    else
        data = unsafe_wrap(Base.Array, ptr, vectorlen(t, off))
        return Array{T, S, typeof(t)}(t, a, data)
    end
end

function structsizeof end

Base.IndexStyle(::Type{<:Array}) = Base.IndexLinear()
Base.size(x::Array) = size(x.data)
Base.@propagate_inbounds function Base.getindex(A::Array{T, S}, i::Integer) where {T, S}
    if T === S
        return A.data[i]
    elseif T <: Struct
        return init(T, bytes(A._tab), A.pos + (i - 1) * structsizeof(T))
    else # T isa Table
        return init(T, bytes(A._tab), indirect(A._tab, A.pos + (i - 1) * 4))
    end
end

Base.@propagate_inbounds function Base.setindex!(A::Array{T, S}, v, i::Integer) where {T, S}
    if T === S
        return setindex!(A.data, v, i)
    else
        error("setindex! not supported for reference/table types")
    end
end

function union(t::Table, off)
    off += pos(t)
    return off + get(t, off, UOffsetT)
end

function union!(t::Table, t2::Table, off)
    off += pos(t)
    t2.pos = off + get(t, off, UOffsetT)
    t2.bytes = bytes(t)
    return
end

"""
GetVOffsetTSlot retrieves the VOffsetT that the given vtable location
points to. If the vtable value is zero, the default value `d`
will be returned.
"""
function getoffsetslot(t::Table, slot, d)
    off = offset(t, slot)
    if off == 0
        return d
    end
    return off
end

"""
`getslot` retrieves the `T` that the given vtable location
points to. If the vtable value is zero, the default value `d`
will be returned.
"""
function getslot(t::Table, slot, d::T) where {T}
    off = offset(t, slot)
    if off == 0
        return d
    end

    return get(t, pos(t) + off, T)
end


text:
export SimpleObservable, stddev

type SimpleObservable <: ScalarObservable
  bins :: Vector{Float64}
  num :: Int64
  sum :: Float64
  sum2 :: Float64
end

SimpleObservable() = SimpleObservable(zeros(0), 0, 0.0, 0.0)

function reset!(obs :: SimpleObservable)
  obs.bins = zeros(0)
  obs.num = 0
  obs.sum = 0.0
  obs.sum2 = 0.0
  return obs
end

count(obs::SimpleObservable) = obs.num

function push!(obs :: SimpleObservable, value)
  push!(obs.bins, value)
  obs.num += 1
  obs.sum += value
  obs.sum2 += value^2
  return obs
end

function mean(obs::SimpleObservable)
  if obs.num > 0
    return obs.sum / obs.num
  else
    return NaN
  end
end

function var(obs::SimpleObservable)
  if obs.num  > 1
    v = (obs.sum2 - obs.sum*obs.sum/obs.num)/(obs.num-1)
    return maxzero(v)
  elseif obs.num == 1
    return Inf
  else
    return NaN
  end
end
stddev(obs::SimpleObservable) = sqrt(var(obs))
stderror(obs::SimpleObservable) = sqrt(var(obs)/count(obs))
function confidence_interval(obs::SimpleObservable, confidence_rate :: Real)
  if count(obs) < 2
    return Inf
  end
  q = 0.5 + 0.5confidence_rate
  correction = quantile( TDist(obs.num - 1), q)
  serr = stderror(obs)
  return correction * serr
end

function confidence_interval(obs::SimpleObservable, confidence_rate_symbol::Symbol = :sigma1)
  n = parsesigma(confidence_rate_symbol)
  return confidence_interval(obs, erf(0.5n*sqrt(2.0)))
end

function merge!(obs::SimpleObservable, other::SimpleObservable)
  append!(obs.bins, other.bins)
  obs.num += other.num
  obs.sum += other.sum
  obs.sum2 += other.sum2
  return obs
end
merge(lhs::SimpleObservable, rhs::SimpleObservable) = merge!(deepcopy(lhs), rhs)

export SimpleObservableSet
typealias SimpleObservableSet MCObservableSet{SimpleObservable}

function merge!(obs::SimpleObservableSet, other::SimpleObservableSet)
  obs_names = Set(keys(obs))
  union!(obs_names, Set(keys(other)))
  for name in obs_names
    if !haskey(obs, name)
      # in other only
      obs[name] = deepcopy(other[name])
    elseif haskey(other, name)
      # in both
      merge!(obs[name], other[name])

    # else
    # in obs only
    # NOTHING to do
    end
  end
  return obs
end
merge(lhs::SimpleObservableSet, rhs::SimpleObservableSet) = merge!(deepcopy(lhs), rhs)


text:
using LinearAlgebra
include("Pauli.jl")

function BlochRotation(axis,phase)

#Axis is the axis of rotation, should be a 1x3 array. It is not normalized,
#therefore being axis=angle*axis
#phase is the general phase

p=[Pauli(1), Pauli(2), Pauli(3)]
a=axis

exp(-phase*im)exp(-(p[1]*axis[1]+p[2]*a[2]+p[3]*a[3])/2*im)

end


text:
using BCBIStats
using Test

my_tests = ["cooccur.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * %s\n", my_test)
    include(my_test)
end

max_stars_repo_name: bcbi/BCBIStats

max_issues_repo_name: bcbi/BCBIStats

max_forks_repo_name: bcbi/BCBIStats


text:
"""
    leftjoin!(df1, df2; on, makeunique=false, source=nothing,
              matchmissing=:error)

Perform a left join of two data frame objects by updating the `df1` with the
joined columns from `df2`.

A left join includes all rows from `df1` and leaves all rows and columns from `df1`
untouched. Note that each row in `df1` must have at most one match in `df2`. Otherwise,
this function would not be able to execute the join in-place since new rows would need to be
added to `df1`.

# Arguments
- `df1`, `df2`: the `AbstractDataFrames` to be joined

# Keyword Arguments
- `on` : A column name to join `df1` and `df2` on. If the columns on which
  `df1` and `df2` will be joined have different names, then a `left=>right`
  pair can be passed. It is also allowed to perform a join on multiple columns,
  in which case a vector of column names or column name pairs can be passed
  (mixing names and pairs is allowed).
- `makeunique` : if `false` (the default), an error will be raised
  if duplicate names are found in columns not joined on;
  if `true`, duplicate names will be suffixed with `_i`
  (`i` starting at 1 for the first duplicate).
- `source` : Default: `nothing`. If a `Symbol` or string, adds indicator
  column with the given name, for whether a row appeared in only `df1` (`"left_only"`)
  or in both (`"both"`). If the name is already in use,
  the column name will be modified if `makeunique=true`.
- `matchmissing` : if equal to `:error` throw an error if `missing` is present
  in `on` columns; if equal to `:equal` then `missing` is allowed and missings are
  matched; if equal to `:notequal` then missings are dropped in `df2` `on` columns;
  `isequal` is used for comparisons of rows for equality

The columns added to `df1` from `df2` will support missing values.

It is not allowed to join on columns that contain `NaN` or `-0.0` in real or
imaginary part of the number. If you need to perform a join on such values use
CategoricalArrays.jl and transform a column containing such values into a
`CategoricalVector`.

See also: [`leftjoin`](@ref).

# Examples
```jldoctest
julia> name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
3×2 DataFrame
 Row │ ID     Name
     │ Int64  String
─────┼──────────────────
   1 │     1  John Doe
   2 │     2  Jane Doe
   3 │     3  Joe Blogs

julia> job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ ID     Job
     │ Int64  String
─────┼───────────────
   1 │     1  Lawyer
   2 │     2  Doctor
   3 │     4  Farmer

julia> leftjoin!(name, job, on = :ID)
3×3 DataFrame
 Row │ ID     Name       Job
     │ Int64  String     String?
─────┼───────────────────────────
   1 │     1  John Doe   Lawyer
   2 │     2  Jane Doe   Doctor
   3 │     3  Joe Blogs  missing

julia> job2 = DataFrame(identifier = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])
3×2 DataFrame
 Row │ identifier  Job
     │ Int64       String
─────┼────────────────────
   1 │          1  Lawyer
   2 │          2  Doctor
   3 │          4  Farmer

julia> leftjoin!(name, job2, on = :ID => :identifier, makeunique=true, source=:source)
3×5 DataFrame
 Row │ ID     Name       Job      Job_1    source
     │ Int64  String     String?  String?  String
─────┼───────────────────────────────────────────────
   1 │     1  John Doe   Lawyer   Lawyer   both
   2 │     2  Jane Doe   Doctor   Doctor   both
   3 │     3  Joe Blogs  missing  missing  left_only
```
"""
function leftjoin!(df1::AbstractDataFrame, df2::AbstractDataFrame;
                   on::Union{<:OnType, AbstractVector}=Symbol[], makeunique::Bool=false,
                   source::Union{Nothing, Symbol, AbstractString}=nothing,
                   matchmissing::Symbol=:error)

    _check_consistency(df1)
    _check_consistency(df2)

    if !is_column_insertion_allowed(df1)
        throw(ArgumentError("leftjoin! is only supported if `df1` is a `DataFrame`, " *
                            "or a SubDataFrame created with `:` as column selector"))
    end

    if on == []
        throw(ArgumentError("Missing join argument 'on'."))
    end

    joiner = DataFrameJoiner(df1, df2, on, matchmissing, :left)

    right_noon_names = names(joiner.dfr, Not(joiner.right_on))
    if !(makeunique || isempty(intersect(right_noon_names, names(df1))))
        throw(ArgumentError("the following columns are present in both " *
                            "left and right data frames but not listed in `on`: " *
                            join(intersect(right_noon_names, names(df1)), ", ") *
                            ". Pass makeunique=true to add a suffix automatically to " *
                            "columns names from the right data frame."))
    end

    left_ixs_inner, right_ixs_inner = find_inner_rows(joiner)

    right_ixs = _map_leftjoin_ixs(nrow(df1), left_ixs_inner, right_ixs_inner)

    # TODO: consider adding threading support in the future
    for colname in right_noon_names
        rcol = joiner.dfr[!, colname] # note that joiner.dfr does not have to be df2
        rcol_joined = compose_joined_rcol!(rcol, similar_missing(rcol, nrow(df1)),
                                          right_ixs)
        # if df1 isa SubDataFrame we must copy columns
        insertcols!(df1, colname => rcol_joined, makeunique=makeunique,
                    copycols=!(df1 isa DataFrame))
    end

    if source !== nothing
        pool = ["left_only", "right_only", "both"]
        invpool = Dict{String, UInt32}("left_only" => 1,
                                       "right_only" => 2,
                                       "both" => 3)
        indicatorcol = PooledArray(PooledArrays.RefArray(UInt32.(2 .* (right_ixs .> 0) .+ 1)),
                                   invpool, pool)

        unique_indicator = source
        if makeunique
            try_idx = 0
            while hasproperty(df1, unique_indicator)
                try_idx += 1
                unique_indicator = Symbol(source, "_", try_idx)
            end
        end

        if hasproperty(df1, unique_indicator)
            throw(ArgumentError("joined data frame already has column " *
                                ":$unique_indicator. Pass makeunique=true to " *
                                "make it unique using a suffix automatically."))
        end
        df1[!, unique_indicator] = indicatorcol
    end
    return df1
end

function _map_leftjoin_ixs(out_len::Int,
                           left_ixs_inner::Vector{Int},
                           right_ixs_inner::Vector{Int})
    right_ixs = zeros(Int, out_len)
    @inbounds for (li, ri) in zip(left_ixs_inner, right_ixs_inner)
        if right_ixs[li] > 0
            throw(ArgumentError("duplicate rows found in right table"))
        end
        right_ixs[li] = ri
    end
    return right_ixs
end

function compose_joined_rcol!(rcol::AbstractVector,
                              rcol_joined::AbstractVector,
                              right_ixs::Vector{Int})
    @assert length(rcol_joined) == length(right_ixs)
    @inbounds for (i, idx) in enumerate(right_ixs)
        if idx > 0
            rcol_joined[i] = rcol[idx]
        end
    end
    return rcol_joined
end


text:
"""
    timenormalize(data, events::AbstractVector{Int}[, len=100])

Normalize the first dimension to lengths of `len` bounded by the events.
"""
function timenormalize(data::AbstractArray{T}, events::AbstractVector{Int}, len::Int=100) where T
    dims = size(data)
    oobevents = events .∈ Ref(axes(data,1))
    any(~, oobevents) &&
        throw(error("events $(events[.~oobevents]) are not valid indices of `data`"))
    res = Array{T}(undef, (length(events)-1)*len, dims[2:end]...)

    fill_normdims!(res, data, events, len)

    return res
end

function fill_normdims!(res, data, events::AbstractVector{Int}, len::Int=100)
    (size(res,2) != size(data,2)) && throw(
        ArgumentError("`res` and `data` must have the same number of columns"))
    for i in 1:size(data,2)
        # Create interpolation object
        itp = interpolate(@view(data[:,i]), BSpline(Cubic(Line(Interpolations.OnGrid()))))

        fill_normstrides!(@view(res[:,i]), itp, events, len)
    end

    nothing
end

function fill_normstrides!(nstr, str, events::AbstractVector{Int}, len::Int=100)
    @assert mod(size(nstr,1),len) === 0
    @assert size(nstr,1) === (length(events)-1)*len

    # Create time normalized strides
    @inbounds for s in 1:length(events)-1
        nstr[(1:len).+(s-1)*len] = str(normtime(events[s], events[s+1], len))
    end

    nothing
end

"""
    timestoindices(t1::Real, t2::Real, times[, len=100])

Find the equivalent indices for `times` when the times `t1` and `t2` match the `begin` and
`end + 1` of an array of length `len`.

This is useful when an array has been timenormalized using `basetimes` as the normalizing events, but when the location of `times` in the normalized array is also needed.

# Examples
```jldoctest
julia> normedarray = timenormalize(collect(.01:.01:3), [132, 246])
100-element Vector{Float64}:
 1.32
 1.3314
 1.3428
 1.3541999999999996
[...]

julia> times = [1.56, 2.0]
2-element Vector{Float64}:
 1.56
 2.0

julia> idx = timestoindices(1.32, 2.46, times)
2-element Vector{Int64}:
 23
 61

julia> isapprox(normedarray[idx], times; atol=.02)
true
```
"""
function timestoindices(t1::T, t2::T, events, len::Int=100) where T <: Real
    t1 < t2 || throw(DomainError("t1 must be less than t2"))
    all(t1 .<= events .<= t2) || throw(DomainError("events must be between normalizing events"))

    nt = normtime(t1, t2, len)
    return [ searchsortedfirst(nt, ev) for ev in events ]
end

"""
    timestoindices(basetimes::AbstractVector{<:Real}, times[, len=100])

Find the equivalent indices for elements in `times` which are between `basetimes[i]` and
`basetimes[i+1]` which match the `begin` and `end + 1` of an array of length `len` for all successive pairs in `basetimes`.
"""
function timestoindices(basetimes::AbstractVector{T}, times::AbstractVector{T}, len::Int=100) where T
    ntimes = similar(times)
    for i in 1:(length(basetimes)-1)
        rel = findall(x -> basetimes[i] ≤ x ≤ basetimes[i+1], times)
        @views ntimes[rel] .= timestoindices(basetimes[i], basetimes[i+1], times[rel], len)
    end
    return ntimes
end

"""
    normtime(t1, t2[, len=100])

Returns a range [t1,t2) with a length of `len`.
"""
function normtime(t1, t2, len::Int=100)
    # st = (t2-t1)/len
    # t1:st:round(T, t2-st)
    range(t1, stop=t2, length=len+1)[1:len]
end

"""
    limitcycle(data, [events::AbstractVector{Int}, len=100])

Calculate the ensemble average and stdev of `data` with a normalized length of `len`
"""
function limitcycle(data, len::Int=100; mean=mean, std=std)
    N = size(data, 1)
    mod(N, len) == 0 ||
        throw(ArgumentError("length of data must be even multiple of len"))

    wid = size(data, 2)
    ensemble_avg = similar(data, (len, wid))
    ensemble_std = similar(data, (len, wid))

    @inbounds @views for dim in 1:wid
        for i in 1:len
            ensemble_avg[i,dim] = mean(data[i:len:N, dim])
            ensemble_std[i,dim] = std(data[i:len:N, dim])
        end
    end

    return (ensemble_avg, ensemble_std)
end

function limitcycle(data, events::AbstractVector{Int}, len::Int=100; mean=mean, std=std)
    normed = timenormalize(data, events, len)
    limitcycle(normed, len; mean, std)
end


text:
@testset "plots.jl" begin
    n = 5000

    y = rand(GeneralizedExtremeValue(1.0, 1.0, 0.1), n)
    fm = MaximumLikelihoodEVA(BlockMaxima(Variable("y", y)), [1.0, 0.0, 0.1])

    x = rand(n)
    yns = rand.(GeneralizedExtremeValue.(x, 1.0, 0.1))
    fmns = MaximumLikelihoodEVA(BlockMaxima(Variable("y", yns), locationcov = [Variable("x", x)]), [0.0, 1.0, 0.0, 0.1])

    @testset "probplot_data(fm)" begin
        df = probplot_data(fm)

        # Returns a dataframe with n values in column Model
        @test length(df[:, :Model]) == n

        # Returns a dataframe with n values in column Empirical
        @test length(df[:, :Empirical]) == n

    end

    @testset "probplot(fm)" begin
        # Plot does not throw
        @test_logs probplot(fm)

    end

    @testset "qqplot_data(fm)" begin
        df = qqplot_data(fm)

        # Returns a dataframe with n values in column Model
        @test length(df[:, :Model]) == n

        # Returns a dataframe with n values in column Empirical
        @test length(df[:, :Empirical]) == n

        # Info for non-stationary model but does not throw
        @test_logs (:info, "The graph is optimized for stationary models and the model provided is not.") qqplot_data(fmns)

    end

    @testset " qqplot(fm)" begin
        # Plot does not throw
        @test_logs qqplot(fm)

    end

    @testset "returnlevelplot_data(fm)" begin
        df = returnlevelplot_data(fm)

        # Returns a dataframe with n values in column Data
        @test length(df[:, :Data]) == n

        # Returns a dataframe with n values in column Period
        @test length(df[:, :Period]) == n

        # Returns a dataframe with n values in column Level
        @test length(df[:, :Level]) == n

    end

    @testset "returnlevelplot(fm)" begin
        # Plot does not throw
        @test_logs returnlevelplot(fm)

    end

    @testset "histplot_data(fm)" begin
        dfs = histplot_data(fm)

        # Returns a dataframe with n values in column Data
        @test length(dfs[:h][:, :Data]) == n

        # Returns a dataframe with 1000 values in column DataRange
        @test length(dfs[:d][:, :DataRange]) == 1000

        # Returns a dataframe with 1000 values in column Density
        @test length(dfs[:d][:, :Density]) == 1000

        # Returns a dictionary with key nbin
        @test haskey(dfs, :nbin)

        # Returns a dictionary with key xmin
        @test haskey(dfs, :xmin)

        # Returns a dictionary with key xmax
        @test haskey(dfs, :xmax)

        # Info for non-stationary model but does not throw
        @test_logs (:info, "The graph is optimized for stationary models and the model provided is not.") histplot_data(fmns)

    end

    @testset "histplot(fm)" begin
        # Plot does not throw
        @test_logs histplot(fm)

    end

    @testset "diagnosticplots(fm)" begin
        # Plots do not throw
        @test_logs diagnosticplots(fm)

    end

    @testset "mrlplot_data(y)" begin

        σ = 1.0
        ξ = 0.1
        pd = GeneralizedPareto(σ, ξ)
        y = rand(pd, 1000)

        @test_throws AssertionError mrlplot_data(y, -3)

        df = mrlplot_data(y, 10)

        # Returns a dictionary with the correct keys
        @test hasproperty(df, :Threshold)
        @test hasproperty(df, :mrl)
        @test hasproperty(df, :lbound)
        @test hasproperty(df, :ubound)

        # Theoretical values
        m = (σ .+ ξ * df[:, :Threshold]) / (1 - ξ)

        @test df[1, :lbound] < m[1] < df[1, :ubound]
        @test df[5, :lbound] < m[5] < df[5, :ubound]

    end

    @testset "mrlplot(y)" begin
        # Plots do not throw
        @test_logs mrlplot([.5; .8; .9])

    end

end


text:
# An implementation of CGS for the solution of the square linear system Ax = b.
#
# This method is described in
#
# Y. Saad, Iterative methods for sparse linear systems.
# PWS Publishing Company, Boston, USA, 1996.
#
# P. Sonneveld, CGS, A Fast Lanczos-Type Solver for Nonsymmetric Linear systems.
# SIAM Journal on Scientific and Statistical Computing, 10(1), pp. 36--52, 1989.
#
# Alexis Montoison, <alexis.montoison@polymtl.ca>
# Montreal, October 2018.

export cgs

"""Solve the consistent linear system Ax = b using conjugate gradient squared algorithm.

From "Iterative Methods for Sparse Linear Systems (Y. Saad)" :

«The method is based on a polynomial variant of the conjugate gradients algorithm.
Although related to the so-called bi-conjugate gradients (BCG) algorithm,
it does not involve adjoint matrix-vector multiplications, and the expected convergence
rate is about twice that of the BCG algorithm.

The Conjugate Gradient Squared algorithm works quite well in many cases.
However, one difficulty is that, since the polynomials are squared, rounding errors
tend to be more damaging than in the standard BCG algorithm. In particular, very
high variations of the residual vectors often cause the residual norms computed
to become inaccurate.

TFQMR and BICGSTAB were developed to remedy this difficulty.»

This implementation allows a right preconditioner M.
"""
function cgs(A :: AbstractLinearOperator, b :: AbstractVector{T};
             M :: AbstractLinearOperator=opEye(size(A,1)),
             atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
             itmax :: Int=0, verbose :: Bool=false) where {T <: Number}

  m, n = size(A)
  m == n || error("System must be square")
  length(b) == m || error("Inconsistent problem size")
  verbose && @printf("CGS: system of size %d\n", n)

  # Initial solution x₀ and residual r₀.
  x = zeros(T, n) # x₀
  r = copy(b)     # r₀
  # Compute ρ₀ = < r₀,r₀ > and residual norm ‖r₀‖₂.
  ρ = @kdot(n, r, r)
  rNorm = sqrt(ρ)
  rNorm == 0 && return x, SimpleStats(true, false, [rNorm], [], "x = 0 is a zero-residual solution")

  iter = 0
  itmax == 0 && (itmax = 2*n)

  rNorms = [rNorm;]
  ε = atol + rtol * rNorm
  verbose && @printf("%5d  %7.1e\n", iter, rNorm)

  # Set up workspace.
  u = copy(r)  # u₀
  p = copy(r)  # p₀
  q = zeros(n) # q₋₁

  # Stopping criterion.
  solved = rNorm ≤ ε
  tired = iter ≥ itmax
  status = "unknown"

  while !(solved || tired)

    y = M * p                 # yₘ = M⁻¹pₘ
    v = A * y                 # vₘ = Ayₘ
    σ = @kdot(n, v, b)        # σₘ = < AM⁻¹pₘ,r₀ >
    α = ρ / σ                 # αₘ = ρₘ / σₘ
    @. q = u - α * v          # qₘ = uₘ - αₘ * AM⁻¹pₘ
    @kaxpy!(n, 1.0, q, u)     # uₘ₊½ = uₘ + qₘ
    z = M * u                 # zₘ = M⁻¹uₘ₊½
    @kaxpy!(n, α, z, x)       # xₘ₊₁ = xₘ + αₘ * M⁻¹(uₘ + qₘ)
    w = A * z                 # wₘ = AM⁻¹(uₘ + qₘ)
    @kaxpy!(n, -α, w, r)      # rₘ₊₁ = rₘ - αₘ * AM⁻¹(uₘ + qₘ)
    ρ_next = @kdot(n, r, b)   # ρₘ₊₁ = < rₘ₊₁,r₀ >
    β = ρ_next / ρ            # βₘ = ρₘ₊₁ / ρₘ
    @. u = r + β * q          # uₘ₊₁ = rₘ₊₁ + βₘ * qₘ
    @kaxpby!(n, 1.0, q, β, p) # pₘ₊₁ = uₘ₊₁ + βₘ * (qₘ + βₘ * pₘ)
    @kaxpby!(n, 1.0, u, β, p)

    # Update ρ.
    ρ = ρ_next # ρₘ ← ρₘ₊₁

    # Update iteration index.
    iter = iter + 1

    # Compute residual norm ‖rₘ‖₂.
    rNorm = @knrm2(n, r)
    push!(rNorms, rNorm)

    # Update stopping criterion.
    solved = rNorm ≤ ε
    tired = iter ≥ itmax
    verbose && @printf("%5d  %7.1e\n", iter, rNorm)
  end
  verbose && @printf("\n")

  status = tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol"
  stats = SimpleStats(solved, false, rNorms, T[], status)
  return (x, stats)
end


text:
using SyslogLogging
using Test
using Logging
using Sockets

function do_logging()
    @debug("hello syslog")
    @info("hello syslog")
    @warn("hello syslog")
    @error("hello syslog")
    @debug("hello syslog", param1="param1", param2=100, dict=Dict("abc"=>12, "def"=>20))
    @info("hello syslog", param1="param1", param2=100, dict=Dict("abc"=>12, "def"=>20))
    @warn("hello syslog", param1="param1", param2=100, dict=Dict("abc"=>12, "def"=>20))
    @error("hello syslog", param1="param1", param2=100, dict=Dict("abc"=>12, "def"=>20))
    nothing
end

function test_single_logger()
    logger = SyslogLogger("sysloglogger")
    with_logger(do_logging, logger)
end

function test_multiple_loggers()
    lck = ReentrantLock()
    logger1 = SyslogLogger("sysloglogger1"; lck=lck)
    logger2 = SyslogLogger("sysloglogger2"; lck=lck)
    @sync begin
        @async with_logger(do_logging, logger1)
        @async with_logger(do_logging, logger2)
    end
end

function test_remote_syslog_constructor()
    for host in ("127.0.0.1", ip"127.0.0.1", "localhost")
        @test_throws Base.IOError SyslogLogger("syslogloggerremote", host=host, port=10000, tcp=true) # ECONNREFUSED
        SyslogLogger("syslogloggerremote", host=host, port=10000, tcp=false)
    end
end

test_remote_syslog_constructor()
test_single_logger()
test_multiple_loggers()


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

# SPDX-License-Identifier: MIT
# Run this script with each new Julia release to update "../julia.spdx.json"

using UUIDs
using Dates
using JSON
using TimeZones
using DataStructures

spdxDocument= "../julia.spdx.json"
spdxData= JSON.parsefile(spdxDocument; dicttype=OrderedDict{String, Any})

# At the moment we can only update a few items automatically with each release.
# These are the crucial elements to make a new version of the SPDX file.
# Any other changes (ex. Adding or removing of external dependencies, updating copyright text, etc.) must be performed manually
spdxData["documentNamespace"]= "https://julialang.org/spdxdocs/julia-spdx-" * string(uuid4())
spdxData["creationInfo"]["created"]=  Dates.format(now(tz"UTC"), "yyyy-mm-ddTHH:MM:SS") * "Z"

for pkg in spdxData["packages"]
    if pkg["SPDXID"] == "SPDXRef-JuliaMain"
        pkg["versionInfo"]= readline("../VERSION")
        pkg["downloadLocation"]= "git+https://github.com/JuliaLang/julia.git@v" * pkg["versionInfo"]
        break
    end
end

open(spdxDocument, "w") do f
    JSON.print(f, spdxData, 4)
end

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
function bessely0(x::BigFloat)
    if x <= zero(x)
        if iszero(x)
            return -Inf
        else
            return throw(DomainError(x, "NaN result for non-NaN input."))
        end
    elseif isinf(x)
        return zero(x)
    end

    xx = abs(x)

    if xx <= 2.0
        z = xx * xx
        p = evalpoly(z, Y0_2N) / evalpoly(z, Y0_2D)
        p = TWOOPI(BigFloat) * log(x) * besselj0(x) + p
        return p
    end

    xinv = inv(xx)
    z = xinv * xinv
    if xinv <= 0.25
        if xinv <= 0.125
            if xinv <= 0.0625
                p = evalpoly(z, P16_IN) / evalpoly(z, P16_ID)
                q = evalpoly(z, Q16_IN) / evalpoly(z, Q16_ID)
            else
                p = evalpoly(z, P8_16N) / evalpoly(z, P8_16D)
                q = evalpoly(z, Q8_16N) / evalpoly(z, Q8_16D)
            end
        elseif xinv <= 0.1875
            p = evalpoly(z, P5_8N) / evalpoly(z, P5_8D)
	        q = evalpoly(z, Q5_8N) / evalpoly(z, Q5_8D)
        else
            p = evalpoly(z, P4_5N) / evalpoly(z, P4_5D)
	        q = evalpoly(z, Q4_5N) / evalpoly(z, Q4_5D)
        end
    else
        if xinv <= 0.375
            if xinv <= 0.3125
                p = evalpoly(z, P3r2_4N) / evalpoly(z, P3r2_4D)
	            q = evalpoly(z, Q3r2_4N) / evalpoly(z, Q3r2_4D)
            else
                p = evalpoly(z, P2r7_3r2N) / evalpoly(z, P2r7_3r2D)
	            q = evalpoly(z, Q2r7_3r2N) / evalpoly(z, Q2r7_3r2D)
            end
        elseif xinv <= 0.4375
            p = evalpoly(z, P2r3_2r7N) / evalpoly(z, P2r3_2r7D)
	        q = evalpoly(z, Q2r3_2r7N) / evalpoly(z, Q2r3_2r7D)
        else
            p = evalpoly(z, P2_2r3N) / evalpoly(z, P2_2r3D)
	        q = evalpoly(z, Q2_2r3N) / evalpoly(z, Q2_2r3D)
        end
    end

    p = 1.0 + z * p
    q = z * xinv * q
    q = q - 0.125 * xinv
    c = cos(x)
    s = sin(x)
    ss = s - c
    cc = s + c
    z = -cos(x + x)
    if s * c < zero(x)
        cc = z / ss
    else
        ss = z / cc
    end
    z = ONEOSQPI(BigFloat) * (p * ss + q * cc) / sqrt(x)
    return z
end


text:
using IICGP
using Test

cfg_filename = string(@__DIR__, "/dualcgp_test.yaml")
game = "assault"

main_cfg, enco_cfg, cont_cfg, reducer, bootstrap = IICGP.dualcgp_config(
    cfg_filename, game
)
enco = IPCGPInd(enco_cfg)
cont = CGPInd(cont_cfg)
@test true


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenLibm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenLibm")
JLLWrappers.@generate_main_file("OpenLibm", UUID("05823500-19ac-5b8b-9628-191a04bc5112"))
end  # module OpenLibm_jll

max_stars_repo_name: giordano/JLLWrappers

max_issues_repo_name: giordano/JLLWrappers

max_forks_repo_name: giordano/JLLWrappers


text:
"This scripts computes the Wigner little-d functions for the rotations"

export dd, du

const β = π/2

"Lower triangle of modes. The counters are bounded as n ≤ 2j, and |m| ≤ n in steps of 2."

"coeff₁ to coeff₄ are the coefficients for the three-term recursions of the Wigner little-d"

coeff₁(n, k) = (2 * (-n/2 - (k - n/2 - 1) * cos(β)) * csc(β)) / √(k * (n - k + 1))

coeff₂(n, k) = √((k - 1) * (n + 2 - k) / (k * (n + 1 - k)))

coeff₃(n, k, l) = (2 * (1 - k + n/2 + (l - n/2 - 1) * cos(β)) * csc(β)) / √(l * (n - l + 1))

coeff₄(n, k) = √((k - 1) * (n + 2 - k) / (k * (n + 1 - k)))

"`dd(n, j)` is a function who depends on the energy lever _n_ and gives an _n × n_ little-d coefficient matrix. n ∈ [0, 2j] for j ≥ 0. Warning: j is tested to be accurate below 28, that is, matrices of 56 points are safe to be considered as accurate."
function dd(n, j)
    vdd = zeros((2 * j + 1, 2 * j + 1))
    vdd[1, 1] = cos(β/2)^n
    vdd[2, 1] = √n * (cos(β) - 1) * csc(β) * vdd[1, 1]
    for k in 2:n
        vdd[k + 1, 1] = coeff₁(n, k) * vdd[k, 1] - coeff₂(n ,k) * vdd[k - 1, 1]
    end
    for k in 1:(n + 1)
        for l in 1:n
            if l == 1
                vdd[k, l + 1] = coeff₃(n, k, l) * vdd[k, l]
            else
                vdd[k, l + 1] = coeff₃(n, k, l) * vdd[k, l] - coeff₄(n, l) * vdd[k, l - 1]
            end
        end
    end
    return vdd
end

# For j = 28 (N = 56), there's no overflow. Need to check

coeff₅(n, k, j) = (2 * ((1/2) * (n - 4 * j) - (k - 1 + (1/2) * (n - 4 * j)) * cos(β)) * csc(β)) / √(k * (4 * j + 1 - k - n))

coeff₆(n, k, j) = √(((k - 1) * (4 * j + 2 - k - n)) / (k * (4 * j + 1 - k - n)))

coeff₇(n, k, l, j) = (2 * (1 - k + (1/2) * (4 * j - n) + (l - 1 + (1/2) * (n - 4 * j)) * cos(β)) * csc(β)) / √(l * (4 * j + 1 - l - n))

coeff₈(n, k, j) = √(((k - 1) * (4 * j + 2 - k - n)) / (k * (4 * j + 1 - k - n)))

"`du(n, j)` is a function who depends on the energy lever _n_ and gives an _n × n_ little-d coefficient matrix. n ∈ [2j, 4j] for j ≥ 0. This modes are counted in reverse, n = 4j is the lowest, and n = 2j is the highest. Warning: j is tested to be accurate below 28, that is, matrices of 56 points are safe to be considered as accurate."
function du(n, j)
    vdu = zeros((2 * j + 1, 2 * j + 1))
    vdu[1, 1] = cos(β/2)^(4 * j - n)
    vdu[2, 1] = √(4 * j - n) * (cos(β) - 1) * csc(β) * vdu[1, 1]
    for k in 2:(4 * j - n)
        vdu[k + 1, 1] = coeff₅(n, k, j) * vdu[k, 1] - coeff₆(n, k, j) * vdu[k - 1, 1]
    end
    for k in 1:(4 * j + 1 - n)
        for l in 1:(4 * j - n)
            if l == 1
                vdu[k, l + 1] = coeff₇(n, k, l, j) * vdu[k, l]
            else
                vdu[k, l + 1] = coeff₇(n, k, l, j) * vdu[k, l] - coeff₈(n, l, j) * vdu[k, l - 1]
            end
        end
    end
    return vdu
end


text:
module WordErrorRate

using Printf: @sprintf

export WER, wer, align, pralign

struct WER
    ref::Vector
    hyp::Vector
    nsub::Int
    nins::Int
    ndel::Int
    refeval::Vector{UInt8}
    hypeval::Vector{UInt8}
end

const CORRECT, SUBSTITUTION, INSERTION, DELETION = 0:3

"Compute the Word Error Rate similar to the NIST sclite package"
function WER(ref::Vector, hyp::Vector; subcost=4, inscost=3, delcost=3)
    nrow = length(ref) + 1
    ncol = length(hyp) + 1
    ## initialization
    cost = zeros(UInt16, nrow, ncol)
    decision = zeros(UInt8, nrow, ncol)
    for i in 2:nrow
        cost[i, 1] = (i-1) * delcost
        decision[i, 1] = DELETION
    end
    for j in 2:ncol
        cost[1, j] = (j-1) * inscost
        decision[1, j] = INSERTION
    end
    ## computation
    for i in 2:nrow
        for j in 2:ncol
            if ref[i-1] == hyp[j-1]
                cost[i, j] = cost[i-1, j-1]
                decision[i, j] = CORRECT
            else
                sid = [cost[i-1, j-1] + subcost, cost[i, j-1] + inscost, cost[i-1, j] + delcost]
                mini = argmin(sid)
                cost[i, j] = sid[mini]
                decision[i, j] = mini
            end
        end
    end
    ## backtrack
    refeval = UInt8[]
    hypeval = UInt8[]
    i, j = nrow, ncol
    nsub = nins = ndel = 0
    while (i > 1) || (j > 1)
        dec = decision[i, j]
        if dec < 2
            nsub += dec
            i -= 1
            j -= 1
            pushfirst!(refeval, dec)
            pushfirst!(hypeval, dec)
        elseif dec == INSERTION ## ins
            j -= 1
            pushfirst!(hypeval, dec)
            nins += 1
        else
            i -= 1
            pushfirst!(refeval, dec)
            ndel += 1
        end
    end

    return WER(ref, hyp, nsub, nins, ndel, refeval, hypeval)
end

wer(w::WER) = (w.nsub + w.nins + w.ndel) / length(w.ref)
wer(ref::Array, hyp::Array; kwargs...) = wer(WER(ref, hyp; kwargs...))

correct(w::WER) = length(w.ref) - w.nsub - w.ndel

function align(w::WER)
    refali = Int[]
    hypali = Int[]
    evalali = Int8[]
    ri = hi = 1
    while ri ≤ length(w.refeval) || hi ≤ length(w.hypeval)
        if ri ≤ length(w.refeval) && w.refeval[ri] == DELETION
            push!(hypali, 0)
            push!(refali, ri)
            push!(evalali, DELETION)
            ri += 1
        elseif hi ≤ length(w.hypeval) && w.hypeval[hi] == INSERTION
            push!(refali, 0)
            push!(hypali, hi)
            push!(evalali, INSERTION)
            hi += 1
        else
            push!(refali, ri)
            push!(hypali, hi)
            push!(evalali, w.refeval[ri])
            ri += 1
            hi += 1
        end
    end
    return refali, hypali, evalali
end

"""print alignment like sclite does"""
function pralign(io::IO, w::WER; ins="*", del="*")
    println(io, "Scores: (#C #S #D #I) ", correct(w), " ", w.nsub, " ", w.ndel, " ", w.nins)
    refali, hypali, evalali = align(w)
    prref = [i > 0 ? w.ref[i] : del^length(string(w.hyp[j])) for (i, j) in zip(refali, hypali)]
    prhyp = [j > 0 ? w.hyp[j] : ins^length(string(w.ref[i])) for (i, j) in zip(refali, hypali)]
    preval = [" SID"[i+1] for i in evalali]
    lengths = [maximum(length(string(s)) for s in x) for x in zip(prref, prhyp, preval)]
    for (name, a) in zip(("REF: ", "HYP: ", "Eval:"), (prref, prhyp, preval))
        print(io, name, " ")
        for (w, l) in zip(a, lengths)
            s = string(w)
            print(io, s, " "^(l+1-length(s)))
        end
        println(io)
    end
end

pralign(w::WER; kwargs...) = pralign(stdout, w; kwargs...)
function pralign(::Type{String}, w::WER; kwargs...)
    io = IOBuffer()
    pralign(io, w; kwargs...)
    return String(take!(io))
end

Base.show(io::IO, w::WER) = print(io, @sprintf("Word Error Rate: %4.2f%%, Ns, Ni, Nd = (%d, %d, %d)", 100wer(w), w.nsub, w.nins, w.ndel))

end


text:
using DifferentialEquations, Plots

A = [
  1. 0 0 -5
  4 -2 4 -3
  -4 0 0 1
  5 -2 2 3
]
u0 = rand(4, 2)
tspan = (0.0, 1.0)
f(t,u) = A*u

prob = ODEProblem(f, u0, tspan)

sol = solve(prob)
plot(sol)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "bsdiff_endsley"
version = v"4.3.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/mendsley/bsdiff/archive/64ad986cb7bfa8b9145a2d48cd95986660b35d53.tar.gz", "1181466689aa224f4a2dd2376820588c67d20f4f0d50055339fcb171fb877a29"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/bsdiff-64ad986cb7bfa8b9145a2d48cd95986660b35d53
./autogen.sh
export CPPFLAGS="-I${prefix}/include"
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# Disable Windows for now, as there are many BSD-isms in the source code
# that we don't want to bother to patch out.  Things like err.h and whatnot.
platforms = filter(!Sys.iswindows, supported_platforms())

# The products that we will ensure are always built
products = [
    ExecutableProduct("bspatch", :bspatch),
    ExecutableProduct("bsdiff", :bsdiff)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # Future versions of bzip2 should allow a more relaxed compat because the
    # soname of the macOS library shouldn't change at every patch release.
    Dependency("Bzip2_jll", v"1.0.6"; compat="=1.0.6"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
### A Pluto.jl notebook ###
# v0.17.2

using Markdown
using InteractiveUtils

# ╔═╡ ab30fd04-508e-11ec-3026-a5b39dae84a8
md"""
## Quickstart
Test
"""

# ╔═╡ Cell order:
# ╟─ab30fd04-508e-11ec-3026-a5b39dae84a8


text:
function bauman_form(X::IntervalBox{N, T}, f::Function, g::SVector{N, Interval{T}}) where{N, T}
    cb = Float64[]

    for i in 1:length(g)
        if g[i] >= 0
            append!(cb, X[i].lo)
        elseif g[i] <= 0
            append!(cb, X[i].hi)
        else
            append!(cb,(g[i].hi * X[i].lo - g[i].lo * X[i].hi)/(g[i].hi - g[i].lo))
        end
    end

    (sum1, sum2)  = (0, 0)

    for i in 1:length(g)
        sum1 = sum1 + (g[i].hi)*(X[i].lo - cb[i])
        sum2 = sum2 + (g[i].lo)*(X[i].hi - cb[i])
    end

    func_val = f(cb)
    static_cb = SVector{length(cb)}(cb)

    lbb = func_val + min(sum1, sum2)

    return (lbb, func_val, static_cb)

end


text:
function inverseSide(side::HorizonSide)
    # Return inverse side
    return HorizonSide(mod(Int(side)+2, 4))
end

function main(r::Robot)
    # Пусть робот окажется в левом верхнем углу
    stepsToWest = moveSide(r, West)
    stepsToNord = moveSide(r, Nord)

    # Двигаемся по часовой стрелке
    sides = (Ost, Sud, West, Nord, Ost)
    for side in sides
        while isborder(r, side) != true
            if ismarker(r)
                followByCount(r, Ost, stepsToWest)
                followByCount(r, Sud, stepsToNord)
                return true
            end
            moveSide(r, side, true)
        end
    end
end

function moveSide(r::Robot, side::HorizonSide, needMark=false)
    # Move untill you find pr
    steps = 0
    while isborder(r, side) != true
        moveStep(r, side, needMark)
        steps += 1
    end
    return steps
end

function moveStep(r::Robot, side::HorizonSide, needMark::Bool)
    # Take step in one direction
    if needMark
        putmarker!(r)
    end
    move!(r, side)
end

function followByCount(r::Robot, side::HorizonSide, count::Int)
    for i = 0:count-1
        moveStep(r, side, false)
    end
end

max_stars_repo_name: TheDarkMoon/mirea-progs

max_issues_repo_name: TheDarkMoon/mirea-progs

max_forks_repo_name: TheDarkMoon/mirea-progs


text:
@define_scenario MaintenanceSchedulingScenario = begin
    ρ::PriceCurve{24,Float64}
    Q̃::Inflows{typeof(Skellefteälven),Float64}

    @expectation begin
        ρ = mean([s.ρ for s in scenarios])
        Q̃ = mean([s.Q̃ for s in scenarios])
        return MaintenanceSchedulingScenario(ρ, Q̃, probability = 1.0)
    end
end

function StochasticPrograms.expected(ξ₁::MaintenanceSchedulingScenario, ξ₂::MaintenanceSchedulingScenario)
    ρ = mean([ξ₁.ρ, ξ₂.ρ])
    Q̃ = mean([ξ₁.Q̃, ξ₂.Q̃])
    expected = MaintenanceSchedulingScenario(ρ, Q̃)
    return ExpectedScenario(expected)
end

function penalty(scenario::MaintenanceSchedulingScenario, t)
    if 8 <= t <= 12 || 17 <= t <= 24
        return 1.15*scenario.ρ[t]
    else
        return 1.1*scenario.ρ[t]
    end
end
function reward(scenario::MaintenanceSchedulingScenario, t)
    if 8 <= t <= 12 || 17 <= t <= 24
        return 0.85*scenario.ρ[t]
    else
        return 0.9*scenario.ρ[t]
    end
end

@sampler RecurrentMaintenanceSchedulingSampler = begin
    date::Date
    plants::PlantCollection
    price_forecaster::Forecaster
    flow_forecaster::Forecaster

    @sample MaintenanceSchedulingScenario begin
        price_curve = forecast(sampler.price_forecaster, month(sampler.date))
        while !all(price_curve .>= 0)
            price_curve = forecast(sampler.price_forecaster, month(sampler.date))
        end
        flows = forecast(sampler.flow_forecaster, week(sampler.date))
        while !all(flows[:, 2] .>= 0)
            flows = forecast(sampler.flow_forecaster, week(sampler.date))
        end
        return MaintenanceSchedulingScenario(PriceCurve(price_curve), Inflows(sampler.plants, flows[:, 2]))
    end
end


text:
lines = readlines("inputs/03.txt")
all_digits = [0,0,0,0,0,0,0,0,0,0,0,0]
for x in lines
  digits = [parse(Int64, l) for l in x]
  for d in 1:12
    all_digits[d] += digits[d]
  end
end
norm = [x/1000 for x in all_digits]

bits  = [x >= .5 ? 1 : 0 for x in norm]
bits2 = [x < .5 ? 1 : 0 for x in norm]

for b in bits
  print(b)
end

for b in bits2
  print(b)
end


text:
using Documenter
using IndexedDims
using Test

@testset "IndexedDims.jl" begin
    doctest(IndexedDims)

    a = IndexedDimsArray([1 2 3; 4 5 6; 7 8 9], 2:4, 2:4)

    # constructor
    @test ==(
        IndexedDimsArray([1 2; 3 4; 5 6]),
        IndexedDimsArray([1 2; 3 4; 5 6], Base.OneTo(2), Base.OneTo(3))
    )

    # standard fallback
    @test a[1] == 1
    @test a[3, 2] == 8

    # fancy
    @test a[:, :] == a
    @test typeof(a[:, :]) == typeof(a)
    @test a[3:4, 2:3] == [4 5; 7 8]
    @test a[3:4, :] == [4 5 6; 7 8 9]
    @test a[3:4, bypass(1:2)] == [4 5; 7 8]
end


text:
# TableLogger

mutable struct TableCol
    id::Symbol
    name::String
    type::DataType
    width::Int
    alignment::Symbol # :left, :center, :right
    b_format::Bool
end

mutable struct TableEntry{T}
    col_id::Symbol
    value::T
end

mutable struct TableSetup
    cols::Vector{TableCol}
    col_idx::Dict{Symbol,Int}
    new_row_criteria::Bool
    diff_criteria::Dict{Symbol,Any}
    last_row::Vector{TableEntry}
end

# SolverOptions
mutable struct ActivityOptions
    decay::Float64
    max_probes::Int
    max_confidence_deviation::Float64
end

mutable struct SolverOptions
    logging::Vector{Symbol}
    table::TableSetup
    time_limit::Float64 # time limit in backtracking in seconds
    seed::Int
    traverse_strategy::Symbol
    branch_strategy::Symbol
    branch_split::Symbol # defines splitting in the middle, or takes smallest, biggest value
    backtrack::Bool
    max_bt_steps::Int
    backtrack_sorting::Bool
    keep_logs::Bool
    rtol::Float64
    atol::Float64
    solution_type::Type
    all_solutions::Bool
    all_optimal_solutions::Bool
    lp_optimizer::Any
    no_prune::Bool
    activity::ActivityOptions
    simplify::Bool
end

# General
"""
Save the changes for a variable with two vectors
- `indices` saves a list of indices that can be extracted by the step_nr
    - i.e [1,1,3,5] would mean that step_nr 1 did no changes to the variable and steps 2 and 3 have two changes each
- `changes`:
    Explanation of the tuple
    - [1] :fix, :rm, :rm_below, :rm_above
    - [2] To which value got it fixed, which value was removed, which value was the upper/lower bound
    - [3] Only if fixed it saves the last ptr to revert the changes otherwise 0
    - [4] How many values got removed (0 for fix)
"""
mutable struct VariableChanges
    indices :: Vector{Int}
    changes :: Vector{Tuple{Symbol, Int, Int, Int}}
    VariableChanges() = new(Int[1,1], Tuple{Symbol, Int, Int, Int}[])
end

mutable struct Variable
    idx::Int
    lower_bound::Int # initial lower and
    upper_bound::Int # upper bound of the variable see min, max otherwise
    first_ptr::Int
    last_ptr::Int
    values::Vector{Int}
    indices::Vector{Int}
    init_vals::Vector{Int} # saves all initial values
    init_val_to_index::Vector{Int} # saves the index in which val appears in init_vals
    offset::Int
    min::Int # the minimum value during the solving process
    max::Int # for initial see lower/upper_bound
    changes::VariableChanges
    has_upper_bound::Bool # must be true to work
    has_lower_bound::Bool # must be true to work
    is_fixed::Bool
    is_integer::Bool # must be true to work
    # branching strategies
    activity::Float64 #  + 1 if variable was used in node, * activity.decay if it wasn't
    in_all_different::Vector{Bool} # saves in which all different constraint this variable is used
                                   # can be used to check if two or more variables are in the same alldifferent constraint
end

mutable struct NumberConstraintTypes
    equality::Int
    inequality::Int
    notequal::Int
    alldifferent::Int
    table::Int
    indicator::Int
    reified::Int
end

mutable struct CSInfo
    pre_backtrack_calls::Int
    backtracked::Bool
    backtrack_fixes::Int
    in_backtrack_calls::Int
    backtrack_reverses::Int
    n_constraint_types::NumberConstraintTypes
end

#====================================================================================
=================== SETS FOR VARIABLES AND CONSTRAINTS ==============================
====================================================================================#

struct Integers <: MOI.AbstractScalarSet
    values::Vector{Int}
end
Integers(vals::Union{UnitRange{Int},StepRange{Int,Int}}) = Integers(collect(vals))
Base.copy(I::Integers) = Integers(I.values)

struct AllDifferent <: JuMP.AbstractVectorSet end
JuMP.moi_set(::AllDifferent, dim) = CPE.AllDifferent(dim)

struct TableSetInternal <: MOI.AbstractVectorSet
    dimension::Int
    table::Array{Int,2}
end
Base.copy(T::TableSetInternal) = TableSetInternal(T.dimension, T.table)

struct TableSet <: JuMP.AbstractVectorSet
    table::Array{Int,2}
end
function JuMP.moi_set(ts::TableSet, dim)
    if size(ts.table, 2) != dim
        throw(ArgumentError("The table provided has $(size(ts.table,2)) columns but the variable vector has $dim elements"))
    end
    TableSetInternal(dim, ts.table)
end

struct GeqSetInternal <: MOI.AbstractVectorSet
    dimension::Int
end
Base.copy(G::GeqSetInternal) = GeqSetInternal(G.dimension)

struct GeqSet <: JuMP.AbstractVectorSet end
JuMP.moi_set(::GeqSet, dim) = GeqSetInternal(dim)

struct AllEqual <: JuMP.AbstractVectorSet end
JuMP.moi_set(::AllEqual, dim) = CPE.AllEqual(dim)

#====================================================================================
====================== TYPES FOR TRAVERSING ========================================
====================================================================================#

abstract type Priority end

struct PriorityDFS{T<:Real} <: Priority
    depth::Int
    bound::T
    neg_idx::Int # the negative backtrack index (negative because of maximizing)
end

function Base.isless(p1::PriorityDFS, p2::PriorityDFS)
    if p1.depth < p2.depth
        return true
    elseif p1.depth == p2.depth
        if p1.bound < p2.bound
            return true
        elseif p1.bound == p2.bound
            return p1.neg_idx < p2.neg_idx
        else
            return false
        end
    end
    return false
end

struct PriorityBFS{T<:Real} <: Priority
    bound::T
    depth::Int
    neg_idx::Int # the negative backtrack index (negative because of maximizing)
end

function Base.isless(p1::PriorityBFS, p2::PriorityBFS)
    if p1.bound < p2.bound
        return true
    elseif p1.bound == p2.bound
        if p1.depth < p2.depth
            return true
        elseif p1.depth == p2.depth
            return p1.neg_idx < p2.neg_idx
        else
            return false
        end
    end
    return false
end

#====================================================================================
====================== TYPES FOR CONSTRAINTS ========================================
====================================================================================#

abstract type Constraint end

"""
    BoundRhsVariable
idx - variable index in the lp Model
lb  - lower bound of that variable
ub  - upper bound of that variable
"""
mutable struct BoundRhsVariable
    idx::Int
    lb::Int
    ub::Int
end

mutable struct MatchingInit
    l_in_len::Int
    matching_l::Vector{Int}
    matching_r::Vector{Int}
    index_l::Vector{Int}
    process_nodes::Vector{Int}
    depths::Vector{Int}
    parents::Vector{Int}
    used_l::Vector{Bool}
    used_r::Vector{Bool}
end

mutable struct SCCInit
    index_ei::Vector{Int}
    ids::Vector{Int}
    low::Vector{Int}
    on_stack::Vector{Bool}
    group_id::Vector{Int}
end

"""
    RSparseBitSet

See https://arxiv.org/pdf/1604.06641.pdf
words[x] will save 64 possibilities in a TableConstraint a `1` at position y will mean that the
words in row (x-1)*64+y of the table are possible a 0 means that they aren't
Similar to `Variable` a whole block of 64 rows can be removed by changing the `indices` and `last_ptr`.
The `mask` saves the current mask to change words
"""
mutable struct RSparseBitSet
    words::Vector{UInt64}
    indices::Vector{Int}
    last_ptr::Int
    mask::Vector{UInt64}
    RSparseBitSet() = new([], [], 1, [])
end

mutable struct TableSupport
    # defines the range for each variable
    # i.e [1,3,7,10] means that the first variable has 2 values, the second 4
    var_start::Vector{Int}
    values::Array{UInt64,2}
    TableSupport() = new([], Array{UInt64,2}(undef, (0, 0)))
end

mutable struct TableResidues
    # defines the range for each variable
    # i.e [1,3,7,10] means that the first variable has 2 values, the second 4
    var_start::Vector{Int}
    values::Vector{Int}
    TableResidues() = new([], [])
end

mutable struct TableBacktrackInfo
    words::Vector{UInt64}
    last_ptr::Int
    indices::Vector{Int}
end

struct IndicatorSet{A,F,S<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet}} <: MOI.AbstractVectorSet
    set::S
    dimension::Int
end
IndicatorSet{A,F}(set::S) where {A,F,S} = IndicatorSet{A,F,S}(set, 1+MOI.dimension(set))
Base.copy(I::IndicatorSet{A,F,S}) where {A,F,S} = IndicatorSet{A,F,S}(I.set, I.dimension)

struct ReifiedSet{A,F,S<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet}} <:
       MOI.AbstractVectorSet
    set::S
    dimension::Int
end
Base.copy(R::ReifiedSet{A,F,S}) where {A,F,S} = ReifiedSet{A,F,S}(R.set, R.dimension)

abstract type AbstractBoolSet{
    F1<:Union{SAF,VAF,MOI.VectorOfVariables},
    F2<:Union{SAF,VAF,MOI.VectorOfVariables},
    F1dim<:Val,
    F2dim<:Val,
    S1<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet},
    S2<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet},
} <: MOI.AbstractVectorSet end

struct BoolSetInternals{S1,S2}
    lhs_set::S1
    rhs_set::S2
    lhs_dimension::Int
    rhs_dimension::Int
    dimension::Int
end

const BOOL_SET_TO_CONSTRAINT = (
    :AndSet => (constraint = :AndConstraint, op = :(&&)),
    :OrSet  => (constraint = :OrConstraint, op = :(||)),
    :XorSet => (constraint = :XorConstraint, op = :(⊻), needs_call = true, specific_constraint = true),
    :XNorSet => (constraint = :XNorConstraint, res_op = :(!), op = :(⊻), needs_call = true, specific_constraint = true)
)

for (set, bool_data) in BOOL_SET_TO_CONSTRAINT
    @eval begin
        struct $set{
            F1<:Union{SAF,VAF,MOI.VectorOfVariables},
            F2<:Union{SAF,VAF,MOI.VectorOfVariables},
            F1dim<:Val,
            F2dim<:Val,
            S1<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet},
            S2<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet}
        } <: AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}
            bsi::BoolSetInternals{S1,S2}
        end
    end
end

function (::Type{BS})(lhs_set::S1, rhs_set::S2) where {S1,S2,F1,F2,BS<:CS.AbstractBoolSet{F1,F2,F1dim,F2dim,S1_,S2_} where {F1dim,F2dim,S1_,S2_}}
    lhs_dim = MOI.dimension(lhs_set)
    rhs_dim = MOI.dimension(rhs_set)
    F1dim = Val{lhs_dim}
    F2dim = Val{rhs_dim}
    internals = BoolSetInternals(lhs_set, rhs_set, lhs_dim, rhs_dim, lhs_dim + rhs_dim)
    return typeof_without_params(BS){F1,F2,F1dim,F2dim,S1,S2}(internals)
end

function Base.copy(A::AbstractBoolSet{F1,F2,F1dim,F2dim,S1,S2}) where {F1,F2,F1dim,F2dim,S1,S2}
    typeof_without_params(A){F1,F2,F1dim,F2dim,S1,S2}(A.bsi)
end

struct ComplementSet{F<:MOI.AbstractFunction, S<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet}} <: MOI.AbstractVectorSet
    set::S
    dimension::Int
end
ComplementSet{F}(set::S) where {F,S} = ComplementSet{F,S}(set, MOI.dimension(set))
Base.copy(A::ComplementSet{F,S}) where {F,S} = ComplementSet{F,S}(A.set, A.dimension)

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

mutable struct ImplementedConstraintFunctions
    init::Bool
    activate::Bool
    update_init::Bool
    finished_pruning::Bool
    restore_pruning::Bool
    single_reverse_pruning::Bool
    reverse_pruning::Bool
    update_best_bound::Bool
end

mutable struct ConstraintInternals{
    FCT<:Union{MOI.AbstractScalarFunction,MOI.AbstractVectorFunction},
    SET<:Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet},
}
    idx::Int
    fct::FCT
    set::SET
    indices::Vector{Int}
    pvals::Vector{Int}
    is_initialized::Bool
    is_activated::Bool
    is_deactivated::Bool # can be deactivated if it's absorbed by other constraints
    bound_rhs::Vector{BoundRhsVariable}# should be set if `update_best_bound` is true
end

#====================================================================================
====================== CONSTRAINTS ==================================================
====================================================================================#

mutable struct BasicConstraint <: Constraint
    std::ConstraintInternals
end

mutable struct EqualConstraint <: Constraint
    std::ConstraintInternals
    first_ptrs::Vector{Int} # for faster apply_changes!
end

mutable struct AllDifferentConstraint <: Constraint
    std::ConstraintInternals
    pval_mapping::Vector{Int}
    vertex_mapping::Vector{Int}
    vertex_mapping_bw::Vector{Int}
    di_ei::Vector{Int}
    di_ej::Vector{Int}
    matching_init::MatchingInit
    scc_init::SCCInit
    # corresponds to `in_all_different`: Saves the constraint idxs where all variables are part of this alldifferent constraint
    sub_constraint_idxs::Vector{Int}
end

abstract type BoolConstraint{C1<:Constraint,C2<:Constraint} <: Constraint end

mutable struct BoolConstraintInternals{C1<:Constraint,C2<:Constraint}
    lhs_activated::Bool
    lhs_activated_in_backtrack_idx::Int
    rhs_activated::Bool
    rhs_activated_in_backtrack_idx::Int
    lhs::C1
    rhs::C2
end

struct XorConstraint{C1,C2} <: BoolConstraint{C1,C2}
    std::ConstraintInternals
    bool_std::BoolConstraintInternals{C1,C2}
    complement_lhs::Union{Nothing, Constraint}
    complement_rhs::Union{Nothing, Constraint}
end

struct XNorConstraint{C1,C2} <: BoolConstraint{C1,C2}
    std::ConstraintInternals
    bool_std::BoolConstraintInternals{C1,C2}
    complement_lhs::Union{Nothing, Constraint}
    complement_rhs::Union{Nothing, Constraint}
end

for (set, bool_data) in BOOL_SET_TO_CONSTRAINT
    get(bool_data, :specific_constraint, false) && continue
    @eval begin
        struct $(bool_data.constraint){C1,C2} <: BoolConstraint{C1,C2}
            std::ConstraintInternals
            bool_std::BoolConstraintInternals{C1,C2}
        end
    end
end

# support for a <= b constraint
mutable struct SingleVariableConstraint <: Constraint
    std::ConstraintInternals
    lhs::Int
    rhs::Int
end

mutable struct LinearConstraint{T<:Real} <: Constraint
    std::ConstraintInternals
    in_all_different::Bool
    is_strict::Bool # for differentiate between < and <=
    is_equal::Bool # for ==
    is_rhs_strong::Bool # saves whether the rhs is already a strong rhs
    rhs::T # combines value - constant
    # same as rhs but for `<` it saves the original value - constant
    # (rhs is then computed to be usable as <=)
    strict_rhs::T
    coeffs::Vector{T}
    mins::Vector{T}
    maxs::Vector{T}
    pre_mins::Vector{T}
    pre_maxs::Vector{T}
end

mutable struct TableConstraint <: Constraint
    std::ConstraintInternals
    current::RSparseBitSet
    supports::TableSupport
    last_sizes::Vector{Int}
    residues::TableResidues
    # holds current, last_ptr and indices from each node
    # maybe it's better to compute some of them to save some space...
    # This is the easy implementation first
    backtrack::Vector{TableBacktrackInfo}
    changed_vars::Vector{Int}
    unfixed_vars::Vector{Int}
    sum_min::Vector{Int}
    sum_max::Vector{Int}
end

mutable struct GeqSetConstraint <: Constraint
    std::ConstraintInternals
    vidx::Int
    greater_than::Vector{Int}
    # saves all constraints where all indices are part of the greater_than
    sub_constraint_idxs::Vector{Int}
end

abstract type ActivatorConstraint <: Constraint end

mutable struct ActivatorConstraintInternals
    activate_on::MOI.ActivationCondition
    activator_in_inner::Bool
    inner_activated::Bool
    inner_activated_in_backtrack_idx::Int
end

mutable struct IndicatorConstraint{C<:Constraint} <: ActivatorConstraint
    std::ConstraintInternals
    act_std::ActivatorConstraintInternals
    inner_constraint::C
end

mutable struct ReifiedConstraint{C<:Constraint, AC<:Union{Constraint,Nothing}} <: ActivatorConstraint
    std::ConstraintInternals
    act_std::ActivatorConstraintInternals
    inner_constraint::C
    complement_constraint::AC
    complement_inner_constraint::Bool
    complement_inner_constraint_in_backtrack_idx::Int
end

#====================================================================================
====================== OBJECTIVES ==================================================
====================================================================================#

abstract type ObjectiveFunction end

mutable struct SingleVariableObjective <: ObjectiveFunction
    fct::MOI.SingleVariable
    vidx::Int # index of the variable
    indices::Vector{Int}
end

# used to designate a feasibility sense
struct NoObjective <: ObjectiveFunction end

mutable struct LinearCombination{T<:Real}
    indices::Vector{Int}
    coeffs::Vector{T}
end

mutable struct LinearCombinationObjective{T<:Real} <: ObjectiveFunction
    fct::MOI.ScalarAffineFunction
    lc::LinearCombination{T}
    constant::T
    indices::Vector{Int} # must exist to update the objective only if one of these changed
end

mutable struct BacktrackObj{T<:Real}
    idx::Int
    step_nr::Int
    parent_idx::Int
    depth::Int
    status::Symbol
    is_feasible::Bool
    vidx::Int
    lb::Int # lb <= var[vidx] <= ub
    ub::Int
    best_bound::T
    primal_start::Vector{Float64}
    solution::Vector{Float64} # holds the solution values of the bound computation
end

function Base.convert(::Type{B}, obj::BacktrackObj{T2}) where {T1,T2,B<:BacktrackObj{T1}}
    return BacktrackObj{T1}(
        obj.idx,
        obj.step_nr,
        obj.parent_idx,
        obj.depth,
        obj.status,
        obj.is_feasible,
        obj.vidx,
        obj.lb,
        obj.ub,
        convert(T1, obj.best_bound),
        obj.primal_start,
        obj.solution,
    )
end

mutable struct TreeLogNode{T<:Real}
    id::Int
    status::Symbol
    feasible::Bool
    best_bound::T
    step_nr::Int
    vidx::Int
    lb::Int
    ub::Int
    var_states::Dict{Int,Vector{Int}}
    var_changes::Dict{Int,Vector{Tuple{Symbol,Int,Int,Int}}}
    activity::Dict{Int,Float64}
    children::Vector{TreeLogNode{T}}
end

mutable struct Solution{T<:Real}
    incumbent::T
    values::Vector{Int}
    backtrack_id::Int # save where the solution was found
    hash::UInt64
end
Solution(incumbent, values, backtrack_id) =
    Solution(incumbent, values, backtrack_id, hash(values))

mutable struct ActivityObj
    nprobes::Int
    is_free::Vector{Bool}
    ActivityObj() = new(0, [false]) # will be overwritten later
end

"""
    BranchVarObj

Determines the next branch variable and stores if still feasible and if solution was found
"""
mutable struct BranchVarObj
    is_feasible::Bool
    is_solution::Bool
    vidx::Int # only relevant if is_feasible && !is_solution
end

struct VarAndVal
    vidx::Int
    lb::Int
    ub::Int
end

mutable struct ConstraintSolverModel{T<:Real}
    lp_model::Union{Nothing,Model} # only used if lp_optimizer is set
    lp_x::Vector{VariableRef}
    init_search_space::Vector{Variable}
    search_space::Vector{Variable}
    init_fixes::Vector{Tuple{Int,Int}}
    subscription::Vector{Vector{Int}}
    constraints::Vector{Constraint}
    root_infeasible_vars::Vector{VarAndVal}
    bt_infeasible::Vector{Int}
    c_backtrack_idx::Int
    c_step_nr::Int
    backtrack_vec::Vector{BacktrackObj{T}}
    backtrack_pq::PriorityQueue
    sense::MOI.OptimizationSense
    objective::ObjectiveFunction
    var_in_obj::Vector{Bool} # saves whether a variable is part of the objective function
    traverse_strategy::Val
    branch_strategy::Val
    branch_split::Val
    in_probing_phase::Bool
    activity_vars::ActivityObj
    best_sol::T # Objective of the best solution
    best_bound::T # Overall best bound
    solutions::Vector{Solution}
    info::CSInfo
    input::Dict{Symbol,Any}
    logs::Vector{TreeLogNode{T}}
    options::SolverOptions
    start_time::Float64
    solve_time::Float64 # seconds spend in solve
end
parametric_type(::ConstraintSolverModel{T}) where {T} = T


text:
struct Pol{T<:Number}
    coeffs::Vector{T}
    function Pol{T}(coeffs) where {T <: Number}
        length(coeffs) == 1 || removeTailZeroes!(coeffs)
        new(coeffs)
    end
end

Pol(coeffs::Vector{T}) where {T<:Number} = Pol{T}(coeffs)
Pol() = Pol([0])
Pol(n::Number) = Pol([n])

rev(p::Pol) = Pol(p.coeffs[end:-1:1])

# Get properties
coeffs(p::Pol) = p.coeffs
coeff(p::Pol, d::Integer) = p[d+1]

lc(p::Pol) = p.coeffs[end]

function lm(p::Pol{T}) where {T<:Number}
    coeffs = zeros(T, length(p.coeffs))
    coeffs[end] = p.coeffs[end]
    return Pol(coeffs)
end

deg(p::Pol) = length(p.coeffs)-1

function variable(x::String, t=Int)
    global _varname = x
    return one(t)*Pol([0, 1])
end

Base.length(p::Pol) = length(p.coeffs)

function subs(p::Pol{T}, x::S) where {T,S<:Number}
    res = lc(p)
    @inbounds for i in length(p)-1:-1:1
        res = x*res+p[i]
    end
    return res
end

# PROMOTION AND CONVERSION
Base.promote_rule(::Type{Pol{T}}, ::Type{S}) where {T, S <: Number} = Pol{promote_type(T, S)}
Base.promote_rule(::Type{Pol{T}}, ::Type{Pol{S}}) where {T, S <: Number} = Pol{promote_type(T, S)}
Base.convert(::Type{Pol{T}}, x::S) where {T,S<:Number} = Pol{promote_type(T, S)}([x])
Base.convert(::Type{Pol{T}}, p::Pol{S}) where {T, S <: Number} = Pol{promote_type(T, S)}(p.coeffs)

Base.zero(p::Pol{T}) where {T<:Number} = Pol(zero(T))
Base.zero(::Type{Pol{T}}) where {T<:Number} = Pol(zero(T))
Base.zero(::Type{Pol}) = Pol(zero(Int))
Base.one(::Type{Pol{T}}) where {T <: Number} = Pol(one(T))

# indexing
function Base.getindex(p::Pol{T}, i::Integer) where {T<:Number}
    i < 1 && throw(BoundsError(p, i))
    i > length(p.coeffs) && return zero(T)
    return p.coeffs[i]
end

# DISPLAYING
function Base.show(io::IO, p::Pol)
    print(io, string(p))
end

function Base.string(p::Pol)
    if iszero(p)
        return "0"
    else
        str = ""
        @inbounds for (i, c) in Iterators.reverse(enumerate(p.coeffs))
            iszero(c) && continue
            str *= c>0 ? "+" : "-"
            str *= (abs(c) != 1 || i==1) ? string(abs(c)) : ""
            str *= i==1 ? "" : _varname;
            str *= i<=2 ? "" : toSuperscript(i-1)
        end
        str = str[1]=='+' ? str[2:end] : str
        return str
    end
end

function toSuperscript(n::Integer)
    exps = ['⁰', '¹', '²', '³', '⁴', '⁵', '⁶', '⁷', '⁸', '⁹']
    dig = reverse(digits(n))
    return join([exps[d+1] for d in dig])
end


text:
# This is coord_conversion.jl
#
# This file contains functions to convert from geographical coordinates to Cartesian coordinates and vice versa.
#
# Author: Marcel Thielmann, 05/2021

function ConvertGeo2Cart(inputdata::GeoData)

     if cmp(inputdata.depth.unit,"km")
        rearth = 3671
     elseif cmp(inputdata.depth.unit,"m")
        rearth = 3671e3
     else
        warn("No proper depth unit is given (km or m), using m")
        rearth = 3671e3
     end

    lon = inputdata.lon.values
    lat = inputdata.lat.values
    depth = inputdata.depth.values

    R = rearth+inputdata.depth.values # radius from the center of the earth

    # compute Cartesian coordinates and assign them to a ValueList variable
    X = ValueList("x",inputdata.depth.unit,R.*cosd(inputdata.lon.values).*cosd(inputdata.lat.values))
    Y = ValueList("y",inputdata.depth.unit,R.*sind(inputdata.lon.values).*cosd(inputdata.lat.values))
    Z = ValueList("z",inputdata.depth.unit,R.*sind(inputdata.lat.values))

    # assign all data to the respective struct
    convdata = ParaviewData(X,Y,Z,inputdata.values)

   return convdata

end