text:
using Core.Compiler: MethodInstance, IncrementalCompact, insert_node_here!,
    userefs, SlotNumber, IRCode, compute_basic_blocks, _methods_by_ftype,
    retrieve_code_info, CodeInfo, SSAValue, finish, complete, non_dce_finish!,
    GotoNode, GotoIfNot, block_for_inst, ReturnNode, Argument, compact!,
    OldSSAValue, NewvarNode, quoted

using Base.Meta

cname(nc, N, name) = Symbol(string("∂⃖", superscript(N), subscript(nc), name))

using Core.Compiler: construct_domtree, scan_slot_def_use, construct_ssa!,
    NewInstruction, effect_free, CFG, BasicBlock, bbidxiter, PhiNode,
    Instruction, StmtRange, cfg_insert_edge!, insert_node!,
    non_effect_free, cfg_delete_edge!, domsort_ssa!

struct ∂ϕNode; end

struct BBEnv
    ctx_obj::Any
    bb_start_idx::Int
end

const debug = false

function check_back(nargs, covecs, msg)
    if length(covecs) != nargs
        error("Reverse for stmt `$msg` returned incorrect number of covectors (`$(typeof(covecs))`)")
    end
end

function new_to_regular(@nospecialize(stmt))
    urs = userefs(stmt)
    for op in urs
        val = op[]
        if isa(val, NewSSAValue)
            op[] = SSAValue(val.id)
        end
    end
    return urs[]
end

function expand_switch(code::Vector{Any}, bb_ranges::Vector{UnitRange{Int}}, slot_map)
    renumber = Vector{SSAValue}(undef, length(code))
    new_code = Vector{Any}()

    for val in values(slot_map)
        push!(new_code, Expr(:(=), val, ChainRulesCore.ZeroTangent()))
    end

    # TODO: This is a terrible data structure for this
    phi_rewrites = Dict{Pair{Int, Int}, Int}()

    # First expand switches into sequences of branches.
    # N.B: Code here isn't necessarily in domsort order, so we
    # must do the expansion and statement rewriting in two passes
    for (bb, range) in enumerate(bb_ranges)
        for i in range
            stmt = code[i]
            renumber[i] = SSAValue(length(new_code)+1)
            if isexpr(stmt, :switch)
                cond = stmt.args[1]
                labels = stmt.args[2]
                dests = stmt.args[3]
                for (label, dest) in zip(labels[1:end-1], dests[1:end-1])
                    push!(new_code, Expr(:call, !=, cond, label))
                    comp = NewSSAValue(length(new_code))
                    push!(new_code, GotoIfNot(comp, dest))
                    phi_rewrites[bb=>dest] = length(new_code)
                end
                push!(new_code, GotoNode(dests[end]))
                phi_rewrites[bb=>dests[end]] = length(new_code)
            else
                push!(new_code, stmt)
                if isa(stmt, GotoNode)
                    phi_rewrites[bb=>stmt.label] = length(new_code)
                elseif isa(stmt, GotoIfNot)
                    phi_rewrites[bb=>stmt.dest] = length(new_code)
                    phi_rewrites[bb=>bb+1] = length(new_code)
                else i == last(range)
                    phi_rewrites[bb=>bb+1] = length(new_code)
                end
                if isa(stmt, PhiNode)
                    # just temporarily remember the old bb here. This is a terrible hack, but oh well
                    push!(stmt.edges, bb)
                end
            end
        end
    end

    # Now rewrite branch targets back to statement indexing
    for i = 1:length(new_code)
        stmt = new_code[i]
        stmt = Core.Compiler.renumber_ssa!(stmt, renumber)
        stmt = new_to_regular(stmt)
        if isa(stmt, GotoNode)
            stmt = GotoNode(renumber[first(bb_ranges[stmt.label])].id)
        elseif isa(stmt, GotoIfNot)
            stmt = GotoIfNot(stmt.cond, renumber[first(bb_ranges[stmt.dest])].id)
        elseif isa(stmt, PhiNode)
            old_phi_bb = pop!(stmt.edges)
            stmt = PhiNode(map(old_pred->Int32(phi_rewrites[old_pred=>old_phi_bb]), stmt.edges),
                stmt.values)
        end
        new_code[i] = stmt
    end

    new_code
end

include("compiler_utils.jl")
include("hacks.jl")

struct PendingCtx
    list::Vector{Int}
end
PendingCtx() = PendingCtx(Int[])

# Split critical edges
# This is absolutely terrible, we really need better tools for this in
# Base
function split_critical_edges!(ir)
    cfg = ir.cfg
    blocks_to_split = Int[]
    edges_to_split = Pair{Int, Int}[]
    for (bb, block) in enumerate(cfg.blocks)
        length(block.preds) <= 1 && continue
        for pred in block.preds
            if length(cfg.blocks[pred].succs) > 1
                if pred+1 == bb
                    # Splitting a fallthrough edge
                    push!(blocks_to_split, bb)
                else
                    push!(edges_to_split, pred=>bb)
                end
            end
        end
    end

    if length(edges_to_split) == 0 && length(blocks_to_split) == 0
        return ir
    end

    for (pred, bb) in edges_to_split
        push!(ir, NewInstruction(GotoNode(bb)))
        push!(ir.cfg, BasicBlock(StmtRange(length(ir.stmts), length(ir.stmts))))
        new_bb = length(ir.cfg.blocks)
        gin = ir.stmts[last(cfg.blocks[pred].stmts)][:inst]
        @assert isa(gin, GotoIfNot)
        ir.stmts[last(cfg.blocks[pred].stmts)][:inst] =
            GotoIfNot(gin.cond, new_bb)
        cfg_delete_edge!(cfg, pred, bb)
        cfg_insert_edge!(cfg, pred, new_bb)
        cfg_insert_edge!(cfg, new_bb, bb)

        for i in cfg.blocks[bb].stmts
            stmt = ir.stmts[i][:inst]
            if isa(stmt, PhiNode)
                map!(stmt.edges, stmt.edges) do edge
                    edge == pred ? new_bb : edge
                end
            elseif stmt !== nothing
                break
            end
        end
    end

    for bb in blocks_to_split
        insert_node!(ir, cfg.blocks[bb].stmts.start,
            non_effect_free(NewInstruction(Expr(:new_bb_marker, bb))))
    end

    ir = compact!(ir)
    cfg = ir.cfg

    bb_rename_offset = Int[1 for _ in 1:length(cfg.blocks)]
    # Now expand basic blocks
    ninserted = 0
    i = 1
    while i < length(ir.stmts)
        if isexpr(ir.stmts[i][:inst], :new_bb_marker)
            bb = ir.stmts[i][:inst].args[1]
            ir.stmts[i][:inst] = nothing
            bbnew = bb + ninserted
            insert!(cfg.blocks, bbnew, BasicBlock(i:i))
            bb_rename_offset[bb] += 1
            bblock = cfg.blocks[bbnew+1]
            cfg.blocks[bbnew+1] = BasicBlock((i+1):last(bblock.stmts),
                bblock.preds, bblock.succs)
            i += 1
            while i <= last(bblock.stmts)
                stmt = ir.stmts[i][:inst]
                i += 1
                if isa(stmt, PhiNode)
                    map!(stmt.edges, stmt.edges) do edge
                        edge == bb-1 ? -bbnew : edge
                    end
                elseif stmt !== nothing
                    break
                end
            end
            continue
        end
        i += 1
    end
    bb_rename_mapping = cumsum(bb_rename_offset)
    # Repair CFG
    for i in 1:length(cfg.blocks)
        bb = cfg.blocks[i]
        cfg.blocks[i] = BasicBlock(bb.stmts, map(x->bb_rename_mapping[x], bb.preds),
            map(x->bb_rename_mapping[x], bb.succs))
    end
    for block in blocks_to_split
        cfg_delete_edge!(cfg, bb_rename_mapping[block-1], bb_rename_mapping[block])
        cfg_insert_edge!(cfg, bb_rename_mapping[block-1], bb_rename_mapping[block]-1)
        cfg_insert_edge!(cfg, bb_rename_mapping[block]-1, bb_rename_mapping[block])
    end

    # Repair IR
    for i in 1:length(ir.stmts)
        stmt = ir.stmts[i][:inst]
        isa(stmt, Union{GotoNode, PhiNode, GotoIfNot}) || continue
        if isa(stmt, GotoNode)
            ir.stmts[i][:inst] = GotoNode(bb_rename_mapping[stmt.label])
        elseif isa(stmt, PhiNode)
            map!(stmt.edges, stmt.edges) do edge
                edge < 0 && return -edge
                return bb_rename_mapping[edge]
            end
        else
            @assert isa(stmt, GotoIfNot)
            ir.stmts[i][:inst] = GotoIfNot(stmt.cond, bb_rename_mapping[stmt.dest])
        end
    end

    ir′ = compact!(domsort_ssa!(ir, construct_domtree(ir.cfg.blocks)))

    return ir′
end

Base.iterate(c::IncrementalCompact, args...) = Core.Compiler.iterate(c, args...)
Base.iterate(p::Core.Compiler.Pair, args...) = Core.Compiler.iterate(p, args...)
Base.iterate(urs::Core.Compiler.UseRefIterator, args...) = Core.Compiler.iterate(urs, args...)
Base.iterate(x::Core.Compiler.BBIdxIter, args...) = Core.Compiler.iterate(x, args...)
Base.getindex(urs::Core.Compiler.UseRefIterator, args...) = Core.Compiler.getindex(urs, args...)
Base.getindex(urs::Core.Compiler.UseRef, args...) = Core.Compiler.getindex(urs, args...)
Base.getindex(c::Core.Compiler.IncrementalCompact, args...) = Core.Compiler.getindex(c, args...)
Base.setindex!(c::Core.Compiler.IncrementalCompact, args...) = Core.Compiler.setindex!(c, args...)
Base.setindex!(urs::Core.Compiler.UseRef, args...) = Core.Compiler.setindex!(urs, args...)
function transform!(ci, meth, nargs, sparams, N)
    n_closures = 2^N - 1

    code = ci.code
    cfg = compute_basic_blocks(code)
    slotnames = Symbol[Symbol("#self#"), :args, ci.slotnames...]
    slotflags = UInt8[(0x00 for i = 1:2)..., ci.slotflags...]
    slottypes = UInt8[(0x00 for i = 1:2)..., ci.slotflags...]

    ir = IRCode(Core.Compiler.InstructionStream(code, Any[],
        Any[nothing for i = 1:length(code)],
        ci.codelocs, UInt8[0 for i = 1:length(code)]), cfg, Core.LineInfoNode[ci.linetable...],
        Any[Any for i = 1:2], Any[], Any[sparams...])

    # SSA conversion
    domtree = construct_domtree(ir.cfg.blocks)
    defuse_insts = scan_slot_def_use(VERSION >= v"1.8.0-DEV.267" ? Int(meth.nargs) : meth.nargs-1, ci, ir.stmts.inst)
    ci.ssavaluetypes = Any[Any for i = 1:ci.ssavaluetypes]
    if VERSION >= v"1.8.0-DEV.267"
        ir = construct_ssa!(ci, ir, domtree, defuse_insts, Any[Any for i = 1:length(slotnames)])
    else
        ir = construct_ssa!(ci, ir, domtree, defuse_insts, nargs, Any[Any for i = 1:length(slotnames)])
    end
    ir = compact!(ir)
    cfg = ir.cfg

    # If we have more than one basic block, canonicalize by creating a single
    # return node in the last basic block.

    if length(cfg.blocks) != 1
        ϕ = PhiNode()

        bb_start = length(ir.stmts)+1
        push!(ir, NewInstruction(ϕ))
        push!(ir, NewInstruction(ReturnNode(SSAValue(length(ir.stmts)))))
        push!(ir.cfg, BasicBlock(StmtRange(bb_start, length(ir.stmts))))
        new_bb_idx = length(ir.cfg.blocks)

        for (bb, i) in bbidxiter(ir)
            bb == new_bb_idx && break
            stmt = ir.stmts[i].inst
            if isa(stmt, ReturnNode)
                push!(ϕ.edges, bb)
                if !isa(stmt.val, SSAValue)
                    push!(ϕ.values, insert_node!(ir, i,
                        non_effect_free(NewInstruction(stmt.val))))
                else
                    push!(ϕ.values, stmt.val)
                end
                ir[i] = NewInstruction(GotoNode(new_bb_idx))
                cfg_insert_edge!(cfg, bb, new_bb_idx)
            end
        end

        ir = compact!(ir)
        ir = split_critical_edges!(ir)

        # If that resulted in the return not being the last block, fix that now.
        # We keep things simple this way, such that the basic blocks in the
        # forward and reverse are simply inverses of each other (i.e. the
        # exist block needs to be last, since the entry block needs to be first
        # in the reverse pass).

        if !isa(ir.stmts[end][:inst], ReturnNode)
            new_bb_idx = length(ir.cfg.blocks)+1
            for (bb, i) in bbidxiter(ir)
                stmt = ir.stmts[i][:inst]
                if isa(stmt, ReturnNode)
                    ir[i] = NewInstruction(GotoNode(new_bb_idx))
                    push!(ir, NewInstruction(stmt))
                    push!(ir.cfg, BasicBlock(StmtRange(length(ir.stmts), length(ir.stmts))))
                    cfg_insert_edge!(ir.cfg, bb, new_bb_idx)
                    break
                end
            end
        end

        cfg = ir.cfg

        # Now add a special control flow marker to every basic block
        # TODO: This lowering of control flow is extremely simplistic.
        # It needs to be improved in the very near future, but let's get
        # something working for now.
        for block in cfg.blocks
            if length(block.preds) != 0
                insert_node!(ir, block.stmts.start,
                    non_effect_free(NewInstruction(Expr(:phi_placeholder, copy(block.preds)))))
            end
        end

        ir = compact!(ir)
        cfg = ir.cfg
    end

    orig_bb_ranges = [first(bb.stmts):last(bb.stmts) for bb in cfg.blocks]

    # Special case: If there's only 1 `nothing` statement in the first BB,
    # we run the risk that compact skips it below. In that case, replace that
    # statement by a GotoNode to prevent that.
    if length(orig_bb_ranges[1]) == 1 && ir.stmts[1][:inst] === nothing
        ir.stmts[1][:inst] = GotoNode(2)
    end

    revs = Any[Any[nothing for i = 1:length(ir.stmts)] for i = 1:n_closures]
    opaque_cis = map(1:n_closures) do nc
        code = Any[] # Will be filled in later
        opaque_ci = ccall(:jl_new_code_info_uninit, Ref{CodeInfo}, ())
        opaque_ci.code = code
        if nc % 2 == 1
            opaque_ci.slotnames = Symbol[Symbol("#self#"), :Δ]
            opaque_ci.slotflags = UInt8[0, 0]
        else
            opaque_ci.slotnames = [Symbol("#oc#"), ci.slotnames...]
            opaque_ci.slotflags = UInt8[0, ci.slotflags...]
        end
        opaque_ci.linetable = Core.LineInfoNode[ci.linetable[1]]
        opaque_ci.inferred = false
        opaque_ci
    end

    nfixedargs = meth.isva ? meth.nargs - 1 : meth.nargs
    meth.isva || @assert nfixedargs == nargs+1

    extra_slotnames = Symbol[]
    extra_slotflags = UInt8[]

    # First go through and assign an accumulation slot to every used SSAValue/Argument
    has_cfg = length(cfg.blocks) != 1
    if has_cfg
        slot_map = Dict{Union{SSAValue, Argument}, SlotNumber}()
        phi_tmp_slot_map = Dict{Int, SlotNumber}()
        phi_uses = Dict{Int, Vector{Pair{SlotNumber, SlotNumber}}}()
        for (bb, i) in Iterators.reverse(bbidxiter(ir))
            first_bb_idx = cfg.blocks[bb].stmts.start
            stmt = ir.stmts[i][:inst]
            for urs in userefs(stmt)
                val = urs[]
                (isa(val, SSAValue) || isa(val, Argument)) || continue
                if !haskey(slot_map, val)
                    sn′ = SlotNumber(length(extra_slotnames) + 3)
                    push!(extra_slotnames, Symbol(string("for_", val)))
                    push!(extra_slotflags, UInt8(0))
                    slot_map[val] = sn′
                end
            end
            if isa(stmt, PhiNode)
                sn′ = SlotNumber(length(extra_slotnames) + 3)
                push!(extra_slotnames, Symbol(string("phi_temp_", i)))
                push!(extra_slotflags, UInt8(0))
                phi_tmp_slot_map[i] = sn′

                for use_n in 1:length(stmt.edges)
                    edge = stmt.edges[use_n]
                    # This assume the absence of critical edges
                    @assert length(cfg.blocks[edge].succs) == 1
                    val = stmt.values[use_n]
                    (isa(val, SSAValue) || isa(val, Argument)) || continue
                    push!(get!(phi_uses, edge, Vector{Pair{SlotNumber, SlotNumber}}()),
                        sn′=>slot_map[val])
                end
            end
        end
    end

    # TODO: Can we use the same method for each 2nd order of the transform
    # (except the last and the first one)
    for nc = 1:2:n_closures
        arg_accums = Union{Nothing, Vector{Any}}[nothing for i = 1:(meth.nargs)]
        accums = Union{Nothing, Vector{Any}}[nothing for i = 1:length(ir.stmts)]

        opaque_ci = opaque_cis[nc]
        code = opaque_ci.code

        function insert_node_rev!(node)
            push!(code, node)
            SSAValue(length(code))
        end

        if has_cfg
            append!(opaque_ci.slotnames, extra_slotnames)
            append!(opaque_ci.slotflags, extra_slotflags)
        end

        is_accumable(val) = isa(val, SSAValue) || isa(val, Argument)
        function accum!(val, accumulant)
            is_accumable(val) || return
            if !has_cfg
                if isa(val, SSAValue)
                    id = val.id
                    if isa(accums[id], Nothing)
                        accums[id] = Any[]
                    end
                    push!(accums[id], accumulant)
                elseif isa(val, Argument)
                    id = val.n
                    if isa(arg_accums[id], Nothing)
                        arg_accums[id] = Any[]
                    end
                    push!(arg_accums[id], accumulant)
                end
            else
                sn = slot_map[val]
                accumed = insert_node_rev!(Expr(:call, accum, sn, accumulant))
                insert_node_rev!(Expr(:(=), sn, accumed))
            end
        end

        function do_accum(for_val)
            if !has_cfg
                this_accums = isa(for_val, SSAValue) ? accums[for_val.id] :
                    arg_accums[for_val.n]
                if this_accums === nothing || isempty(this_accums)
                    return ChainRulesCore.ZeroTangent()
                elseif length(this_accums) == 1
                    return this_accums[]
                else
                    return insert_node_rev!(Expr(:call, accum, this_accums...))
                end
            else
                return get(slot_map, for_val, ChainRulesCore.ZeroTangent())
            end
        end

        to_back_bb(i) = length(cfg.blocks) - i + 1

        current_env = nothing
        ctx_map = Dict{Int, Any}()
        function retrieve_ctx_obj(current_env, i)
            if current_env === nothing
                return insert_node_rev!(Expr(:call, getfield, Argument(1),
                    i - first(orig_bb_ranges[end]) + 1))
            elseif isa(current_env, BBEnv)
                bbidx = i - current_env.bb_start_idx + 1
                @assert bbidx > 0
                return insert_node_rev!(Expr(:call, getfield, current_env.ctx_obj,
                    bbidx))
            end
            error()
        end

        function access_ctx_map(dest)
            if !haskey(ctx_map, dest)
                ctx_map[dest] = PendingCtx()
            end
            val = ctx_map[dest]
            if isa(val, PendingCtx)
                rval = insert_node_rev!(error)
                push!(val.list, rval.id)
                return rval
            end
            return val
        end

        bb_ranges = Vector{UnitRange{Int}}(undef, length(cfg.blocks))

        for (bb, i) in Iterators.reverse(bbidxiter(ir))
            first_bb_idx = cfg.blocks[bb].stmts.start
            last_bb_idx = cfg.blocks[bb].stmts.stop
            stmt = ir.stmts[i][:inst]

            if has_cfg && i == last_bb_idx
                if haskey(phi_uses, bb)
                    for (accumulant, accumulator) in phi_uses[bb]
                        accumed = insert_node_rev!(Expr(:call, accum, accumulator, accumulant))
                        insert_node_rev!(Expr(:(=), accumulator, accumed))
                    end
                end
                if !isa(stmt, Union{GotoNode, GotoIfNot, ReturnNode})
                    current_env = BBEnv(access_ctx_map(bb+1),
                        first(ir.cfg.blocks[bb].stmts))
                end
            end

            if isa(stmt, Core.ReturnNode)
                accum!(stmt.val, Argument(2))
                current_env = nothing
            elseif isexpr(stmt, :call)
                Δ = do_accum(SSAValue(i))
                callee = retrieve_ctx_obj(current_env, i)
                vecs = call = insert_node_rev!(Expr(:call, callee, Δ))
                if nc != n_closures
                    vecs = insert_node_rev!(Expr(:call, getfield, call, 1))
                    revs[nc+1][i] = insert_node_rev!(Expr(:call, getfield, call, 2))
                end
                if debug
                    insert_node_rev!(Expr(:call, check_back, length(stmt.args), vecs, string(stmt)))
                end
                for (j, arg) in enumerate(stmt.args)
                    if is_accumable(arg)
                        accum!(arg, insert_node_rev!(Expr(:call, getfield, vecs, j)))
                    end
                end
            elseif isexpr(stmt, :new)
                Δ = do_accum(SSAValue(i))
                newT = retrieve_ctx_obj(current_env, i)
                if nc != n_closures
                    revs[nc+1][i] = newT
                end
                # No gradient accumulation for zero-argument structs
                if length(stmt.args) != 1
                    # TODO: Use newT here?
                    #canon = insert_node_rev!(Expr(:call, ChainRulesCore.canonicalize, Δ))
                    #nt = insert_node_rev!(Expr(:call, ChainRulesCore.backing, canon))
                    nt = Δ
                    for (j, arg) in enumerate(stmt.args)
                        j == 1 && continue
                        if is_accumable(arg)
                            accum!(arg, insert_node_rev!(Expr(:call, lifted_getfield, nt, j-1)))
                        end
                    end
                end
            elseif isexpr(stmt, :splatnew)
                @assert length(stmt.args) == 2
                Δ = do_accum(SSAValue(i))
                newT = retrieve_ctx_obj(current_env, i)
                if nc != n_closures
                    revs[nc+1][i] = newT
                end
                #canon = insert_node_rev!(Expr(:call, ChainRulesCore.canonicalize, Δ))
                #nt = insert_node_rev!(Expr(:call, ChainRulesCore.backing, canon))
                nt = Δ
                arg = stmt.args[2]
                if is_accumable(stmt.args[2])
                    accum!(stmt.args[2], nt)
                end
            elseif isa(stmt, GlobalRef) || isexpr(stmt, :static_parameter) || isexpr(stmt, :throw_undef_if_not) || isexpr(stmt, :loopinfo)
                # We drop gradients for globals and static parameters
            elseif isexpr(stmt, :inbounds)
                # Nothing to do
            elseif isa(stmt, PhiNode)
                Δ = do_accum(SSAValue(i))
                @assert length(ir.cfg.blocks[bb].preds) >= 1
                insert_node_rev!(Expr(:(=), phi_tmp_slot_map[i], Δ))
            elseif isa(stmt, GotoNode)
                current_env = BBEnv(access_ctx_map(stmt.label),
                    first(ir.cfg.blocks[bb].stmts))
            elseif isa(stmt, GotoIfNot)
                current_env = BBEnv(insert_node_rev!(PhiNode(Int32.(map(to_back_bb, Int32[bb+1, stmt.dest])),
                    Any[access_ctx_map(bb+1),
                        access_ctx_map(stmt.dest)])), first(ir.cfg.blocks[bb].stmts))
            elseif isexpr(stmt, :phi_placeholder)
                @assert i == first_bb_idx
                tup = retrieve_ctx_obj(current_env, i)
                if length(stmt.args[1]) > 1
                    branch = insert_node_rev!(Expr(:call, getfield, tup, 1))
                    ctx = insert_node_rev!(Expr(:call, getfield, tup, 2))
                    insert_node_rev!(Expr(:switch, branch, collect(1:length(stmt.args[1])),
                        map(to_back_bb, stmt.args[1])))
                else
                    ctx = tup
                    insert_node_rev!(GotoNode(to_back_bb(stmt.args[1][1])))
                end
                if haskey(ctx_map, bb)
                    for item in (ctx_map[bb]::PendingCtx).list
                        code[item] = ctx
                    end
                end
                ctx_map[bb] = ctx
            elseif isa(stmt, Nothing)
                # Nothing to do
            else
                error((N, meth, stmt))
            end
            if has_cfg && haskey(slot_map, SSAValue(i))
                # Reset accumulator slots
                insert_node_rev!(Expr(:(=), slot_map[SSAValue(i)], ChainRulesCore.ZeroTangent()))
            end
            if i == first_bb_idx && bb != 1
                @assert isexpr(stmt, :phi_placeholder)
            end
            if i == first_bb_idx
                back_bb = to_back_bb(bb)
                bb_ranges[back_bb] = (back_bb == 1 ? 1 : (last(bb_ranges[back_bb - 1])+1)):length(code)
            end
        end

        ret_tuple = arg_tuple = insert_node_rev!(Expr(:call, tuple, [do_accum(Argument(i)) for i = 1:(nfixedargs)]...))
        if meth.isva
            ret_tuple = arg_tuple = insert_node_rev!(Expr(:call, Core._apply_iterate, iterate, Core.tuple, arg_tuple,
                do_accum(Argument(nfixedargs+1))))
        end
        if nc != n_closures
            lno = LineNumberNode(1, :none)
            next_oc = insert_node_rev!(Expr(:new_opaque_closure, Tuple{(Any for i = 1:nargs+1)...}, meth.isva, Union{}, Any,
                Expr(:opaque_closure_method, cname(nc+1, N, meth.name), Int(meth.nargs), lno, opaque_cis[nc+1]), revs[nc+1]...))
            ret_tuple = insert_node_rev!(Expr(:call, tuple, arg_tuple, next_oc))
        end
        insert_node_rev!(Core.ReturnNode(ret_tuple))
        bb_ranges[end] = first(bb_ranges[end]):length(code)

        if has_cfg
            code = opaque_ci.code = expand_switch(code, bb_ranges, slot_map)
        end

        opaque_ci.codelocs = Int32[0 for i=1:length(code)]
        opaque_ci.ssavaluetypes = length(code)
        opaque_ci.ssaflags = UInt8[0 for i=1:length(code)]
    end

    for nc = 2:2:n_closures
        fwds = Any[nothing for i = 1:length(ir.stmts)]

        opaque_ci = opaque_cis[nc]
        code = opaque_ci.code
        function insert_node_here!(node)
            push!(code, node)
            SSAValue(length(code))
        end

        for i in 1:length(ir.stmts)
            stmt = ir.stmts[i][:inst]
            if isa(stmt, Expr)
                stmt = copy(stmt)
            end
            urs = userefs(stmt)
            for op in urs
                val = op[]
                if isa(val, Argument)
                    op[] = Argument(val.n + 1)
                elseif isa(val, SSAValue)
                    op[] = fwds[val.id]
                else
                    op[] = ZeroTangent()
                end
            end
            stmt = urs[]

            if isexpr(stmt, :call)
                callee = insert_node_here!(Expr(:call, getfield, Argument(1), i))
                pushfirst!(stmt.args, callee)
                call = insert_node_here!(stmt)
                prim = insert_node_here!(Expr(:call, getfield, call, 1))
                dual = insert_node_here!(Expr(:call, getfield, call, 2))
                fwds[i] = prim
                revs[nc+1][i] = dual
            elseif isa(stmt, ReturnNode)
                lno = LineNumberNode(1, :none)
                next_oc = insert_node_here!(Expr(:new_opaque_closure, Tuple{Any}, false, Union{}, Any,
                    Expr(:opaque_closure_method, cname(nc+1, N, meth.name), 1, lno, opaque_cis[nc + 1]), revs[nc+1]...))
                ret_tup = insert_node_here!(Expr(:call, tuple, stmt.val, next_oc))
                insert_node_here!(ReturnNode(ret_tup))
            elseif isexpr(stmt, :new)
                newT = insert_node_here!(Expr(:call, getfield, Argument(1), i))
                revs[nc+1][i] = newT
                fnames = insert_node_here!(Expr(:call, fieldnames, newT))
                nt = insert_node_here!(Expr(:call, Core.apply_type, NamedTuple, fnames))
                stmt.head = :call
                stmt.args[1] = Core.tuple
                thetup = insert_node_here!(stmt)
                thent = insert_node_here!(Expr(:call, nt, thetup))
                ntT = insert_node_here!(Expr(:call, typeof, thent))
                compT = insert_node_here!(Expr(:call, Core.apply_type, Tangent, newT, ntT))
                fwds[i] = insert_node_here!(Expr(:new, compT, thent))
            elseif isexpr(stmt, :splatnew)
                error()
            elseif isa(stmt, GlobalRef)
                fwds[i] = ZeroTangent()
            elseif isa(stmt, Union{GotoNode, GotoIfNot})
                return :(error("Control flow support not fully implemented yet for higher-order reverse mode (TODO)"))
            elseif !isa(stmt, Expr)
                @show stmt
                error()
            else
                fwds[i] = insert_node_here!(stmt)
            end
        end

        opaque_ci.codelocs = Int32[0 for i=1:length(code)]
        opaque_ci.ssavaluetypes = length(code)
        opaque_ci.ssaflags = UInt8[0 for i=1:length(code)]
    end

    # TODO: This is absolutely aweful, but the best we can do given the data structures we have
    has_terminator = [isa(ir.stmts[last(range)].inst, Union{GotoNode, GotoIfNot}) for range in orig_bb_ranges]

    compact = IncrementalCompact(ir)

    arg_mapping = Any[]
    for argno in 1:nfixedargs
        push!(arg_mapping, insert_node_here!(compact,
            NewInstruction(Expr(:call, getfield, Argument(2), argno), Any, Int32(0))))
    end

    if meth.isva
        # Extract the rest of the arguments and make a tuple out of them
        ssas = map((nfixedargs+1):(nargs+1)) do i
            insert_node_here!(compact,
                NewInstruction(Expr(:call, getfield, Argument(2), i), Any, Int32(0)))
        end
        push!(arg_mapping,
            insert_node_here!(compact,
                NewInstruction(Expr(:call, tuple, ssas...), Any, Int32(0))))
    end

    rev = revs[1]
    active_bb = 1

    phi_nodes = Any[PhiNode() for _ = 1:length(orig_bb_ranges)]

    for ((old_idx, idx), stmt) in compact
        # remap arguments
        urs = userefs(stmt)
        compact[idx] = nothing
        for op in urs
            val = op[]
            if isa(val, Argument)
                op[] = arg_mapping[val.n]
            end
            if isexpr(val, :static_parameter)
                op[] = quoted(sparams[val.args[1]])
            end
        end
        compact[idx] = stmt = urs[]
        # f(args...) -> ∂⃖{N}(args...)
        orig_stmt = stmt
        if isexpr(stmt, :(=))
            stmt = stmt.args[2]
        end
        if isexpr(stmt, :call)
            compact[idx] = Expr(:call, ∂⃖{N}(), stmt.args...)
            if isexpr(orig_stmt, :(=))
                orig_stmt.args[2] = stmt
                stmt = orig_stmt
            end
            compact.ssa_rename[compact.idx-1] = insert_node_here!(compact,
                NewInstruction(Expr(:call, getfield, SSAValue(idx), 1), Any, compact.result[idx][:line]),
                true)
            rev[old_idx] = insert_node_here!(compact,
                NewInstruction(Expr(:call, getfield, SSAValue(idx), 2), Any, compact.result[idx][:line]),
                true)
        elseif isexpr(stmt, :static_parameter)
            stmt = quoted(sparams[stmt.args[1]])
            if isexpr(orig_stmt, :(=))
                orig_stmt.args[2] = stmt
                stmt = orig_stmt
            end
            compact[idx] = stmt
        elseif isexpr(stmt, :new) || isexpr(stmt, :splatnew)
            rev[old_idx] = stmt.args[1]
        elseif isexpr(stmt, :phi_placeholder)
            compact[idx] = phi_nodes[active_bb]
            # TODO: This is a base julia bug
            push!(compact.late_fixup, idx)
            rev[old_idx] = SSAValue(idx)
        elseif isa(stmt, Core.ReturnNode)
            lno = LineNumberNode(1, :none)
            compact[idx] = Expr(:new_opaque_closure, Tuple{Any}, false, Union{}, Any,
                Expr(:opaque_closure_method, cname(1, N, meth.name), 1, lno, opaque_cis[1]), rev[orig_bb_ranges[end]]...)
            argty = insert_node_here!(compact,
                NewInstruction(Expr(:call, typeof, stmt.val), Any, compact.result[idx][:line]), true)
            applyty = insert_node_here!(compact,
                NewInstruction(Expr(:call, Core.apply_type, OpticBundle, argty), Any, compact.result[idx][:line]),
                true)
            retval = insert_node_here!(compact,
                NewInstruction(Expr(:new, applyty, stmt.val, SSAValue(idx)), Any, compact.result[idx][:line]),
                true)
            compact.ssa_rename[compact.idx-1] = insert_node_here!(compact,
                NewInstruction(Core.ReturnNode(retval), Any, compact.result[idx][:line]),
                true)
        end

        succs = cfg.blocks[active_bb].succs

        if old_idx == last(orig_bb_ranges[active_bb])
            if length(succs) != 0
                override = false
                if has_terminator[active_bb]
                    terminator = compact[idx]
                    compact[idx] = nothing
                    override = true
                end
                function terminator_insert_node!(node)
                    if override
                        compact[idx] = node.stmt
                        override = false
                        return SSAValue(idx)
                    else
                        return insert_node_here!(compact, node, true)
                    end
                end
                tup = terminator_insert_node!(
                    effect_free(NewInstruction(Expr(:call, tuple, rev[orig_bb_ranges[active_bb]]...), Any, Int32(0))))
                for succ in succs
                    preds = cfg.blocks[succ].preds
                    if length(preds) == 1
                        val = tup
                    else
                        selector = findfirst(==(active_bb), preds)
                        val = insert_node_here!(compact, effect_free(NewInstruction(Expr(:call, tuple, selector, tup), Any, Int32(0))), true)
                    end
                    pn = phi_nodes[succ]
                    push!(pn.edges, active_bb)
                    push!(pn.values, val)
                end
                if has_terminator[active_bb]
                    insert_node_here!(compact, NewInstruction(terminator, Any, Int32(0)), true)
                end
            end
            active_bb += 1
        end
    end

    non_dce_finish!(compact)
    ir = complete(compact)
    ir = compact!(ir)
    Core.Compiler.verify_ir(ir)

    Core.Compiler.replace_code_newstyle!(ci, ir, nargs+1)
    ci.ssavaluetypes = length(ci.code)
    ci.slotnames = slotnames
    ci.slotflags = slotflags
    ci.slottypes = slottypes

    ci
end


text:
module ForwardPopGenSimulations

# package code goes here

include("chromosomes.jl")
end # module

ext: jl


text:
using PartitionedArrays

function throw_assert(parts)
  nparts=length(parts)
  map_parts(parts) do part
      @assert rand(1:nparts) != part
  end
end

function test_exception(parts)
  throw_assert(parts)
end


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

import TowardsNotebook
TowardsNotebook.notebook()

max_stars_repo_name: Jblew/towards-julia-notebook

max_issues_repo_name: Jblew/towards-julia-notebook

max_forks_repo_name: Jblew/towards-julia-notebook


text:
@inline function energy_integral(bound_low::Real, bound_up::Real,
                                 norm::Real, slope::Real, ρ::Real)

    # energy integral (eq.21 M01)
    en = 4π * c_light * norm * bound_low^4 /  ρ
    energy = en * ( (bound_up/bound_low)^(4.0 - slope ) - 1.0 ) / ( 4.0 - slope )

    if ( 4.0 - slope_soft ) < slope < ( 4.0 + slope_soft )
        slope_var = (slope - 4.0)/slope_soft
        energy2 = en * log10(bound_up/bound_low)
        if slope_var != 0.0
            energy = energy * slope_var + energy2 * ( 1.0 - slope_var )
        else
            energy = energy2
        end
    end

    return energy
end

# function calculateCREnergyInCGS(CR_N, CR_S, CR_Cut, ρ;
#                                 pmin::Real=10.0, pmax::Real=1.e7,
#                                 mc::Real=0.0, SelectBin::Int=-1,
#                                 units::GadgetPhysical=GadgetPhysical(),
#                                 verbose=true)
#     # calculates engergy per bin in cgs units.

#     if mc == 0.0
#         error("Error! No mc specified! Please provide m_particle * c_light to continue!")
#         return 1
#     end

#     CR_E = 0.0
#     c_light = 2.9979e10
#     ρ *= units.m_unit / units.l_unit^3  # convert rho in cgs units

#     CR_N .*= 1.e20

#     Nbins = size(CR_N, 1)
#     bin_width = log10(pmax/pmin)/Nbins
#     bin = 10.0.^collect(log10(pmin):bin_width:log10(pmax)) .* mc
#     above_cut = findall( bin .> CR_Cut/mc )
#     if length(above_cut) > 0
#         bin[above_cut] .= CR_Cut/mc
#     end

#     if SelectBin == -1
#         CR_E = zeros(length(CR_N))
#         for i = 1:length(CR_N)
#             CR_E[i] = energy_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
#         end
#     elseif length(SelectBin) == 1
#         CR_E = energy_integral(bin[SelectBin], bin[SelectBin+1], CR_N[SelectBin], CR_S[SelectBin], ρ)
#     elseif length(SelectBin) > 1
#         CR_E = zeros(length(SelectBin))
#         if i ∈ SelectBin
#             CR_E[i-SelectBin[1]+1] = energy_integral(bin[i], bin[i+1], CR_N[i], CR_S[i], ρ)
#         end
#     end
#     return CR_E
# end


text:
# if the environment is not defined, use this
if ! haskey(ENV, "GENIE_ENV")
  ENV["GENIE_ENV"] = "dev"
end
if ! haskey(ENV, "NWORKERS") && in("s", ARGS)
  ENV["NWORKERS"] = 1
end

max_stars_repo_name: anirudhmurali/hello-julia-genie

max_issues_repo_name: anirudhmurali/hello-julia-genie

max_forks_repo_name: anirudhmurali/hello-julia-genie


text:
export MessageCreate,
    MessageUpdate,
    MessageDelete,
    MessageDeleteBulk,
    MessageReactionAdd,
    MessageReactionRemove,
    MessageReactionRemoveAll

"""
Sent when a [`Message`](@ref) is sent.
"""
struct MessageCreate <: AbstractEvent
    message::Message
end
@boilerplate MessageCreate :docs :mock
MessageCreate(; kwargs...) = MessageCreate(Message(; kwargs...))
MessageCreate(d::Dict{Symbol, Any}) = MessageCreate(; d...)

"""
Sent when a [`Message`](@ref) is updated.
"""
struct MessageUpdate <: AbstractEvent
    message::Message
end
@boilerplate MessageUpdate :docs :mock
MessageUpdate(; kwargs...) = MessageUpdate(Message(; kwargs...))
MessageUpdate(d::Dict{Symbol, Any}) = MessageUpdate(; d...)

"""
Sent when a [`Message`](@ref) is deleted.
"""
struct MessageDelete <: AbstractEvent
    id::Snowflake
    channel_id::Snowflake
    guild_id::Optional{Snowflake}
end
@boilerplate MessageDelete :constructors :docs :mock

"""
Sent when multiple [`Message`](@ref)s are deleted in bulk.
"""
struct MessageDeleteBulk <: AbstractEvent
    ids::Vector{Snowflake}
    channel_id::Snowflake
    guild_id::Optional{Snowflake}
end
@boilerplate MessageDeleteBulk :constructors :docs :mock

"""
Sent when a [`Reaction`](@ref) is added to a [`Message`](@ref).
"""
struct MessageReactionAdd <: AbstractEvent
    user_id::Snowflake
    channel_id::Snowflake
    message_id::Snowflake
    guild_id::Optional{Snowflake}
    emoji::Emoji
end
@boilerplate MessageReactionAdd :constructors :docs :mock

"""
Sent when a [`Reaction`](@ref) is removed from a [`Message`](@ref).
"""
struct MessageReactionRemove <: AbstractEvent
    user_id::Snowflake
    channel_id::Snowflake
    message_id::Snowflake
    guild_id::Optional{Snowflake}
    emoji::Emoji
end
@boilerplate MessageReactionRemove :constructors :docs :mock

"""
Sent when all [`Reaction`](@ref)s are removed from a [`Message`](@ref).
"""
struct MessageReactionRemoveAll <: AbstractEvent
    channel_id::Snowflake
    message_id::Snowflake
    guild_id::Optional{Snowflake}
end
@boilerplate MessageReactionRemoveAll :constructors :docs :mock


text:
import FrankWolfe
import LinearAlgebra

n = Int(1e6);

xpi = rand(n);
total = sum(xpi);
const xp = xpi ./ total;

const f(x) = norm(x - xp)^2

function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

function cf(x, xp)
    return @. norm(x - xp)^2
end

function cgrad(storage, x, xp)
    @. storage = 2 * (x - xp)
end

lmo_prob = FrankWolfe.ProbabilitySimplexOracle(1);
x0 = FrankWolfe.compute_extreme_point(lmo_prob, zeros(n));

FrankWolfe.benchmarkOracles(f, grad!, lmo_prob, n; k=100, T=Float64)

FrankWolfe.benchmarkOracles(
    x -> cf(x, xp),
    (storage, x) -> cgrad(storage, x, xp),
    lmo_prob,
    n;
    k=100,
    T=Float64,
)


text:
# ---
# title: 31. Next Permutation
# id: problem31
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/next-permutation/description/>
# hidden: true
# ---
#
# Implement **next permutation** , which rearranges numbers into the
# lexicographically next greater permutation of numbers.
#
# If such an arrangement is not possible, it must rearrange it as the lowest
# possible order (i.e., sorted in ascending order).
#
# The replacement must be **[in place](http://en.wikipedia.org/wiki/In-
# place_algorithm)** and use only constant extra memory.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,3]
#     Output: [1,3,2]
#
#
# **Example 2:**
#
#
#
#     Input: nums = [3,2,1]
#     Output: [1,2,3]
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,1,5]
#     Output: [1,5,1]
#
#
# **Example 4:**
#
#
#
#     Input: nums = [1]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 100`
#   * `0 <= nums[i] <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
struct PoincareInvariant2nd{DT,ND,NC,NV,ET,ΩT,ϑT}
    equ::ET
    ω::ΩT
    D²ϑ::ϑT
    Δt::DT
    nx::Int
    ny::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::Vector{DT}
    J::Vector{DT}
    K::Vector{DT}
    L::Vector{DT}
end

function PoincareInvariant2nd(f_equ::Function, f_surface::Function, ω::ΩT, D²ϑ::ϑT, Δt::TT, nd::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΩT,ϑT}

    if get_config(:verbosity) > 1
        println()
        println("Second Euler-Poincaré Integral Invariant")
        println("========================================")
        println()
        println(" nx    = ", nx)
        println(" ny    = ", ny)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute Chebyshev points
    c = points(Chebyshev(0..1)^2, nx*ny)

    # compute initial conditions
    q₀ = zeros(DT, (nd, length(c)))

    for i in 1:length(c)
        q₀[:,i] .= f_surface(c[i][1], c[i][2])
    end

    # initialise euation
    equ = f_equ(q₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I = zeros(DT, nt+1)
    J = zeros(DT, nt+1)
    K = zeros(DT, nt+1)
    L = zeros(DT, nt+1)

    # get size of coefficient and value vectors
    SC = Chebyshev(0..1)^2
    SU = Ultraspherical(1, 0..1)^2
    Dx = Derivative(SC, [1,0])
    Dy = Derivative(SC, [0,1])

    fq = Fun(Dx * Fun(SC, ApproxFun.transform(SC, q₀[1,:])), SU)
    nc = ncoefficients(fq)
    nv = length(values(fq))

    # initialise Poincare invariant
    PoincareInvariant2nd{DT,nd,nc,nv,typeof(equ),ΩT,ϑT}(equ, ω, D²ϑ, DT(Δt), nx, ny, ntime, nsave, nt, I, J, K, L)
end

function PoincareInvariant2nd(f_equ::Function, f_surface::Function, ω::ΩT, Δt::TT, nd::Int, nx::Int, ny::Int, ntime::Int, nsave::Int=1, DT=Float64) where {TT,ΩT}
    PoincareInvariant2nd(f_equ, f_surface, ω, (), Δt, nd, nx, ny, ntime, nsave, DT)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant2nd{DT}, sol::Solution) where {DT}

    local verbosity = get_config(:verbosity)

    # local q = permutedims(sol.q.d, [3,1,2])
    # if isdefined(sol, :p)
    #     local p = permutedims(sol.p.d, [3,1,2])
    # end
    # if isdefined(sol, :λ)
    #     local λ = permutedims(sol.λ.d, [3,1,2])
    # end

    verbosity ≤ 1 ? prog = Progress(size(sol.q.d,2), 5) : nothing

    for i in 1:size(sol.q.d,2)
        verbosity > 1 ? println("      it = ", i-1) : nothing
        pinv.I[i] = compute_noncanonical_invariant(pinv, sol.t.t[i], sol.q.d[:,i,:])
        verbosity > 1 ? println("           I_q = ", pinv.I[i], ",   ε_q = ", (pinv.I[i]-pinv.I[1])/pinv.I[1]) : nothing

        if isdefined(sol, :p)
            pinv.J[i] = compute_canonical_invariant(pinv, sol.q.d[:,i,:], sol.p.d[:,i,:])
            verbosity > 1 ? println("           I_p = ", pinv.J[i], ",   ε_p = ", (pinv.J[i]-pinv.J[1])/pinv.J[1]) : nothing
        end

        if isdefined(sol, :λ)
            pinv.K[i] = compute_noncanonical_correction(pinv, sol.t.t[i], sol.q.d[:,i,:], sol.λ.d[:,i,:])
            pinv.L[i] = pinv.I[i] - pinv.Δt^2 * pinv.K[i]
            verbosity > 1 ? println("           I_λ = ", pinv.L[i], ",   ε_λ = ", (pinv.L[i]-pinv.L[1])/pinv.L[1]) : nothing
            verbosity > 1 ? println("           K_λ = ", pinv.K[i]) : nothing
        end

        verbosity ≤ 1 ? next!(prog) : nothing
    end

    return (pinv.I, pinv.J, pinv.K, pinv.L)
end

function CommonFunctions.write_to_hdf5(pinv::PoincareInvariant2nd{DT}, sol::Solution, output_file::String) where {DT}
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5

        write(h5, "t", sol.t.t)
        write(h5, "I", pinv.I)

        isdefined(sol, :p) ? write(h5, "J", pinv.J) : nothing
        isdefined(sol, :λ) ? write(h5, "K", pinv.L) : nothing
        isdefined(sol, :λ) ? write(h5, "L", pinv.L) : nothing

    end
end


text:
module Covid19Analysis

using Logging, LoggingExtras

function main()
  Base.eval(Main, :(const UserApp = Covid19Analysis))

  include(joinpath("..", "genie.jl"))

  Base.eval(Main, :(const Genie = Covid19Analysis.Genie))
  Base.eval(Main, :(using Genie))
end; main()

end


text:
"""
Given a graph `g` and a rooted tree `t` within `g`, we can decompose `g` into cycles:

Every non-tree edge `(u,v)` of `g` defines a cycle, i.e. the
unique path `v --> u` within the tree `t` plus the edge `(u, v)` to close it.

This module contains functions to compute within the cycle space.
"""
module Cycle

import GraphIdx
import GraphIdx: IncidenceIndex, Weights, Graph
import GraphIdx.Tree: ChildrenIndex, root_node, lowest_common_ancestors, dfs_finish

"""
All kind of indices to represant a cycle basis corresponding to spanning tree
"""
mutable struct CycleBasis
    pi::Vector{Int}
    cidx::ChildrenIndex
    non_tree_edges::Vector{Pair{Int,Int}}
    non_tree_enum::Vector{Int}
    tree_enum::Vector{Int}
    non_tree_idx::IncidenceIndex
    lca::Vector{Int}
end

GraphIdx.num_nodes(cb::CycleBasis) =
    length(cb.tree_enum)

GraphIdx.num_edges(cb::CycleBasis) =
    GraphIdx.num_nodes(cb) + length(cb.non_tree_enum)

"""
    CycleBasis(graph)
    CyccleBasis(m, n)

If no tree is provided, just allocate the space.
"""
CycleBasis(g::Graph)::CycleBasis =
    CycleBasis(GraphIdx.num_edges(g), GraphIdx.num_nodes(g))

function CycleBasis(m::Integer, n::Integer)
    non_tree_edges = Vector{Pair{Int,Int}}()
    non_tree_enum = Vector{Int}()
    sizehint!(non_tree_edges, m - n + 1)
    sizehint!(non_tree_enum, m - n + 1)
    tree_enum = Vector{Int}(undef, n)
    return CycleBasis(
        [],                     # pi
        ChildrenIndex([], []),  # cidx
        non_tree_edges,
        non_tree_enum,
        tree_enum,
        IncidenceIndex([], []), # non_tree_idx
        [],                     # lca
    )
end

"""
    CycleBasis(graph, pi, [cidx])

Construct a cycle basis given a tree with parents `pi` and `cidx::ChildrenIndex`.
"""
function CycleBasis(g::Graph, pi::Vector{Int}, cidx::ChildrenIndex)::CycleBasis
    n = GraphIdx.num_nodes(g)
    m = GraphIdx.num_edges(g)
    non_tree = Vector{Pair{Int,Int}}()
    non_tree_enum = Vector{Int}()
    sizehint!(non_tree, m - n + 1)
    sizehint!(non_tree_enum, m - n + 1)
    tree_enum = Vector{Int}(undef, n)
    tree_enum[root_node(cidx)] = 0
    GraphIdx.enumerate_edges(g) do i::Int, u::Int, v::Int
        if pi[u] == v
            tree_enum[u] = i
        elseif pi[v] == u
            tree_enum[v] = i
        else
            push!(non_tree, u => v)
            push!(non_tree_enum, i)
        end
    end
    nidx = IncidenceIndex(n, non_tree)
    lca = lowest_common_ancestors(cidx, pi, nidx)
    CycleBasis(
        pi,
        cidx,
        non_tree,
        non_tree_enum,
        tree_enum,
        nidx,
        lca,
    )
end

CycleBasis(g::Graph, pi::Vector{Int}) =
    CycleBasis(g, pi, ChildrenIndex(pi))

"""
    enumerate_non_tree(func, cb)

For every cycle call function `func(ei, u, v, r)` whereby
- `ei` is the index of the corresponding non-tree edge `(u, v)`
- `r` is the root of the cycle, i.e. the lowest common
  ancestor of `u` and `v` within the tree

Notice that for every edge `(u, v)` the function is called twice:
With `(u, v)` and `(v, u)`.
"""
@inline function enumerate_non_tree(func::F, cb::CycleBasis) where {F<:Function}
    n = GraphIdx.num_nodes(cb)
    for v = 1:n
        for (u, i) in cb.non_tree_idx[v]
            r = cb.lca[i]
            func(i, u, v, r)
        end
    end
end

"""
    extract_rotate(cb, λ)

Return the weights
TODO: finish
"""
function extract_rotate(
    cb::CycleBasis, lam::Weights{F}
)::NTuple{2, Vector{F}} where {F}
    n = GraphIdx.num_nodes(cb)
    tlam = Vector{Float64}(undef, n)
    tlam0 = Vector{Float64}(undef, n)
    ldiff = zeros(Float64, n)
    stack = Vector{Int}()
    sizehint!(stack, n)

    for (i, ei) in enumerate(cb.tree_enum)
        tlam[i] = ei > 0 ? lam[ei] : NaN
    end
    enumerate_non_tree(cb) do i, _, v, r
        let lami = lam[cb.non_tree_enum[i]]
            ldiff[v] += lami
            ldiff[r] -= lami
        end
    end
    tlam0 .= tlam
    for v in dfs_finish(cb.pi, stack)
        tlam[v] += ldiff[v]
        ldiff[cb.pi[v]] += ldiff[v]
        ldiff[v] = 0
    end
    return tlam0, tlam
end

function GraphIdx.GraphViz.dot(
    graph::GraphIdx.Graph, cycles::CycleBasis, tree_label::F = i -> ""
) where {F<:Function}
    open(`dot -Tx11 -Kneato`, "w") do io::IO
        GraphIdx.GraphViz.dot(io, graph) do io::IO
            println(io, "{ edge [color=black]")
            for (i, v) in enumerate(cycles.pi)
                if i != v
                    println(io, i, " -> ", v, tree_label(i))
                end
            end
            println(io, "}")
            Cycle.enumerate_non_tree(cycles) do ei::Int, u::Int, v::Int, r::Int
                if u < v
                    local col = 1 + (ei % 8)
                    tstyle = " [ color=$(col) ]"
                    println(io, u, " -> ", v , " [style=dotted, color=", col, "]")
                    for i in [u, v]
                        while i != r
                            local pi = cycles.pi[i]
                            println(io, i, " -> ", pi, tstyle)
                            i = pi
                        end
                    end
                end
            end
        end
    end
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
using Distributions
using Iterators
using Gurobi
using JLD

using EmpiricalBayes

nreps = 50
true_dists = [MixtureModel([ Normal(-0.3,.5), Normal(1.05,.5)])]
ys = collect(range(0, stop=3, length=7))
minimax_σs = [0.1, 0.5, 1.0]
ns = [20_000; 100_000]

combs = []

for (true_dist, y, σ, n) in product(true_dists, ys, minimax_σs, ns)
    push!(combs, Dict(:dist=>true_dist, :y=>y, :σ=>σ, :n=>n) )
end

i = parse(Int64, ARGS[1])

comb = combs[i]

n = comb[:n]
σ = comb[:σ]
true_dist = comb[:dist]
y = comb[:y]

marginal_grid = collect(range(-10, stop=10, length=1001));
prior_grid = collect(range(-3, stop=3, length=51));
d_true = NormalConvolutionProblem(true_dist, marginal_grid)
f = BinnedMarginalDensity(d_true)
m = n
ds = MixingNormalConvolutionProblem(Normal, 0.5, prior_grid, marginal_grid)

sim_array = []

for l = 1:nreps
    @show l
     Xs = rand(d_true,m)
     tst = donoho_ci(Xs, marginal_grid, prior_grid, y, σ, d_true)
     tst2 = donoho_ci(Xs, marginal_grid, prior_grid, y, σ, NPMLE)
     push!(sim_array, (tst, tst2))
end

res = (i, comb, sim_array)

save("/scratch/users/ignat/sims/Apr5/mysim_$i.jld", "res", res)


text:
using Plots
using StatsPlots
using CSV
using Query
using DataFrames
using IncrementalPruning
using SARSOP
using POMDPSolve
using ContinuousObservationToyPOMDPs
using POMDPs
using POMDPPolicies
using POMDPModelTools
using LinearAlgebra
using QMDP

n = 200
maxwidth = 41
widthstep = 4
depth = 3
p_correct = 0.85

fname = joinpath(dirname(@__FILE__()), "data", "width_41_Mon_7_Oct_2019_16_55.csv")

# pgfplots()
pyplot()

color_d = Dict("wait" => "blue",
               "listen" => "red")
color_f(a) = get(color_d, string(a), "black")

label_d = Dict("wait" => "Wait",
               "listen" => "Listen",
              )
label_f(a) = label_d[string(a)]

results = CSV.read(fname)
# @show results

plotted = results |> @filter(_.solver=="POWSS") |> @filter(_.action in ("wait", "listen")) |> DataFrame

p = plot(xlabel="Width (C)", ylabel="Estimated Q Value", size=(400,300), legend=(0.6,0.45))
for a in ("wait", "listen")
    df = plotted |> @filter(_.action == a) |> DataFrame
    plot!(p, df[!,:width], df[!,:mean],
          ribbon=df[!,:std],
          color=color_f(a),
          label="POWSS "*label_f(a)
        )
end

dm = TimedDOTigerPOMDP(horizon=depth, p_correct=p_correct)
op  = solve(PruneSolver(), dm)
b = initialize_belief(updater(op), initialstate_distribution(dm))
optpairs = Dict(first(item)=>last(item) for item in zip(ordered_actions(dm), actionvalues(op, b)))

a = :wait
Sz = Vector{Set{Vector{Float64}}}(undef, length(observations(dm)))
for (zind, z) in enumerate(ordered_observations(dm))
    V = Set(IncrementalPruning.dpval(α,a,z,dm) for α in op.alphas)
    Sz[zind] = IncrementalPruning.filtervec(V, PruneSolver().optimizer_factory)
end
Sa = Set([IncrementalPruning.AlphaVec(α,a) for α in IncrementalPruning.incprune(Sz, PruneSolver().optimizer_factory)])

bvec = zeros(length(states(dm)))
bvec[stateindex(dm, (:left, 0))] = 0.5
bvec[stateindex(dm, (:right, 0))] = 0.5
@show optpairs[:wait] = maximum(dot(bvec, α.alpha) for α in Sa)

for a in [:wait, :listen]
    plot!(p, [1, maximum(results[!, :width])], fill(optpairs[a], 2),
          color=color_f(a),
          linestyle=:dash,
          label="Optimal "*label_f(a)
         )
end
savefig(p, "fig/powss_convergence_plot.pdf")

# POSS/QMDP
plotted = results |> @filter(_.solver=="POSS") |>  DataFrame

for a in ("wait", "listen")
    df = plotted |> @filter(_.action == a) |> DataFrame
    plot!(p, df[!,:width], df[!,:mean],
          color=color_f(a),
          label="POSS "*label_f(a),
          linestyle=:dot
        )
end
@show savefig(p, "fig/convergence_plot.pdf")


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using Documenter
using Documenter.Writers: HTMLWriter
using DocumenterTools.Themes: compile

# theme extension
for theme in ["light", "dark"]
    mktemp(@__DIR__) do path, io
        write(io, join([
            read(joinpath(HTMLWriter.ASSETS_THEMES, "documenter-$(theme).css"), String),
            read(joinpath(@__DIR__, "src/assets/theme-$(theme).css"), String)
        ], "\n"))
        compile(path, joinpath(@__DIR__, "src/assets/themes/documenter-$(theme).css"))
    end
end

# outline
contents = [
    "Idealni pripad" => "./idealni_pripad.md",
]

# make docs options
makedocs(;
    authors = "Petr Vala",
    sitename = "Jak se budi princezny",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        collapselevel = 1,
    ),
    pages = [
        "Home" => "index.md",
        "Obsah" => contents,
    ],
)

deploydocs(;
    repo="github.com/petrmvala/JakSeBudiPrincezny.git",
    devbranch = "main",
)

max_stars_repo_name: petrmvala/jak-se-budi-princezny

max_issues_repo_name: petrmvala/jak-se-budi-princezny

max_forks_repo_name: petrmvala/jak-se-budi-princezny


text:
"""
    SplitIntervalBridge{T}

The `SplitIntervalBridge` splits a constraint ``l ≤ ⟨a, x⟩ + α ≤ u`` into the constraints ``⟨a, x⟩ + α ≥ l`` and ``⟨a, x⟩ + α ≤ u``.
"""
struct SplitIntervalBridge{T, F<:MOI.AbstractScalarFunction} <: AbstractBridge
    lower::CI{F, MOI.GreaterThan{T}}
    upper::CI{F, MOI.LessThan{T}}
end
function SplitIntervalBridge{T, F}(model, f::F, s::MOI.Interval{T}) where {T, F}
    lower = MOI.add_constraint(model, f, MOI.GreaterThan(s.lower))
    upper = MOI.add_constraint(model, f, MOI.LessThan(s.upper))
    return SplitIntervalBridge{T, F}(lower, upper)
end

MOI.supports_constraint(::Type{SplitIntervalBridge{T}}, ::Type{<:MOI.AbstractScalarFunction}, ::Type{MOI.Interval{T}}) where T = true
function added_constraint_types(::Type{SplitIntervalBridge{T, F}}) where {T, F}
    return [(F, MOI.GreaterThan{T}), (F, MOI.LessThan{T})]
end
function concrete_bridge_type(::Type{<:SplitIntervalBridge},
                              F::Type{<:MOI.AbstractScalarFunction},
                              ::Type{MOI.Interval{T}}) where T
    return SplitIntervalBridge{T, F}
end

# Attributes, Bridge acting as an model
MOI.get(b::SplitIntervalBridge{T, F}, ::MOI.NumberOfConstraints{F, MOI.LessThan{T}}) where {T, F} = 1
MOI.get(b::SplitIntervalBridge{T, F}, ::MOI.NumberOfConstraints{F, MOI.GreaterThan{T}}) where {T, F} = 1
MOI.get(b::SplitIntervalBridge{T, F}, ::MOI.ListOfConstraintIndices{F, MOI.GreaterThan{T}}) where {T, F} = [b.lower]
MOI.get(b::SplitIntervalBridge{T, F}, ::MOI.ListOfConstraintIndices{F, MOI.LessThan{T}}) where {T, F} = [b.upper]

# Indices
function MOI.delete(model::MOI.ModelLike, c::SplitIntervalBridge)
    MOI.delete(model, c.lower)
    MOI.delete(model, c.upper)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal, c::SplitIntervalBridge)
    # lower and upper should give the same value
    return MOI.get(model, attr, c.lower)
end
function MOI.get(model::MOI.ModelLike, a::MOI.ConstraintDual, c::SplitIntervalBridge)
    # Should be nonnegative
    lower_dual = MOI.get(model, MOI.ConstraintDual(), c.lower)
    # Should be nonpositive
    upper_dual = MOI.get(model, MOI.ConstraintDual(), c.upper)
    return lower_dual > -upper_dual ? lower_dual : upper_dual
end

function MOI.get(model::MOI.ModelLike, ::MOI.ConstraintBasisStatus, c::SplitIntervalBridge)
    lower_stat = MOI.get(model, MOI.ConstraintBasisStatus(), c.lower)
    upper_stat = MOI.get(model, MOI.ConstraintBasisStatus(), c.upper)
    if lower_stat == MOI.NONBASIC_AT_LOWER
        Compat.@warn("GreaterThan constraints should not have basis status:" *
            " NONBASIC_AT_LOWER, instead use NONBASIC.")
    end
    if upper_stat == MOI.NONBASIC_AT_UPPER
        Compat.@warn("LessThan constraints should not have basis status:" *
            " NONBASIC_AT_UPPER, instead use NONBASIC.")
    end
    if lower_stat == MOI.NONBASIC
        return MOI.NONBASIC_AT_LOWER
    end
    if upper_stat == MOI.NONBASIC
        return MOI.NONBASIC_AT_UPPER
    end
    if lower_stat != upper_stat
        Compat.@warn("Basis status of lower (`$lower_stat`) and upper (`$upper_stat`) constraint are inconsistent," *
            " both should be basic or super basic.")
    end
    return lower_stat
end

# Constraints
function MOI.modify(model::MOI.ModelLike, c::SplitIntervalBridge, change::MOI.AbstractFunctionModification)
    MOI.modify(model, c.lower, change)
    MOI.modify(model, c.upper, change)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintFunction,
                  c::SplitIntervalBridge{T, F}, func::F) where {T, F}
    MOI.set(model, MOI.ConstraintFunction(), c.lower, func)
    MOI.set(model, MOI.ConstraintFunction(), c.upper, func)
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet, c::SplitIntervalBridge, change::MOI.Interval)
    MOI.set(model, MOI.ConstraintSet(), c.lower, MOI.GreaterThan(change.lower))
    MOI.set(model, MOI.ConstraintSet(), c.upper, MOI.LessThan(change.upper))
end


text:
import Flux: gpu
import Flux.Tracker: data

function build_encoder(seqLen::Integer, dims::Integer)
    positions = 1:seqLen
    columns = 1:dims
    factor = (1e5 ^ (2 / dims)) .^ columns
    factor = reshape(factor, size(factor, 1), 1)

    even = sin.(positions / factor[1:2:end, :])
    odd = cos.(positions / factor[2:2:end, :])

    encoder = zeros(Float32, seqLen, dims)
    encoder[:, 1:2:end] = even
    encoder[:, 2:2:end] = odd
    return encoder
end

struct PositionEncoder{T}
    E::AbstractArray{T, 2}
end

PositionEncoder(seqLen::Integer, dims::Integer) = PositionEncoder(data(gpu(build_encoder(seqLen, dims))))

function (m::PositionEncoder)(x::AbstractArray{T, 3}, mask::AbstractArray{T, 3}) where T
    maxLen = size(m.E, 1)
    x = x[1:min(size(x, 1), maxLen), :, :]
    return (m.E .* mask) .+ x
end


text:
function LocalOrbitalFrame(arg0::Frame, arg1::LOFType, arg2::PVCoordinatesProvider, arg3::JString)
    return LocalOrbitalFrame((Frame, LOFType, PVCoordinatesProvider, JString), arg0, arg1, arg2, arg3)
end


text:
module ClimaCore

include("interface.jl")
include("Utilities/Utilities.jl")
include("RecursiveApply/RecursiveApply.jl")
include("DataLayouts/DataLayouts.jl")
include("Geometry/Geometry.jl")
include("Domains/Domains.jl")
include("Meshes/Meshes.jl")
include("Topologies/Topologies.jl")
include("Spaces/Spaces.jl")
include("Fields/Fields.jl")
include("Operators/Operators.jl")

end # module


text:
function unit(x::AbstractVector, entry::Integer)
	unit(eltype(x), length(x), entry)
end

function unit(length::Integer, entry::Integer)
	unit(Float64, length, entry)
end

function unit(T::Type, length::Integer, entry::Integer)
	u = zeros(T, length)
	u[entry] = one(T)
	return u
end

@doc """
	eigval1(A::Matrix) -> Vector

Returns a unit eigenvector of `A` for the eigenvalue 1 if this eigenspace is 1D.
Otherwise, return an error.
"""->
function eigval1(A::DenseMatrix{Float64})
	E = eigfact(A)
	# Types are hardcoded; otherwise instability occurs
	isreal(E[:values]) || throw(DomainError())
	vals = E[:values]::Vector{Float64}
	vecs = E[:vectors]::Matrix{Float64}

	# Find index of eigenvalue 1
	eval1_index = Vector{Int64}()
	for i in 1:length(vals)
		if isapprox(vals[i], 1.0)
			push!(eval1_index, i)
		end
	end

	if isempty(eval1_index)
		error("1 is not an eigenvalue")
	elseif length(eval1_index) >= 2
		error("The eigenspace of 1 is more than 1D")
	end

	return vecs[:,eval1_index[]]
end

function van_moment(wavename::String)
	lowername = lowercase( wavename )
	WT.vanishingmoments( WT.eval(parse(lowername)) )::Int64
end

@doc """
	support(wavename) -> DaubSupport

The support of the Daubechies scaling function `wavename`.
"""->
function support(wavename::String)
	vm = van_moment(wavename)
	return DaubSupport(-vm+1, vm)
end

@doc """
	isinside(x, S::DaubSupport) -> Bool

Returns `true` if `x` is inside `S`.
"""->
@inline isinside(x, S::DaubSupport) = left(S) <= x <= right(S)

# Convert between values and indices of a vector with the integers in the support S
@inline x2index(x::Integer, S::DaubSupport) = x + 1 - left(S)
@inline index2x(idx::Integer, S::DaubSupport) = idx - 1 + left(S)

# Convert between values and indices of a vector with dyadic rationals
# at resolution R in the support S
@inline x2index(x, S::DaubSupport, R::Integer) = Int( (x-left(S))*2^R ) + 1
@inline index2x(idx::Integer, S::DaubSupport, R::Integer) = (idx - 1)*2.0^(-R) + left(S)


text:
export MachOSegments, MachOSegment, MachOSegmentRef

"""
    MachOSegments

`Segments` type that contains all load commands within a Mach-O file that
contain segment commands.
"""
struct MachOSegments{H <: MachOHandle, S <: MachOLoadCmdRef} <: Segments{H}
    handle::H
    segments::Vector{S}
end

function Segments(lcs::MachOLoadCmds)
    return MachOSegments(handle(lcs), findall(lcs, [MachOSegmentCmd]))
end
function Segments(oh::MachOHandle)
    # Sub out to our MachOLoadCmds-based constructor
    return Segments(MachOLoadCmds(oh))
end

handle(segments::MachOSegments) = segments.handle
lastindex(segments::MachOSegments) = lastindex(segments.segments)
function getindex(segs::MachOSegments{H}, idx) where {H <: MachOHandle}
    seg = MachOSegment(segs.segments[idx])
    return MachOSegmentRef(segs, seg, UInt32(idx))
end

"""
    MachOSegment

MachO `Segment` type which is returned from a `MachOSegmentRef` when it is
dereferenced.  Note that this type exists only for the purposes of type
conformity; all the real work is done within `MachOSegmentCmd`; but since that
already inherits from `MachOLoadCmd`, it cannot also inherit from `Segment`.
Thus, this wrapper type was born to bridge the type hierarchy.
"""
struct MachOSegment{H <: MachOHandle, S <: MachOSegmentCmd{H}} <: Segment{H}
    seg::MachOLoadCmdRef{H,S}
end

Segment(seg::MachOLoadCmdRef{H,MachOSegmentCmd{H}}) where {H <: MachOHandle} = MachOSegment(seg)

MachOLoadCmdRef(seg::MachOSegment) = seg.seg
MachOLoadCmd(seg::MachOSegment) = deref(MachOLoadCmdRef(seg))

is64bit(seg::MachOSegment) = isa(seg.seg, MachOSegment64Cmd)
segment_name(seg::MachOSegment) = segment_name(seg.seg)
segment_offset(seg::MachOSegment) = segment_offset(seg.seg)
segment_file_size(seg::MachOSegment) = segment_file_size(seg.seg)
segment_memory_size(seg::MachOSegment) = segment_memory_size(seg.seg)
segment_num_sections(seg::MachOSegment) = segment_num_sections(seg.seg)

"""
    MachOSegmentRef

MachO `SegmentRef` type which is returned from a `MachOSegments` when it is
indexed into.  Note that this type exists only for the purposes of type
conformity; all the real work is done within `MachOSegmentCmd`; but since that
already inherits from `MachOLoadCmd`, and we needed a `MachOSegment` type that
would inherit from `Segment`, we figured might as well make a `MachOSegmentRef`
as well to keep things nice and symmetric.
"""
struct MachOSegmentRef{H <: MachOHandle} <: SegmentRef{H}
    segments::MachOSegments{H}
    segment::MachOSegment{H}
    idx::UInt32
end
function SegmentRef(segs::MachOSegments, seg::MachOSegment, idx)
    return MachOSegmentRef(segs, seg, idx)
end

Segments(seg::MachOSegmentRef) = seg.segments
segment_number(seg::MachOSegmentRef) = seg.idx
deref(seg::MachOSegmentRef) = seg.segment
@derefmethod MachOLoadCmdRef(seg::MachOSegmentRef)
@derefmethod segment_num_sections(seg::MachOSegmentRef)
@derefmethod segment_name(seg::MachOSegmentRef)
@derefmethod is64bit(seg::MachOSegmentRef)

# Printing
function show(io::IO, seg::Union{MachOSegment,MachOSegmentRef})
    print(io, "MachO Segment", is64bit(seg) ? " (64 bit)" : "")

    if !get(io, :compact, false)
        fsz = string(segment_file_size(seg), base=16)
        msz = string(segment_memory_size(seg), base=16)

        println(io)
        println(io, "       Name: $(segment_name(seg))")
        println(io, "       Size: (File: 0x$(fsz),   Mem: 0x$(msz))")
        println(io, "     Offset: 0x$(string(segment_offset(seg), base=16))")
        print(io,   "   Sections: $(segment_num_sections(seg))")
    else
        print(io, " $(segment_name(seg))")
    end
end


text:
abstract type Hessian{T} end

initialize!(::Hessian) = nothing

compute_hessian!(h::AbstractMatrix, x::AbstractVector, hessian::Hessian) = hessian(h,x)

function compute_hessian(x, hessian::Hessian)
    h = alloc_h(x)
    hessian(h,x)
    return h
end

function check_hessian(H::AbstractMatrix)
    println("Condition Number of Hessian: ", cond(H))
    println("Determinant of Hessian:      ", det(H))
    println("minimum(|Hessian|):          ", minimum(abs.(H)))
    println("maximum(|Hessian|):          ", maximum(abs.(H)))
    println()
end

function print_hessian(H::AbstractMatrix)
    display(H)
    println()
end

struct HessianUser{T, HT} <: Hessian{T}
    H!::HT
end

HessianUser(H!::HT, ::AbstractVector{T}) where {T,HT} = HessianUser{T,HT}(H!)

function (hes::HessianUser{T})(H::AbstractMatrix{T}, x::AbstractVector{T}) where {T}
    hes.H!(H, x)
end

struct HessianAD{T, FT, HT <: AbstractMatrix, CT <: ForwardDiff.HessianConfig} <: Hessian{T}
    F::FT
    H::HT
    Hconfig::CT

    function HessianAD{T}(F::FT, H::HT, Hconfig::CT) where {T, FT, HT, CT}
        new{T, FT, HT, CT}(F, H, Hconfig)
    end
end

function HessianAD(F::Callable, x::AbstractVector{T}) where {T}
    Hconfig = ForwardDiff.HessianConfig(F, x)
    HessianAD{T}(F, alloc_h(x), Hconfig)
end

HessianAD(F::MultivariateObjective, x) = HessianAD(F.F, x)

HessianAD{T}(F, nx::Int) where {T} = HessianAD{T}(F, zeros(T, nx))

function (hes::HessianAD{T})(H::AbstractMatrix{T}, x::AbstractVector{T}) where {T}
    ForwardDiff.hessian!(H, hes.F, x, hes.Hconfig)
end

function (hes::HessianAD{T})(x::AbstractVector{T}) where {T}
    ForwardDiff.hessian!(hes.H, hes.F, x, hes.Hconfig)
end

function compute_hessian_ad!(H::AbstractMatrix{T}, x::AbstractVector{T}, F::FT) where {T, FT}
    hes = HessianAD(F, x)
    hes(H,x)
end

initialize!(H::HessianAD, x) = H(x)
update!(H::HessianAD, x::AbstractVector) = H(x)

Base.inv(H::HessianAD) = inv(H.H)

Base.:\(H::HessianAD, b) = H.H \ b

LinearAlgebra.ldiv!(x, H::HessianAD, b) = x .= H \ b
# LinearAlgebra.ldiv!(x, H::HessianAD, b) = LinearAlgebra.ldiv!(x, H.H, b)
# TODO: Make this work!

function Hessian(ForH, x::AbstractVector{T}; mode = :autodiff, kwargs...) where {T}
    if mode == :autodiff
        Hparams = HessianAD(ForH, x)
    else
        Hparams = HessianUser(ForH, x)
    end
    return Hparams
end

Hessian(H!, F, x::AbstractVector; kwargs...) where {T} = Hessian(H!, nx; mode = :user, kwargs...)

Hessian(H!::Nothing, F, x::AbstractVector; kwargs...) where {T} = Hessian(F, nx;  mode = :autodiff, kwargs...)

Hessian{T}(ForH, nx::Int; kwargs...) where {T} = Hessian(ForH, zeros(T,nx); kwargs...)

Hessian{T}(H!, F, nx::Int; kwargs...) where {T} = Hessian(H!, nx; kwargs...)

Hessian{T}(H!::Nothing, F, nx::Int; kwargs...) where {T} = Hessian(F, nx; kwargs...)

function compute_hessian!(h::AbstractMatrix, x::AbstractVector, ForH; kwargs...)
    hessian = Hessian(ForH, x; kwargs...)
    hessian(h,x)
end


text:
module SparseMatrixCSRTests

using SparseArrays
using SparseMatricesCSR
using LinearAlgebra
using Gridap.Algebra
using Test

maxnz=5
maxrows=5
maxcols=5
maxrowsorcols=7
int_types=(Int32,Int)
float_types=(Float32,Float64)
Bi_types=(0,1)

for Ti in int_types
  for Tv in float_types
    for Bi in Bi_types
      I = Vector{Ti}()
      J = Vector{Ti}()
      V = Vector{Tv}()
      for (ik, jk, vk) in zip(rand(1:maxrows, maxnz), rand(1:maxcols, maxnz), rand(1:Tv(maxnz), maxnz-1))
        push_coo!(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,ik,jk,vk)
      end
      push_coo!(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows,maxcols,maxnz)
      finalize_coo!(SparseMatrixCSR{Bi,Tv,Ti},I,J,V,maxrows,maxcols)
      CSC = sparse(I, J, V, maxrows,maxcols)
      CSR = sparsecsr(Val(Bi),I,J,V,maxrows,maxcols)

      @test is_entry_stored(SparseMatrixCSR{Bi,Tv,Ti},1,1)
      @test is_entry_stored(SparseMatrixCSR{Bi,Tv,Ti},1,2)
      @test is_entry_stored(SparseMatrixCSR{Bi,Tv,Ti},2,1)

      _l = 10
      _I, _J, _V = allocate_coo_vectors(SparseMatrixCSR{Bi,Tv,Ti},_l)
      @test length(_I) == _l
      @test length(_J) == _l
      @test length(_V) == _l
      @test eltype(_I) == Ti
      @test eltype(_J) == Ti
      @test eltype(_V) == Tv

      V2 = copy(V)
      V2 .= 0
      CSR2 = sparsecsr(Val(Bi),I, J, V2,maxrows,maxcols)
      copy_entries!(CSR2,CSR)
      @test CSR2 == CSR
      copy_entries!(CSR2,CSR2)
      @test CSR2 == CSR

      @test CSC == CSR
      @test nnz(CSC) == count(!iszero, CSC) == nnz(CSR) == count(!iszero, CSR)

      TCSC = sparse(J, I, V, maxrows, maxcols)
      TCSR = sparsecsr(Val(Bi), J, I, V, maxrows, maxcols)

      vold = getindex(CSR,maxrows,maxcols)
      add_entry!(+,CSR,1,maxrows,maxcols)
      @test getindex(CSR,maxrows,maxcols) == vold+1

      fill_entries!(CSR,0)
      @test all(x->x==0, nonzeros(CSR))

    end
  end
end

end # module


text:
# m6.5s.jl

using Pkg, DrWatson

using MonteCarloMeasurements
using StanSample
using StatisticalRethinking

df = CSV.read(sr_datadir("milk.csv"), DataFrame; delim=';');
scale!(df, [:kcal_per_g, :perc_fat, :perc_lactose])

stan6_5 = "
data{
    int <lower=1> N;              // Sample size
    vector[N] K;
    vector[N] F;
    vector[N] L;
}
parameters{
    real a;
    real bL;
    real bF;
    real<lower=0> sigma;
}
model{
    vector[N] mu;
    sigma ~ exponential( 1 );
    a ~ normal( 0 , 0.2 );
    bL ~ normal( 0 , 0.5 );
    bF ~ normal( 0 , 0.5 );
    mu = a + bL * L + bF * F;
    K ~ normal( mu , sigma );
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m6_5s = SampleModel("m6.5s", stan6_5);

# Input data for cmdstan

m6_5_data = Dict("N" => size(df, 1), "L" => df[:, :perc_lactose_s],
    "F" => df[:, :perc_fat_s], "K" => df[!, :kcal_per_g_s]);

# Sample using StanSample

rc6_5s = stan_sample(m6_5s, data=m6_5_data);

if success(rc6_5s)
    part6_5s = read_samples(m6_5s, :particles)
    part6_5s |> display
end

# End of m6.5s.jl


text:
#
# The roll-up combinator.
#

function RollUpQuery(qbase::Query, qks::Vector{Query})
    for qk in qks
        @assert fits(output(qbase), input(qk))
        @assert isplain(qk)
    end
    q = RecordQuery(
            qbase >>
            RecordQuery(ostub(qbase), qks))
    q = q >>
        Query(
            RollUpSig(),
            Input(domain(output(q))),
            Output([
                Output([
                    domain(output(qbase)),
                    (output(qk) |> setoptional() |> decorate(:nullrev => true) for qk in qks)...]) |> setoptional(isoptional(qbase)) |> setplural()]))
    q >>
    Query(
        GroupSig(),
        Input(domain(output(q))),
        Output([output(qbase) |> setoptional(false), (output(qk) |> setoptional(true) for qk in qks)...])
            |> setoptional(isoptional(qbase))
            |> setplural(isplural(qbase) && !isempty(qks)))

end

RollUpQuery(qbase::Query, qks::Query...) =
    RollUpQuery(qbase, collect(Query, qks))

immutable RollUpSig <: AbstractPrimitive
end

ev(sig::RollUpSig, ::Input, oty::Output, iflow::InputFlow) =
    OutputFlow(oty, ev(sig, values(iflow), fields(domain(iflow))...))

ev(::RollUpSig, dv::DataVector, oty::Output) =
    rollup_impl(column(dv, 1), fields(oty))

function rollup_impl{OPT,PLU}(dvcol::Column{OPT,PLU}, fs::Vector{Output})
    offs = offsets(dvcol)
    dv = values(dvcol)
    width = length(dv.cols)
    groupoffs = rollup_offsets_impl(width, offs)
    fs = Vector{Column}(width)
    fs[1] = rollup_values_impl(1, width, offs, values(dv, 1))
    for k = 2:width
        fs[k] = rollup_values_impl(k, width, offs, values(dv, k))
    end
    return PlainColumn(
                DataVector(
                    length(dvcol),
                    Column{OPT,true}(
                        groupoffs,
                        DataVector(width*length(dv), fs))))
end

function rollup_offsets_impl(width::Int, offs::AbstractVector{Int})
    offs′ = Vector{Int}(length(offs))
    for k = 1:endof(offs)
        offs′[k] = width * (offs[k] - 1) + 1
    end
    return offs′
end

function rollup_values_impl(k::Int, width::Int, offs::AbstractVector{Int}, vals::AbstractVector)
    len = length(vals)
    offs′ = Vector{Int}(len * width + 1)
    offs′[1] = 1
    n = 1
    m = 1
    idxs = Vector{Int}(len * (width-k+1))
    for i = 1:endof(offs)-1
        l = offs[i]
        r = offs[i+1]
        for p = 1:(width-k+1)
            for j = l:r-1
                idxs[n] = j
                n += 1
                offs′[m+1] = n
                m += 1
            end
        end
        for p = (width-k+2):width
            for j = l:r-1
                offs′[m+1] = n
                m += 1
            end
        end
    end
    return Column{k!=1,false}(offs′, vals[idxs])
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
@eval struct SiteType{T}
  (f::Type{<:SiteType})() = $(Expr(:new, :f))
end

# Note that the complicated definition of
# SiteType above is a workaround for performance
# issues when creating parameterized types
# in Julia 1.4 and 1.5-beta. Ideally we
# can just use the following in the future:
# struct SiteType{T}
# end

"""
SiteType is a parameterized type which allows
making Index tags into Julia types. Use cases
include overloading functions such as `op`,
`siteinds`, and `state` which generate custom
operators, Index arrays, and IndexVals associated
with Index objects having a certain tag.

To make a SiteType type, you can use the string
macro notation: `SiteType"MyTag"`

To make a SiteType value or object, you can use
the notation: `SiteType("MyTag")`

There are currently a few built-in site types
recognized by `ITensors.jl`. The system is easily extensible
by users. To add new operators to an existing site type,
you can follow the instructions [here](http://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_extending).
To create new site types, you can follow the instructions
[here](https://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_basic) and
[here](https://itensor.org/docs.cgi?vers=julia&page=formulas/sitetype_qns).

The current built-in site types are:

- `SiteType"S=1/2"` (or `SiteType"S=½"`)
- `SiteType"S=1"`
- `SiteType"Fermion"`
- `SiteType"tJ"`
- `SiteType"Electron"`

# Examples

Tags on indices get turned into SiteTypes internally, and then
we search for overloads of functions like `op` and `siteind`.
For example:
```julia
julia> s = siteind("S=1/2")
(dim=2|id=862|"S=1/2,Site")

julia> @show op("Sz", s);
op(s, "Sz") = ITensor ord=2
Dim 1: (dim=2|id=862|"S=1/2,Site")'
Dim 2: (dim=2|id=862|"S=1/2,Site")
NDTensors.Dense{Float64,Array{Float64,1}}
 2×2
 0.5   0.0
 0.0  -0.5

julia> @show op("Sx", s);
op(s, "Sx") = ITensor ord=2
Dim 1: (dim=2|id=862|"S=1/2,Site")'
Dim 2: (dim=2|id=862|"S=1/2,Site")
NDTensors.Dense{Float64,Array{Float64,1}}
 2×2
 0.0  0.5
 0.5  0.0

julia> @show op("Sy", s);
op(s, "Sy") = ITensor ord=2
Dim 1: (dim=2|id=862|"S=1/2,Site")'
Dim 2: (dim=2|id=862|"S=1/2,Site")
NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}}
 2×2
 0.0 + 0.0im  -0.0 - 0.5im
 0.0 + 0.5im   0.0 + 0.0im

julia> s = siteind("Electron")
(dim=4|id=734|"Electron,Site")

julia> @show op("Nup", s);
op(s, "Nup") = ITensor ord=2
Dim 1: (dim=4|id=734|"Electron,Site")'
Dim 2: (dim=4|id=734|"Electron,Site")
NDTensors.Dense{Float64,Array{Float64,1}}
 4×4
 0.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0
```

Many operators are available, for example:

- `SiteType"S=1/2"`: `"Sz"`, `"Sx"`, `"Sy"`, `"S+"`, `"S-"`, ...
- `SiteType"Electron"`: `"Nup"`, `"Ndn"`, `"Nupdn"`, `"Ntot"`, `"Cup"`, `"Cdagup"`, `"Cdn"`, `"Cdagup"`, `"Sz"`, `"Sx"`, `"Sy"`, `"S+"`, `"S-"`, ...
- ...

You can view the source code for the internal SiteType definitions
and operators that are defined [here](https://github.com/ITensor/ITensors.jl/tree/master/src/physics/site_types).
"""
SiteType(s::AbstractString) = SiteType{Tag(s)}()

SiteType(t::Tag) = SiteType{t}()

tag(::SiteType{T}) where {T} = T

macro SiteType_str(s)
  SiteType{Tag(s)}
end

# Keep TagType defined for backwards
# compatibility; will be deprecated later
const TagType = SiteType
macro TagType_str(s)
  TagType{Tag(s)}
end

#---------------------------------------
#
# op system
#
#---------------------------------------

@eval struct OpName{Name}
  (f::Type{<:OpName})() = $(Expr(:new, :f))
end

# Note that the complicated definition of
# OpName above is a workaround for performance
# issues when creating parameterized types
# in Julia 1.4 and 1.5-beta. Ideally we
# can just use the following in the future:
# struct OpName{Name}
# end

"""
OpName is a parameterized type which allows
making strings into Julia types for the purpose
of representing operator names.
The main use of OpName is overloading the
`ITensors.op!` method which generates operators
for indices with certain tags such as "S=1/2".

To make a OpName type, you can use the string
macro notation: `OpName"MyTag"`.

To make an OpName value or object, you can use
the notation: `OpName("myop")`
"""
OpName(s::AbstractString) = OpName{Symbol(s)}()
OpName(s::Symbol) = OpName{s}()
name(::OpName{N}) where {N} = N

macro OpName_str(s)
  OpName{Symbol(s)}
end

# Default implementations of op and op!
op(::OpName, ::SiteType; kwargs...) = nothing
op(::OpName, ::SiteType, ::Index...; kwargs...) = nothing
op(::OpName, ::SiteType, ::SiteType,
   sitetypes_inds::Union{SiteType, Index}...; kwargs...) = nothing
op!(::ITensor, ::OpName, ::SiteType, ::Index...; kwargs...) = nothing
op!(::ITensor, ::OpName, ::SiteType, ::SiteType,
    sitetypes_inds::Union{SiteType, Index}...; kwargs...) = nothing

# Deprecated version, for backwards compatibility
op(::SiteType, ::Index, ::AbstractString; kwargs...) = nothing

function _sitetypes(ts::TagSet)
  Ntags = length(ts)
  return SiteType[SiteType(ts[n]) for n in 1:Ntags]
end

_sitetypes(i::Index) = _sitetypes(tags(i))

"""
    op(opname::String, s::Index; kwargs...)

Return an ITensor corresponding to the operator
named `opname` for the Index `s`. The operator
is constructed by calling an overload of either
the `op` or `op!` methods which take a `SiteType`
argument that corresponds to one of the tags of
the Index `s` and an `OpName"opname"` argument
that corresponds to the input operator name.

Operator names can be combined using the `"*"`
symbol, for example `"S+*S-"` or `"Sz*Sz*Sz"`.
The result is an ITensor made by forming each operator
then contracting them together in a way corresponding
to the usual operator product or matrix multiplication.

The `op` system is used by the AutoMPO
system to convert operator names into ITensors,
and can be used directly such as for applying
operators to MPS.

# Example
```julia
s = Index(2, "Site,S=1/2")
Sz = op("Sz", s)
```
"""
function op(name::AbstractString,
            s::Index...;
            kwargs...)

  name = strip(name)

  # TODO: filter out only commons tags
  # if there are multiple indices
  commontags_s = commontags(s...)

  # Interpret operator names joined by *
  # as acting sequentially on the same site
  starpos = findfirst("*", name)
  if !isnothing(starpos)
    op1 = name[1:starpos.start-1]
    op2 = name[starpos.start+1:end]
    return product(op(op1, s...; kwargs...),
                   op(op2, s...; kwargs...))
  end

  common_stypes  = _sitetypes(commontags_s)
  push!(common_stypes,SiteType("Generic"))
  opn = OpName(name)

  #
  # Try calling a function of the form:
  #    op(::OpName, ::SiteType, ::Index; kwargs...)
  #
  for st in common_stypes
    res = op(opn, st, s...; kwargs...)
    if !isnothing(res)
      return res
    end
  end

  # otherwise try calling a function of the form:
  #    op!(::ITensor, ::OpName, ::SiteType, ::Index; kwargs...)
  #
  Op = emptyITensor(prime.(s)..., dag.(s)...)
  for st in common_stypes
    op!(Op, opn, st, s...; kwargs...)
    if !isempty(Op)
      return Op
    end
  end

  #
  # otherwise try calling a function of the form:
  #    op(::OpName, ::SiteType; kwargs...)
  # which returns a Julia matrix
  #
  for st in common_stypes
    op_mat = op(opn, st; kwargs...)
    if !isnothing(op_mat)
      rs = reverse(s)
      return itensor(op_mat, prime.(rs)..., dag.(rs)...)
    end
  end

  if length(s) > 1
    # No overloads for common tags found. It might be a
    # case of making an operator with mixed site types,
    # searching for overloads like:
    #   op(::OpName,
    #      ::SiteType...,
    #      ::Index...;
    #      kwargs...)
    #   op!(::ITensor, ::OpName,
    #       ::SiteType...,
    #       ::Index...;
    #       kwargs...)
    stypes = _sitetypes.(s)

    for st in Iterators.product(stypes...)
      res = op(opn, st..., s...; kwargs...)
      if !isnothing(res)
        return res
      end
    end

    Op = emptyITensor(prime.(s)..., dag.(s)...)
    for st in Iterators.product(stypes...)
      op!(Op, opn, st..., s...; kwargs...)
      if !isempty(Op)
        return Op
      end
    end
    error("Older op interface does not support multiple indices with mixed site types. You may want to overload `op(::OpName, ::SiteType..., ::Index...)` or `op!(::ITensor, ::OpName, ::SiteType..., ::Index...) for the operator \"$name\" and Index tags $(tags.(s)).")
  end

  #
  # otherwise try calling a function of the form:
  #   op(::SiteType, ::Index, ::AbstractString)
  #
  # (Note: this version is for backwards compatibility
  #  after version 0.1.10, and may be eventually
  #  deprecated)
  #
  for st in common_stypes
    res = op(st, s[1], name; kwargs...)
    if !isnothing(res)
      return res
    end
  end

  throw(ArgumentError("Overload of \"op\" or \"op!\" functions not found for operator name \"$name\" and Index tags: $(commontags_s))"))
end

# For backwards compatibility, version of `op`
# taking the arguments in the other order:
op(s::Index,
   opname::AbstractString;
   kwargs...) = op(opname, s; kwargs...)

# To ease calling of other op overloads,
# allow passing a string as the op name
op(opname::AbstractString,t::SiteType) = op(OpName(opname),t)

"""
    op(opname::String,sites::Vector{<:Index},n::Int; kwargs...)

Return an ITensor corresponding to the operator
named `opname` for the n'th Index in the array
`sites`.

# Example

```julia
s = siteinds("S=1/2", 4)
Sz2 = op("Sz", s, 2)
```
"""
op(opname::AbstractString, s::Vector{<:Index},
   ns::Vararg{Int, N}; kwargs...) where {N} =
  op(opname, ntuple(n -> s[ns[n]], Val(N))...; kwargs...)

op(s::Vector{ <: Index}, opname::AbstractString,
   ns::Int...; kwargs...) =
  op(opname, s, ns...; kwargs...)

op(s::Vector{ <: Index}, opname::AbstractString,
   ns::Tuple{Vararg{Int}}, kwargs::NamedTuple) =
  op(opname, s, ns...; kwargs...)

op(s::Vector{ <: Index}, opname::AbstractString,
   ns::Int, kwargs::NamedTuple) =
  op(opname, s, ns; kwargs...)

# This version helps with call like `op.(Ref(s), os)` where `os`
# is a vector of tuples.
op(s::Vector{ <: Index}, os::Tuple{String, Vararg}) =
  op(s, os...)

# Here, Ref is used to not broadcast over the vector of indices
# TODO: consider overloading broadcast for `op` with the example
# here: https://discourse.julialang.org/t/how-to-broadcast-over-only-certain-function-arguments/19274/5
# so that `Ref` isn't needed.
ops(s::Vector{ <: Index},
    os::AbstractArray{ <: Tuple{String, Vararg}}) =
  op.(Ref(s), os)

ops(os::Vector{ <: Tuple{String, Vararg}}, s::Vector{ <: Index}) =
  op.(Ref(s), os)

#---------------------------------------
#
# state system
#
#---------------------------------------

@eval struct StateName{Name}
  (f::Type{<:StateName})() = $(Expr(:new, :f))
end

StateName(s::AbstractString) = StateName{SmallString(s)}()
StateName(s::SmallString) = StateName{s}()
name(::StateName{N}) where {N} = N

macro StateName_str(s)
  StateName{SmallString(s)}
end

state(::SiteType, ::StateName) = nothing
state(::SiteType, ::AbstractString) = nothing

function state(s::Index,
               name::AbstractString)::IndexVal
  stypes  = _sitetypes(s)
  sname = StateName(name)

  # Try calling state(::SiteType"Tag",::StateName"Name")
  for st in stypes
    res = state(st,sname)
    !isnothing(res) && return s(res)
  end

  # Try calling state(::SiteType"Tag","Name")
  for st in stypes
    res = state(st,name)
    !isnothing(res) && return s(res)
  end

  throw(ArgumentError("Overload of \"state\" function not found for Index tags $(tags(s))"))
end

state(s::Index,n::Integer) = s[n]

state(sset::Vector{<:Index},j::Integer,st) = state(sset[j],st)

#---------------------------------------
#
# siteind system
#
#---------------------------------------

space(st::SiteType; kwargs...) = nothing

space(st::SiteType, n::Int; kwargs...) =
  space(st; kwargs...)

function space_error_message(st::SiteType)
  return "Overload of \"space\",\"siteind\", or \"siteinds\" functions not found for Index tag: $(tag(st))"
end

function siteind(st::SiteType; addtags = "", kwargs...)
  sp = space(st; kwargs...)
  isnothing(sp) && return nothing
  return Index(sp, "Site, $(tag(st)), $addtags")
end

function siteind(st::SiteType, n; kwargs...)
  s = siteind(st; kwargs...)
  !isnothing(s) && return addtags(s, "n=$n")
  sp = space(st, n; kwargs...)
  isnothing(sp) && error(space_error_message(st))
  return Index(sp, "Site, $(tag(st)), n=$n")
end

siteind(tag::String; kwargs...) =
  siteind(SiteType(tag); kwargs...)

siteind(tag::String, n; kwargs...) =
  siteind(SiteType(tag), n; kwargs...)

# Special case of `siteind` where integer (dim) provided
# instead of a tag string
siteind(d::Integer,n::Integer; kwargs...) = Index(d,"Site,n=$n")

#---------------------------------------
#
# siteinds system
#
#---------------------------------------

siteinds(::SiteType, N; kwargs...) = nothing

"""
  siteinds(tag::String, N::Integer; kwargs...)

Create an array of `N` physical site indices of type `tag`.
Keyword arguments can be used to specify quantum number conservation,
see the `space` function corresponding to the site type `tag` for
supported keyword arguments.
"""
function siteinds(tag::String,
                  N::Integer; kwargs...)
  st = SiteType(tag)

  si = siteinds(st,N;kwargs...)
  if !isnothing(si)
    return si
  end

  return [siteind(st,j; kwargs...) for j=1:N]
end

"""
  siteinds(f::Function, N::Integer; kwargs...)

Create an array of `N` physical site indices where the site type at site `n` is given
by `f(n)` (`f` should return a string).
"""
function siteinds(f::Function,
                  N::Integer; kwargs...)
  [siteind(f(n),n; kwargs...) for n=1:N]
end

# Special case of `siteinds` where integer (dim)
# provided instead of a tag string
function siteinds(d::Integer,
                  N::Integer; kwargs...)
  return [siteind(d,n; kwargs...) for n=1:N]
end

#---------------------------------------
#
# has_fermion_string system
#
#---------------------------------------

has_fermion_string(::OpName, ::SiteType) = nothing

function has_fermion_string(opname::AbstractString,
                            s::Index;
                            kwargs...)::Bool
  opname = strip(opname)

  # Interpret operator names joined by *
  # as acting sequentially on the same site
  starpos = findfirst("*", opname)
  if !isnothing(starpos)
    op1 = opname[1:starpos.start-1]
    op2 = opname[starpos.start+1:end]
    return xor(has_fermion_string(op1,s; kwargs...),
               has_fermion_string(op2,s; kwargs...))
  end

  Ntags = length(tags(s))
  stypes = _sitetypes(s)
  opn = OpName(opname)
  for st in stypes
    res = has_fermion_string(opn, st)
    !isnothing(res) && return res
  end
  return false
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of row number column.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Show row number" begin
    expected = """
┌─────┬────────┬────────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├─────┼────────┼────────┼────────┼────────┤
│   1 │ 1      │  false │  1.0   │      1 │
│   2 │ 2      │   true │  2.0   │      2 │
│   3 │ 3      │  false │  3.0   │      3 │
│   4 │ 4      │   true │  4.0   │      4 │
│   5 │ 5      │  false │  5.0   │      5 │
│   6 │ 6      │   true │  6.0   │      6 │
└─────┴────────┴────────┴────────┴────────┘
"""
    result = pretty_table(String, data;
                          alignment       = [:l,:r,:c,:r],
                          show_row_number = true)
    @test result == expected

    expected = """
┌─────┬────────┬────────┬────────┬────────┐
│ Row │ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├─────┼────────┼────────┼────────┼────────┤
│  1  │ 1      │  false │  1.0   │      1 │
│  2  │ 2      │   true │  2.0   │      2 │
│  3  │ 3      │  false │  3.0   │      3 │
│  4  │ 4      │   true │  4.0   │      4 │
│  5  │ 5      │  false │  5.0   │      5 │
│  6  │ 6      │   true │  6.0   │      6 │
└─────┴────────┴────────┴────────┴────────┘
"""

    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          row_number_alignment = :c,
                          show_row_number = true)
    @test result == expected

    expected = """
┌───┬────────┬────────┬────────┬────────┐
│ # │ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├───┼────────┼────────┼────────┼────────┤
│ 1 │ 1      │  false │  1.0   │      1 │
│ 2 │ 2      │   true │  2.0   │      2 │
│ 3 │ 3      │  false │  3.0   │      3 │
│ 4 │ 4      │   true │  4.0   │      4 │
│ 5 │ 5      │  false │  5.0   │      5 │
│ 6 │ 6      │   true │  6.0   │      6 │
└───┴────────┴────────┴────────┴────────┘
"""

    result = pretty_table(String, data;
                          alignment = [:l,:r,:c,:r],
                          row_number_alignment = :c,
                          row_number_column_title = "#",
                          show_row_number = true)
    @test result == expected
end